1// The linearized plate as matrices: the stiffness K (the Hessian of the
2// discrete energy of reference.ts at the stress-free state) in symmetric
3// banded storage, and the lumped mass M as a diagonal. A flat plate has one
4// unknown per node (the transverse displacement w; the in-plane motion
5// decouples and is neither struck nor heard). A plate with a profile couples
6// w to the in-plane displacements u and v through the slopes of the profile,
7// and has three unknowns per node, interleaved so that the band stays narrow.
9import {
10 type Plate,
11 bendingStiffness,
12 cellThickness,
13 insideCell,
14 insideNode,
15 membraneStiffness,
16 nodeMass,
17 nodeWeight,
18 shearStiffness,
19} from "../params.ts";
20import { shapeDerivatives } from "./reference.ts";
22export interface LinearSystem {
23 /** Number of unknowns. */
24 n: number;
25 /** Half-bandwidth: K[i, j] = 0 for |i - j| > bw. */
26 bw: number;
27 /** Lower band of K, row-major: band[i * (bw + 1) + d] = K[i, i - d] for d = 0..bw. */
28 band: Float64Array;
29 /** Diagonal mass of every unknown [kg]. */
30 mass: Float64Array;
31 /** 1 for a flat plate (w only), 3 for a plate with a profile (u, v, w interleaved). */
32 fields: 1 | 3;
33 /** Grid node of unknown i (its node index j * nx + i on the grid). */
34 dofNode: Int32Array;
35 /** First unknown of a grid node, or -1 outside the plate. */
36 nodeDof: Int32Array;
37}
39/** Whether the profile vanishes, in which case only the transverse motion is assembled. */
40export function isFlatPlate(pl: Plate): boolean {
41 for (let k = 0; k < pl.w0.length; k++) if (pl.w0[k] !== 0) return false;
42 return true;
43}
45/**
46 * Assembles K and M for the plate. `fields` may force the three-field
47 * system for a flat plate (the tests use it); by default a flat plate gets
48 * the one-field system.
49 */
50export function assemble(pl: Plate, fields?: 1 | 3): LinearSystem {
51 const g = pl.grid, m = pl.mat;
52 const { nx, ny, hx, hy } = g;
53 const nf: 1 | 3 = fields ?? (isFlatPlate(pl) ? 1 : 3);
55 // Number the plate nodes with the shorter grid axis innermost so that the
56 // 13-point bending stencil, which reaches two rows away, gives the
57 // narrowest band.
58 const rowMajor = nx <= ny;
59 const nodeDof = new Int32Array(g.n).fill(-1);
60 const nodeList: number[] = [];
61 const outer = rowMajor ? ny : nx, inner = rowMajor ? nx : ny;
62 for (let o = 0; o < outer; o++) {
63 for (let q = 0; q < inner; q++) {
64 const i = rowMajor ? q : o, j = rowMajor ? o : q;
65 if (!insideNode(g, i, j) || nodeWeight(g, i, j) === 0) continue;
66 nodeDof[j * nx + i] = nodeList.length * nf;
67 nodeList.push(j * nx + i);
68 }
69 }
70 const nNodes = nodeList.length;
71 const n = nNodes * nf;
72 const dofNode = new Int32Array(n);
73 for (let p = 0; p < nNodes; p++) for (let c = 0; c < nf; c++) dofNode[p * nf + c] = nodeList[p];
74 const bw = nf * (2 * inner) + (nf - 1);
75 const band = new Float64Array(n * (bw + 1));
76 const mass = new Float64Array(n);
77 for (let p = 0; p < nNodes; p++) {
78 const k = nodeList[p];
79 const mk = nodeMass(pl, k % nx, Math.floor(k / nx));
80 for (let c = 0; c < nf; c++) mass[p * nf + c] = mk;
81 }
83 const add = (r: number, c: number, v: number): void => {
84 if (v === 0) return;
85 if (r < c) { const t = r; r = c; c = t; }
86 const d = r - c;
87 if (d > bw) throw new Error(`assembly outside the band (${d} > ${bw})`);
88 band[r * (bw + 1) + d] += v;
89 };
90 const wdof = (k: number) => nodeDof[k] + (nf - 1); // the w unknown of grid node k
91 const id = (i: number, j: number) => j * nx + i;
92 const nu = m.nu;
94 // ---- Bending: hx hy [ sum_k w_k D_k (a a' + b b' + nu (a b' + b a')) + 2 (1 - nu) sum_c D_c c c' ] ----
95 const stencilA = (i: number, j: number): [number[], number[]] | null => {
96 if (!(insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j))) return null;
97 return [[wdof(id(i - 1, j)), wdof(id(i, j)), wdof(id(i + 1, j))], [1 / (hx * hx), -2 / (hx * hx), 1 / (hx * hx)]];
98 };
99 const stencilB = (i: number, j: number): [number[], number[]] | null => {
100 if (!(insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1))) return null;
101 return [[wdof(id(i, j - 1)), wdof(id(i, j)), wdof(id(i, j + 1))], [1 / (hy * hy), -2 / (hy * hy), 1 / (hy * hy)]];
102 };
103 const outer2 = (s: [number[], number[]], t: [number[], number[]], f: number): void => {
104 for (let a = 0; a < s[0].length; a++) for (let b = 0; b < t[0].length; b++) {
105 const r = s[0][a], c = t[0][b];
106 // Add the full outer product once per (r, c) pair. The band stores
107 // only r >= c, and the symmetric counterpart is added when the loop
108 // reaches it, so skip r < c here unless the product is symmetric in
109 // (s, t), in which case each unordered pair appears twice.
110 if (r >= c) add(r, c, f * s[1][a] * t[1][b]);
111 }
112 };
113 for (let j = 0; j < ny; j++) {
114 for (let i = 0; i < nx; i++) {
115 const wgt = nodeWeight(g, i, j);
116 if (wgt === 0) continue;
117 const Dk = bendingStiffness(m, pl.h[id(i, j)]);
118 const f = hx * hy * wgt * Dk;
119 const A = stencilA(i, j), B = stencilB(i, j);
120 if (A) outer2(A, A, f);
121 if (B) outer2(B, B, f);
122 if (A && B) { outer2(A, B, f * nu); outer2(B, A, f * nu); }
123 }
124 }
125 for (let cj = 0; cj < ny - 1; cj++) {
126 for (let ci = 0; ci < nx - 1; ci++) {
127 if (!insideCell(g, ci, cj)) continue;
128 const Dc = bendingStiffness(m, cellThickness(pl, ci, cj));
129 const C: [number[], number[]] = [
130 [wdof(id(ci, cj)), wdof(id(ci + 1, cj)), wdof(id(ci, cj + 1)), wdof(id(ci + 1, cj + 1))],
131 [1 / (hx * hy), -1 / (hx * hy), -1 / (hx * hy), 1 / (hx * hy)],
132 ];
133 outer2(C, C, 2 * (1 - nu) * Dc * hx * hy);
134 }
135 }
137 // ---- Membrane, linearized about the profile (three-field system only) ----
138 if (nf === 3) {
139 const q = (hx * hy) / 4;
140 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
141 const Bm = new Float64Array(3 * 12);
142 const dofs = new Int32Array(12);
143 const Ke = new Float64Array(12 * 12);
144 for (let cj = 0; cj < ny - 1; cj++) {
145 for (let ci = 0; ci < nx - 1; ci++) {
146 if (!insideCell(g, ci, cj)) continue;
147 const hc = cellThickness(pl, ci, cj);
148 const Cm = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
149 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
150 for (let b = 0; b < 4; b++) {
151 dofs[b] = nodeDof[nn[b]]; dofs[4 + b] = nodeDof[nn[b]] + 1; dofs[8 + b] = nodeDof[nn[b]] + 2;
152 }
153 Ke.fill(0);
154 for (let gp = 0; gp < 4; gp++) {
155 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
156 let w0x = 0, w0y = 0;
157 for (let b = 0; b < 4; b++) { w0x += dNdx[b] * pl.w0[nn[b]]; w0y += dNdy[b] * pl.w0[nn[b]]; }
158 // Linear strains: exx = u_x + w0_x w_x, eyy = v_y + w0_y w_y, gxy = u_y + v_x + w0_x w_y + w0_y w_x.
159 Bm.fill(0);
160 for (let b = 0; b < 4; b++) {
161 Bm[0 * 12 + b] = dNdx[b]; Bm[0 * 12 + 8 + b] = w0x * dNdx[b];
162 Bm[1 * 12 + 4 + b] = dNdy[b]; Bm[1 * 12 + 8 + b] = w0y * dNdy[b];
163 Bm[2 * 12 + b] = dNdy[b]; Bm[2 * 12 + 4 + b] = dNdx[b]; Bm[2 * 12 + 8 + b] = w0x * dNdy[b] + w0y * dNdx[b];
164 }
165 // Ke += q B' C B with C = [[Cm, nu Cm, 0], [nu Cm, Cm, 0], [0, 0, Gs]].
166 for (let a = 0; a < 12; a++) {
167 const b0 = Bm[a], b1 = Bm[12 + a], b2 = Bm[24 + a];
168 const s0 = Cm * (b0 + nu * b1), s1 = Cm * (b1 + nu * b0), s2 = Gs * b2;
169 for (let c = 0; c <= a; c++) Ke[a * 12 + c] += q * (s0 * Bm[c] + s1 * Bm[12 + c] + s2 * Bm[24 + c]);
170 }
171 }
172 for (let a = 0; a < 12; a++) for (let c = 0; c <= a; c++) {
173 if (a === c) add(dofs[a], dofs[c], Ke[a * 12 + c]);
174 else if (dofs[a] !== dofs[c]) add(dofs[a], dofs[c], Ke[a * 12 + c]);
175 }
176 }
177 }
178 }
180 return { n, bw, band, mass, fields: nf, dofNode, nodeDof };
181}
183/** y = K x from the banded storage. */
184export function stiffnessMatvec(sys: LinearSystem, x: Float64Array, y: Float64Array): void {
185 const { n, bw, band } = sys;
186 y.fill(0);
187 for (let i = 0; i < n; i++) {
188 const base = i * (bw + 1);
189 const dmax = Math.min(bw, i);
190 let s = band[base] * x[i];
191 for (let d = 1; d <= dmax; d++) {
192 const v = band[base + d];
193 if (v === 0) continue;
194 s += v * x[i - d];
195 y[i - d] += v * x[i];
196 }
197 y[i] += s;
198 }
199}
201/** K[i, j] from the banded storage (zero outside the band). */
202export function stiffnessEntry(sys: LinearSystem, i: number, j: number): number {
203 if (i < j) { const t = i; i = j; j = t; }
204 const d = i - j;
205 return d > sys.bw ? 0 : sys.band[i * (sys.bw + 1) + d];
206}
208/** Scatters an unknown vector to grid arrays (u, v, w), zero outside the plate; u and v are zero for the one-field system. */
209export function toGrid(sys: LinearSystem, x: ArrayLike<number>, nGrid: number): { u: Float64Array; v: Float64Array; w: Float64Array } {
210 const u = new Float64Array(nGrid), v = new Float64Array(nGrid), w = new Float64Array(nGrid);
211 for (let i = 0; i < sys.n; i++) {
212 const k = sys.dofNode[i];
213 if (sys.fields === 1) w[k] = x[i];
214 else {
215 const c = i % 3;
216 if (c === 0) u[k] = x[i]; else if (c === 1) v[k] = x[i]; else w[k] = x[i];
217 }
218 }
219 return { u, v, w };
220}
222/** Gathers grid arrays into an unknown vector. */
223export function fromGrid(sys: LinearSystem, f: { u: Float64Array; v: Float64Array; w: Float64Array }): Float64Array {
224 const x = new Float64Array(sys.n);
225 for (let i = 0; i < sys.n; i++) {
226 const k = sys.dofNode[i];
227 if (sys.fields === 1) x[i] = f.w[k];
228 else {
229 const c = i % 3;
230 x[i] = c === 0 ? f.u[k] : c === 1 ? f.v[k] : f.w[k];
231 }
232 }
233 return x;
234}