1ee087eCymbal designer: scripted design, thickness field, mount, single strikeJeremy Magland 1// Float64 reference implementation of the spatial discretization. It mirrors
2// step.wgsl and is used for the time-step estimate, the tests, and the
3// GPU-versus-CPU self-test. It is not used in the time loop.
5import {
6 type Grid,
7 type Material,
8 type Patch,
9 type Plate,
10 bendingStiffness,
11 cellThickness,
12 insideCell,
13 insideNode,
14 membraneStiffness,
15 mountPatch,
16 nodeMass,
17 nodeWeight,
18 shearStiffness,
19} from "../params.ts";
21export interface Fields {
22 u: Float64Array;
23 v: Float64Array;
24 w: Float64Array;
25}
27export function fields(n: number): Fields {
28 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
29}
31/** A plate of uniform thickness on the given grid (for the tests and the time-step estimate). */
32export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
33 const hs = new Float64Array(g.n);
34 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (insideNode(g, i, j)) hs[j * g.nx + i] = h;
35 return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), mount: null, loss: { sigma0: 0, sigma1: 0, sigmaIn: 0 } };
36}
38const G3 = 1 / Math.sqrt(3);
40/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
41function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
42 const xi = gp & 1 ? G3 : -G3;
43 const eta = gp & 2 ? G3 : -G3;
44 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
45 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
46 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
47 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
48 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
49 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
50 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
51 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
52}
54/** Per-node spring stiffness of the mount [N/m] (zero away from the mount). */
55export function mountStiffness(pl: Plate, patch: Patch | null): Float64Array {
56 const g = pl.grid;
57 const K = new Float64Array(g.n);
58 if (!pl.mount || !patch || patch.weight === 0) return K;
59 for (const k of patch.nodes) K[k] = (pl.mount.stiffness * nodeWeight(g, k % g.nx, Math.floor(k / g.nx))) / patch.weight;
60 return K;
61}
63/** Per-node dashpot rate of the mount [N s/m] (zero away from the mount). */
64export function mountDamping(pl: Plate, patch: Patch | null): Float64Array {
65 const g = pl.grid;
66 const C = new Float64Array(g.n);
67 if (!pl.mount || !patch || patch.weight === 0) return C;
68 for (const k of patch.nodes) C[k] = (pl.mount.damping * nodeWeight(g, k % g.nx, Math.floor(k / g.nx))) / patch.weight;
69 return C;
70}
72/**
73 * Undamped accelerations (out.u, out.v, out.w) of every node for the
74 * displacement fields f, measured from the stress-free shape pl.w0. Bending
75 * forces are the gradient of the discrete Kirchhoff energy with the local
76 * bending stiffness D(h) at the nodes and at the cells (natural free-edge
77 * conditions); membrane forces come from bilinear cells with 2x2 Gauss
78 * integration, the cell's mean thickness, and the Marguerre form of the von
79 * Karman strains; the mount adds a spring to ground over its patch. Nodes
80 * outside the plate get zero.
81 */
82export function accelerations(pl: Plate, f: Fields, out: Fields): void {
83 const g = pl.grid, m = pl.mat, w0 = pl.w0;
84 const { nx, ny, hx, hy } = g;
85 const { u, v, w } = f;
86 const hx2 = hx * hx, hy2 = hy * hy;
87 const nu = m.nu;
88 const id = (i: number, j: number) => j * nx + i;
89 const Dn = (i: number, j: number) => bendingStiffness(m, pl.h[id(i, j)]);
90 const Dc = (ci: number, cj: number) => bendingStiffness(m, cellThickness(pl, ci, cj));
92 // A second difference exists where the node and both of its neighbors belong to the plate.
93 const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
94 const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
95 const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
96 const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
97 // Weighted bending moments at nodes; zero where the second difference is undefined.
98 const Ma = (i: number, j: number): number => {
99 if (!hasA(i, j)) return 0;
100 return nodeWeight(g, i, j) * Dn(i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
101 };
102 const Mb = (i: number, j: number): number => {
103 if (!hasB(i, j)) return 0;
104 return nodeWeight(g, i, j) * Dn(i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
105 };
106 // Weighted mixed second difference at a cell center; zero for cells outside the plate.
107 const Cc = (ci: number, cj: number): number => {
108 if (!insideCell(g, ci, cj)) return 0;
109 return (Dc(ci, cj) * (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)])) / (hx * hy);
110 };
112 const q = (hx * hy) / 4; // area per Gauss point
113 const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
114 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
116 // Adds the membrane force of cell (ci, cj) on its local node a to F = [Fu, Fv, Fw].
117 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
118 if (!insideCell(g, ci, cj)) return;
119 const hc = cellThickness(pl, ci, cj);
120 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
121 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
122 for (let b = 0; b < 4; b++) {
123 U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0[nn[b]];
124 }
125 for (let gp = 0; gp < 4; gp++) {
126 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
127 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
128 for (let b = 0; b < 4; b++) {
129 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
130 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
131 wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
132 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
133 }
134 // Marguerre strains: total slope (w0 + w) minus the stress-free slope of w0.
135 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
136 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
137 const gxy = uy + vx + wx * wy - w0x * w0y;
138 const Nxx = C * (exx + nu * eyy);
139 const Nyy = C * (eyy + nu * exx);
140 const Nxy = Gs * gxy;
141 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
142 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
143 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
144 }
145 };
147 const Km = mountStiffness(pl, mountPatch(pl));
148 const F = new Float64Array(3);
149 for (let j = 0; j < ny; j++) {
150 for (let i = 0; i < nx; i++) {
151 const k = id(i, j);
152 const mass = nodeMass(pl, i, j);
153 if (mass === 0) {
154 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
155 continue;
156 }
157 const bend =
158 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
159 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
160 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
161 F[0] = -Km[k] * u[k]; F[1] = -Km[k] * v[k]; F[2] = -hx * hy * bend - Km[k] * w[k];
162 cellForce(i, j, 0, F);
163 cellForce(i - 1, j, 1, F);
164 cellForce(i - 1, j - 1, 2, F);
165 cellForce(i, j - 1, 3, F);
166 out.u[k] = F[0] / mass;
167 out.v[k] = F[1] / mass;
168 out.w[k] = F[2] / mass;
169 }
170 }
171}
173/**
174 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
175 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
176 * transverse component, -2 sigmaIn v on the in-plane components, and the
177 * dashpot of the mount over its patch.
178 */
179export function addDamping(pl: Plate, vel: Fields, out: Fields): void {
180 const g = pl.grid;
181 const { sigma0, sigma1, sigmaIn } = pl.loss;
182 const { nx, ny, hx, hy } = g;
183 const id = (i: number, j: number) => j * nx + i;
184 const Cm = mountDamping(pl, mountPatch(pl));
185 for (let j = 0; j < ny; j++) {
186 for (let i = 0; i < nx; i++) {
187 const k = id(i, j);
188 const mass = nodeMass(pl, i, j);
189 if (mass === 0) continue;
190 const c = vel.w[k];
191 let lap = 0;
192 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
193 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
194 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
195 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
196 const cm = Cm[k] / mass;
197 out.u[k] -= (2 * (sigma0 + sigmaIn) + cm) * vel.u[k];
198 out.v[k] -= (2 * (sigma0 + sigmaIn) + cm) * vel.v[k];
199 out.w[k] += -(2 * sigma0 + cm) * c + 2 * sigma1 * lap;
200 }
201 }
202}
204/** Discrete potential energy (bending + membrane + mount) whose negative gradient is the force in `accelerations`. */
205export function potentialEnergy(pl: Plate, f: Fields): number {
206 const g = pl.grid, m = pl.mat, w0 = pl.w0;
207 const { nx, ny, hx, hy } = g;
208 const { u, v, w } = f;
209 const nu = m.nu;
210 const id = (i: number, j: number) => j * nx + i;
211 let Ub = 0;
212 for (let j = 0; j < ny; j++) {
213 for (let i = 0; i < nx; i++) {
214 let a = 0, b = 0;
215 if (insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
216 if (insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
217 Ub += nodeWeight(g, i, j) * bendingStiffness(m, pl.h[id(i, j)]) * (a * a + b * b + 2 * nu * a * b);
218 }
219 }
220 for (let cj = 0; cj < ny - 1; cj++) {
221 for (let ci = 0; ci < nx - 1; ci++) {
222 if (!insideCell(g, ci, cj)) continue;
223 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
224 Ub += 2 * (1 - nu) * bendingStiffness(m, cellThickness(pl, ci, cj)) * c * c;
225 }
226 }
227 Ub *= 0.5 * hx * hy;
229 let Um = 0;
230 const q = (hx * hy) / 4;
231 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
232 for (let cj = 0; cj < ny - 1; cj++) {
233 for (let ci = 0; ci < nx - 1; ci++) {
234 if (!insideCell(g, ci, cj)) continue;
235 const hc = cellThickness(pl, ci, cj);
236 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
237 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
238 for (let gp = 0; gp < 4; gp++) {
239 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
240 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
241 for (let b = 0; b < 4; b++) {
242 const z0 = w0[nn[b]];
243 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
244 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
245 wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
246 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
247 }
248 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
249 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
250 const gxy = uy + vx + wx * wy - w0x * w0y;
251 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
252 }
253 }
254 }
256 let Us = 0;
257 const Km = mountStiffness(pl, mountPatch(pl));
258 for (let k = 0; k < g.n; k++) if (Km[k] > 0) Us += 0.5 * Km[k] * (u[k] * u[k] + v[k] * v[k] + w[k] * w[k]);
259 return Ub + Um + Us;
260}
262/**
263 * Largest angular frequency of the linearized discrete system, by power
264 * iteration on a small full grid with the same spacing and a uniform
265 * thickness h (the largest of the plate, which bounds the bending frequency;
266 * the in-plane wave speed does not depend on the thickness). The extreme mode
267 * is a grid-scale mode, so neither the plate size nor its outline matters.
268 * `inPlaneMass` is the factor by which the in-plane inertia is scaled.
269 */
270export function maxAngularFrequency(g: Grid, m: Material, h: number, inPlaneMass = 1, nSmall = 12, iters = 300): number {
271 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
272 const pl = uniformPlate(gs, m, h);
273 const n = gs.n;
274 const x = fields(n), y = fields(n);
275 let seed = 12345;
276 const rnd = () => {
277 seed = (seed * 1664525 + 1013904223) >>> 0;
278 return seed / 4294967296 - 0.5;
279 };
280 for (let k = 0; k < n; k++) {
281 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
282 }
283 const eps = 1e-9; // keeps the nonlinear terms negligible
284 let lambda = 0;
285 for (let it = 0; it < iters; it++) {
286 let norm = 0;
287 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
288 norm = Math.sqrt(norm);
289 for (let k = 0; k < n; k++) {
290 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
291 }
292 accelerations(pl, x, y);
293 if (inPlaneMass !== 1) {
294 for (let k = 0; k < n; k++) { y.u[k] /= inPlaneMass; y.v[k] /= inPlaneMass; }
295 }
296 let dot = 0;
297 for (let k = 0; k < n; k++) dot += x.u[k] * y.u[k] + x.v[k] * y.v[k] + x.w[k] * y.w[k];
298 lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1} K
299 for (let k = 0; k < n; k++) {
300 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
301 }
302 }
303 return Math.sqrt(Math.max(lambda, 0));
304}