5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 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 type StandRates,
12 bendingStiffness,
13 cellThickness,
14 insideCell,
15 insideNode,
16 membraneStiffness,
17 nodeMass,
18 nodeWeight,
19 shearStiffness,
20 standRates,
21} from "../params.ts";
22import { shapeDerivatives } from "./reference.ts";
24export interface LinearSystem {
25 /** Number of unknowns. */
26 n: number;
27 /** Half-bandwidth: K[i, j] = 0 for |i - j| > bw. */
28 bw: number;
29 /** Lower band of K, row-major: band[i * (bw + 1) + d] = K[i, i - d] for d = 0..bw. */
30 band: Float64Array;
31 /** Diagonal mass of every unknown [kg]. */
32 mass: Float64Array;
33 /** 1 for a flat plate (w only), 3 for a plate with a profile (u, v, w interleaved). */
34 fields: 1 | 3;
35 /** Grid node of unknown i (its node index j * nx + i on the grid). */
36 dofNode: Int32Array;
37 /** First unknown of a grid node, or -1 outside the plate. */
38 nodeDof: Int32Array;
39 /** Whether the plate is on a stand (no rigid-body modes at zero frequency). */
40 onStand: boolean;
41}
43/** Whether the profile vanishes, in which case only the transverse motion is assembled. */
44export function isFlatPlate(pl: Plate): boolean {
45 for (let k = 0; k < pl.w0.length; k++) if (pl.w0[k] !== 0) return false;
46 return true;
47}
49/**
50 * Assembles K and M for the plate. `fields` may force the three-field
51 * system for a flat plate (the tests use it); by default a flat plate gets
52 * the one-field system.
53 */
54export function assemble(pl: Plate, fields?: 1 | 3, sr: StandRates | null = standRates(pl)): LinearSystem {
55 const g = pl.grid, m = pl.mat;
56 const { nx, ny, hx, hy } = g;
57 const nf: 1 | 3 = fields ?? (isFlatPlate(pl) ? 1 : 3);
59 // Number the plate nodes with the shorter grid axis innermost so that the
60 // 13-point bending stencil, which reaches two rows away, gives the
61 // narrowest band.
62 const rowMajor = nx <= ny;
63 const nodeDof = new Int32Array(g.n).fill(-1);
64 const nodeList: number[] = [];
65 const outer = rowMajor ? ny : nx, inner = rowMajor ? nx : ny;
66 for (let o = 0; o < outer; o++) {
67 for (let q = 0; q < inner; q++) {
68 const i = rowMajor ? q : o, j = rowMajor ? o : q;
69 if (!insideNode(g, i, j) || nodeWeight(g, i, j) === 0) continue;
70 nodeDof[j * nx + i] = nodeList.length * nf;
71 nodeList.push(j * nx + i);
72 }
73 }
74 const nNodes = nodeList.length;
75 const n = nNodes * nf;
76 const dofNode = new Int32Array(n);
77 for (let p = 0; p < nNodes; p++) for (let c = 0; c < nf; c++) dofNode[p * nf + c] = nodeList[p];
78 const bw = nf * (2 * inner) + (nf - 1);
79 const band = new Float64Array(n * (bw + 1));
80 const mass = new Float64Array(n);
81 for (let p = 0; p < nNodes; p++) {
82 const k = nodeList[p];
83 const mk = nodeMass(pl, k % nx, Math.floor(k / nx));
84 for (let c = 0; c < nf; c++) mass[p * nf + c] = mk;
85 }
87 const add = (r: number, c: number, v: number): void => {
88 if (v === 0) return;
89 if (r < c) { const t = r; r = c; c = t; }
90 const d = r - c;
91 if (d > bw) throw new Error(`assembly outside the band (${d} > ${bw})`);
92 band[r * (bw + 1) + d] += v;
93 };
94 const wdof = (k: number) => nodeDof[k] + (nf - 1); // the w unknown of grid node k
95 const id = (i: number, j: number) => j * nx + i;
96 const nu = m.nu;
98 // ---- 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' ] ----
99 const stencilA = (i: number, j: number): [number[], number[]] | null => {
100 if (!(insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j))) return null;
101 return [[wdof(id(i - 1, j)), wdof(id(i, j)), wdof(id(i + 1, j))], [1 / (hx * hx), -2 / (hx * hx), 1 / (hx * hx)]];
102 };
103 const stencilB = (i: number, j: number): [number[], number[]] | null => {
104 if (!(insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1))) return null;
105 return [[wdof(id(i, j - 1)), wdof(id(i, j)), wdof(id(i, j + 1))], [1 / (hy * hy), -2 / (hy * hy), 1 / (hy * hy)]];
106 };
107 const outer2 = (s: [number[], number[]], t: [number[], number[]], f: number): void => {
108 for (let a = 0; a < s[0].length; a++) for (let b = 0; b < t[0].length; b++) {
109 const r = s[0][a], c = t[0][b];
110 // Add the full outer product once per (r, c) pair. The band stores
111 // only r >= c, and the symmetric counterpart is added when the loop
112 // reaches it, so skip r < c here unless the product is symmetric in
113 // (s, t), in which case each unordered pair appears twice.
114 if (r >= c) add(r, c, f * s[1][a] * t[1][b]);
115 }
116 };
117 for (let j = 0; j < ny; j++) {
118 for (let i = 0; i < nx; i++) {
119 const wgt = nodeWeight(g, i, j);
120 if (wgt === 0) continue;
121 const Dk = bendingStiffness(m, pl.h[id(i, j)]);
122 const f = hx * hy * wgt * Dk;
123 const A = stencilA(i, j), B = stencilB(i, j);
124 if (A) outer2(A, A, f);
125 if (B) outer2(B, B, f);
126 if (A && B) { outer2(A, B, f * nu); outer2(B, A, f * nu); }
127 }
128 }
129 for (let cj = 0; cj < ny - 1; cj++) {
130 for (let ci = 0; ci < nx - 1; ci++) {
131 if (!insideCell(g, ci, cj)) continue;
132 const Dc = bendingStiffness(m, cellThickness(pl, ci, cj));
133 const C: [number[], number[]] = [
134 [wdof(id(ci, cj)), wdof(id(ci + 1, cj)), wdof(id(ci, cj + 1)), wdof(id(ci + 1, cj + 1))],
135 [1 / (hx * hy), -1 / (hx * hy), -1 / (hx * hy), 1 / (hx * hy)],
136 ];
137 outer2(C, C, 2 * (1 - nu) * Dc * hx * hy);
138 }
139 }
141 // ---- Membrane, linearized about the profile (three-field system only) ----
142 if (nf === 3) {
143 const q = (hx * hy) / 4;
144 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
145 const Bm = new Float64Array(3 * 12);
146 const dofs = new Int32Array(12);
147 const Ke = new Float64Array(12 * 12);
148 for (let cj = 0; cj < ny - 1; cj++) {
149 for (let ci = 0; ci < nx - 1; ci++) {
150 if (!insideCell(g, ci, cj)) continue;
151 const hc = cellThickness(pl, ci, cj);
152 const Cm = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
153 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
154 for (let b = 0; b < 4; b++) {
155 dofs[b] = nodeDof[nn[b]]; dofs[4 + b] = nodeDof[nn[b]] + 1; dofs[8 + b] = nodeDof[nn[b]] + 2;
156 }
157 Ke.fill(0);
158 for (let gp = 0; gp < 4; gp++) {
159 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
160 let w0x = 0, w0y = 0;
161 for (let b = 0; b < 4; b++) { w0x += dNdx[b] * pl.w0[nn[b]]; w0y += dNdy[b] * pl.w0[nn[b]]; }
162 // 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.
163 Bm.fill(0);
164 for (let b = 0; b < 4; b++) {
165 Bm[0 * 12 + b] = dNdx[b]; Bm[0 * 12 + 8 + b] = w0x * dNdx[b];
166 Bm[1 * 12 + 4 + b] = dNdy[b]; Bm[1 * 12 + 8 + b] = w0y * dNdy[b];
167 Bm[2 * 12 + b] = dNdy[b]; Bm[2 * 12 + 4 + b] = dNdx[b]; Bm[2 * 12 + 8 + b] = w0x * dNdy[b] + w0y * dNdx[b];
168 }
169 // Ke += q B' C B with C = [[Cm, nu Cm, 0], [nu Cm, Cm, 0], [0, 0, Gs]].
170 for (let a = 0; a < 12; a++) {
171 const b0 = Bm[a], b1 = Bm[12 + a], b2 = Bm[24 + a];
172 const s0 = Cm * (b0 + nu * b1), s1 = Cm * (b1 + nu * b0), s2 = Gs * b2;
173 for (let c = 0; c <= a; c++) Ke[a * 12 + c] += q * (s0 * Bm[c] + s1 * Bm[12 + c] + s2 * Bm[24 + c]);
174 }
175 }
176 for (let a = 0; a < 12; a++) for (let c = 0; c <= a; c++) {
177 if (a === c) add(dofs[a], dofs[c], Ke[a * 12 + c]);
178 else if (dofs[a] !== dofs[c]) add(dofs[a], dofs[c], Ke[a * 12 + c]);
179 }
180 }
181 }
182 }
184 // ---- Stand: the pads as springs to ground ----
185 if (sr) {
186 for (const k of sr.nodes) {
187 if (nodeDof[k] < 0) continue;
188 add(wdof(k), wdof(k), sr.kw[k]);
189 if (nf === 3) { add(nodeDof[k], nodeDof[k], sr.ku[k]); add(nodeDof[k] + 1, nodeDof[k] + 1, sr.ku[k]); }
190 }
191 }
193 return { n, bw, band, mass, fields: nf, dofNode, nodeDof, onStand: !!sr && sr.nodes.length > 0 };
194}
196/** y = K x from the banded storage. */
197export function stiffnessMatvec(sys: LinearSystem, x: Float64Array, y: Float64Array): void {
198 const { n, bw, band } = sys;
199 y.fill(0);
200 for (let i = 0; i < n; i++) {
201 const base = i * (bw + 1);
202 const dmax = Math.min(bw, i);
203 let s = band[base] * x[i];
204 for (let d = 1; d <= dmax; d++) {
205 const v = band[base + d];
206 if (v === 0) continue;
207 s += v * x[i - d];
208 y[i - d] += v * x[i];
209 }
210 y[i] += s;
211 }
212}
214/** K[i, j] from the banded storage (zero outside the band). */
215export function stiffnessEntry(sys: LinearSystem, i: number, j: number): number {
216 if (i < j) { const t = i; i = j; j = t; }
217 const d = i - j;
218 return d > sys.bw ? 0 : sys.band[i * (sys.bw + 1) + d];
219}
221/** Scatters an unknown vector to grid arrays (u, v, w), zero outside the plate; u and v are zero for the one-field system. */
222export function toGrid(sys: LinearSystem, x: ArrayLike<number>, nGrid: number): { u: Float64Array; v: Float64Array; w: Float64Array } {
223 const u = new Float64Array(nGrid), v = new Float64Array(nGrid), w = new Float64Array(nGrid);
224 for (let i = 0; i < sys.n; i++) {
225 const k = sys.dofNode[i];
226 if (sys.fields === 1) w[k] = x[i];
227 else {
228 const c = i % 3;
229 if (c === 0) u[k] = x[i]; else if (c === 1) v[k] = x[i]; else w[k] = x[i];
230 }
231 }
232 return { u, v, w };
233}
235/** Gathers grid arrays into an unknown vector. */
236export function fromGrid(sys: LinearSystem, f: { u: Float64Array; v: Float64Array; w: Float64Array }): Float64Array {
237 const x = new Float64Array(sys.n);
238 for (let i = 0; i < sys.n; i++) {
239 const k = sys.dofNode[i];
240 if (sys.fields === 1) x[i] = f.w[k];
241 else {
242 const c = i % 3;
243 x[i] = c === 0 ? f.u[k] : c === 1 ? f.v[k] : f.w[k];
244 }
245 }
246 return x;
247}