5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// Float64 reference implementation of the spatial discretization: the
2// forces of the geometrically nonlinear plate on the grid, the loss terms,
3// the discrete energy, and the largest frequency of the linearized system.
4// It mirrors the WGSL kernel of the full simulation and is used for the
5// time-step estimate, the tests, and the GPU-versus-CPU self-test; the
6// linear model assembles the same discretization as sparse matrices in
7// assemble.ts. It is not used in any time loop.
9import {
10 type Grid,
11 type Material,
12 type Plate,
13 type StandRates,
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 14 ROOM_TEMPERATURE,
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 15 bendingStiffness,
16 cellThickness,
17 insideCell,
18 insideNode,
19 membraneStiffness,
20 nodeMass,
21 nodeWeight,
22 shearStiffness,
23 standRates,
24} from "../params.ts";
25import { type SigmaFit } from "../loss.ts";
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 26import { feltStepConstants } from "../felt.ts";
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 27
28export interface Fields {
29 u: Float64Array;
30 v: Float64Array;
31 w: Float64Array;
32}
34export function fields(n: number): Fields {
35 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
36}
38/** A plate of uniform thickness on the given grid (for the tests and the time-step estimate). */
39export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
40 const hs = new Float64Array(g.n);
41 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;
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 42 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), stand: null, temperature: ROOM_TEMPERATURE };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 43}
45const G3 = 1 / Math.sqrt(3);
47/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
48export function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
49 const xi = gp & 1 ? G3 : -G3;
50 const eta = gp & 2 ? G3 : -G3;
51 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
52 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
53 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
54 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
55 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
56 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
57 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
58 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
59}
61/**
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 62 * Relaxed (static) force of the pads at a node with transverse displacement
63 * w: the felt's relaxed stiffness kInf times the stiffness at the bounce
64 * frequency, the rest of the felt's force being carried by the arms of its
65 * chain (addPadMemory). Without lift-off the pads are a linear spring of
66 * rate kInf kw. With lift-off each pad pushes only while compressed: the
67 * bottom pad by (preloadBottom - w), the top pad by (preloadTop + w); the
68 * constant part of the two, which balances the weight of the plate, is
69 * subtracted so that w = 0 is the rest state.
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 70 */
71export function padForce(sr: StandRates, k: number, w: number, liftoff: boolean): number {
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 72 const s = sr.chain.kInf;
73 if (!liftoff) return -s * sr.kw[k] * w;
74 const kb = s * sr.kwBottom[k];
75 const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 76 const bottom = kb * Math.max(sr.preloadBottom - w, 0);
77 const top = hasTop ? kb * Math.max(sr.preloadTop + w, 0) : 0;
78 const rest = kb * sr.preloadBottom - (hasTop ? kb * sr.preloadTop : 0);
79 return bottom - top - rest;
80}
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 82/** Number of pads in contact with the node at transverse displacement w (both pads without lift-off). */
83export function padsInContact(sr: StandRates, k: number, w: number, liftoff: boolean): number {
84 const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
85 if (!liftoff) return hasTop ? 2 : 1;
86 return (sr.preloadBottom - w > 0 ? 1 : 0) + (hasTop && sr.preloadTop + w > 0 ? 1 : 0);
87}
89/**
90 * The extensions of the arms of the felt's chain at every node, for the
91 * three fields: pad[(field * nArms + arm) * n + k] with field 0, 1, 2 for
92 * u, v, w. Zero is the relaxed felt, which is the initial state; the layout
93 * is that of the GPU buffer.
94 */
95export type PadState = Float64Array;
96export function padState(pl: Plate, sr: StandRates | null): PadState {
97 return new Float64Array(3 * (sr ? sr.chain.tau.length : 0) * pl.grid.n);
98}
100/**
101 * Adds to the accelerations `out` the force of the arms of the felt's
102 * chain, -k_node sum_i k_i p_i on each field, with the transverse force
103 * scaled by the number of pads in contact (the shear of the felt is applied
104 * whether or not a pad has lifted, as its static part is).
105 */
106export function addPadMemory(pl: Plate, sr: StandRates | null, f: Fields, pad: PadState, liftoff: boolean, out: Fields): void {
107 if (!sr) return;
108 const g = pl.grid, n = g.n, ch = sr.chain, nArms = ch.tau.length;
109 for (const k of sr.nodes) {
110 const mass = nodeMass(pl, k % g.nx, Math.floor(k / g.nx));
111 if (mass === 0) continue;
112 let su = 0, sv = 0, sw = 0;
113 for (let a = 0; a < nArms; a++) {
114 su += ch.k[a] * pad[a * n + k];
115 sv += ch.k[a] * pad[(nArms + a) * n + k];
116 sw += ch.k[a] * pad[(2 * nArms + a) * n + k];
117 }
118 const nc = padsInContact(sr, k, f.w[k], liftoff);
119 out.u[k] -= (sr.ku[k] * su) / mass;
120 out.v[k] -= (sr.ku[k] * sv) / mass;
121 out.w[k] -= (nc * sr.kwBottom[k] * sw) / mass;
122 }
123}
125/**
126 * Advances the arm extensions over one step of dt in which the fields moved
127 * by `dw` (the increments of u, v, w) from the displacements `f`: exact for
128 * a linear motion within the step. At a node that had left every pad at the
129 * start of the step the transverse arms relax at once.
130 */
131export function stepPadMemory(pl: Plate, sr: StandRates | null, dt: number, dw: Fields, f: Fields, pad: PadState, liftoff: boolean): void {
132 if (!sr) return;
133 const n = pl.grid.n, ch = sr.chain, nArms = ch.tau.length;
134 const { decay, growth } = feltStepConstants(ch, dt);
135 for (const k of sr.nodes) {
136 const lifted = padsInContact(sr, k, f.w[k], liftoff) === 0;
137 for (let a = 0; a < nArms; a++) {
138 pad[a * n + k] = pad[a * n + k] * decay[a] + dw.u[k] * growth[a];
139 pad[(nArms + a) * n + k] = pad[(nArms + a) * n + k] * decay[a] + dw.v[k] * growth[a];
140 const iw = (2 * nArms + a) * n + k;
141 pad[iw] = lifted ? 0 : pad[iw] * decay[a] + dw.w[k] * growth[a];
142 }
143 }
144}
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 146/**
147 * Undamped accelerations (out.u, out.v, out.w) of every node for the
148 * displacement fields f, measured from the stress-free shape pl.w0. Bending
149 * forces are the gradient of the discrete Kirchhoff energy with the local
150 * bending stiffness D(h) at the nodes and at the cells (natural free-edge
151 * conditions); membrane forces come from bilinear cells with 2x2 Gauss
152 * integration, the cell's mean thickness, and the Marguerre form of the von
153 * Karman strains; the pads of the stand add springs to ground over their
154 * annulus. Nodes outside the plate get zero.
155 */
156export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false, sr: StandRates | null = standRates(pl)): void {
157 const g = pl.grid, m = pl.mat, w0 = pl.w0;
158 const { nx, ny, hx, hy } = g;
159 const { u, v, w } = f;
160 const hx2 = hx * hx, hy2 = hy * hy;
161 const nu = m.nu;
162 const id = (i: number, j: number) => j * nx + i;
163 const Dn = (i: number, j: number) => bendingStiffness(m, pl.h[id(i, j)]);
164 const Dc = (ci: number, cj: number) => bendingStiffness(m, cellThickness(pl, ci, cj));
166 const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
167 const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
168 const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
169 const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
170 const Ma = (i: number, j: number): number => {
171 if (!hasA(i, j)) return 0;
172 return nodeWeight(g, i, j) * Dn(i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
173 };
174 const Mb = (i: number, j: number): number => {
175 if (!hasB(i, j)) return 0;
176 return nodeWeight(g, i, j) * Dn(i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
177 };
178 const Cc = (ci: number, cj: number): number => {
179 if (!insideCell(g, ci, cj)) return 0;
180 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);
181 };
183 const q = (hx * hy) / 4;
184 const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
185 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
187 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
188 if (!insideCell(g, ci, cj)) return;
189 const hc = cellThickness(pl, ci, cj);
190 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
191 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
192 for (let b = 0; b < 4; b++) {
193 U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0[nn[b]];
194 }
195 for (let gp = 0; gp < 4; gp++) {
196 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 197 let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 198 for (let b = 0; b < 4; b++) {
199 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
200 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 201 dwx += dNdx[b] * Wn[b]; dwy += dNdy[b] * Wn[b];
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 202 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
203 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 204 const wx = w0x + dwx, wy = w0y + dwy;
205 // See step.wgsl: every term carries the deviation, so that the single
206 // precision of the GPU has nothing large to cancel.
207 const exx = ux + dwx * (w0x + 0.5 * dwx);
208 const eyy = vy + dwy * (w0y + 0.5 * dwy);
209 const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 210 const Nxx = C * (exx + nu * eyy);
211 const Nyy = C * (eyy + nu * exx);
212 const Nxy = Gs * gxy;
213 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
214 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
215 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
216 }
217 };
219 const F = new Float64Array(3);
220 for (let j = 0; j < ny; j++) {
221 for (let i = 0; i < nx; i++) {
222 const k = id(i, j);
223 const mass = nodeMass(pl, i, j);
224 if (mass === 0) {
225 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
226 continue;
227 }
228 const bend =
229 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
230 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
231 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
232 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
233 if (sr) {
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 234 F[0] -= sr.chain.kInf * sr.ku[k] * u[k]; F[1] -= sr.chain.kInf * sr.ku[k] * v[k];
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 235 F[2] += padForce(sr, k, w[k], liftoff);
236 }
237 cellForce(i, j, 0, F);
238 cellForce(i - 1, j, 1, F);
239 cellForce(i - 1, j - 1, 2, F);
240 cellForce(i, j - 1, 3, F);
241 out.u[k] = F[0] / mass;
242 out.v[k] = F[1] / mass;
243 out.w[k] = F[2] / mass;
244 }
245 }
246}
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 248/** The loss terms of the full simulation on the plate itself: the fitted sigma model on the transverse motion and an in-plane rate. The felt's loss is in its chain (addPadMemory). */
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 249export interface FullLoss {
250 sigma0: number; // [1/s], applied as 2 sigma0 on all components
251 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
252 sigmaIn: number; // [1/s], extra rate on the in-plane motion
253}
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 255/** The loss terms of the full simulation from the fitted sigma model. */
256export function fullLoss(fit: SigmaFit, inPlaneRate: number): FullLoss {
257 return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 258}
260/**
261 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
262 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 263 * transverse component, and -2 sigmaIn v on the in-plane components.
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 264 */
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 265export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields): void {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 266 const g = pl.grid;
267 const { sigma0, sigma1, sigmaIn } = loss;
268 const { nx, ny, hx, hy } = g;
269 const id = (i: number, j: number) => j * nx + i;
270 for (let j = 0; j < ny; j++) {
271 for (let i = 0; i < nx; i++) {
272 const k = id(i, j);
273 const mass = nodeMass(pl, i, j);
274 if (mass === 0) continue;
275 const c = vel.w[k];
276 let lap = 0;
277 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
278 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
279 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
280 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 281 out.u[k] -= 2 * (sigma0 + sigmaIn) * vel.u[k];
282 out.v[k] -= 2 * (sigma0 + sigmaIn) * vel.v[k];
283 out.w[k] += -2 * sigma0 * c + 2 * sigma1 * lap;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 284 }
285 }
286}
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 288/** Discrete potential energy (bending + membrane + the relaxed pads, linear) whose negative gradient is the force in `accelerations`. */
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 289export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = standRates(pl)): number {
290 const g = pl.grid, m = pl.mat, w0 = pl.w0;
291 const { nx, ny, hx, hy } = g;
292 const { u, v, w } = f;
293 const nu = m.nu;
294 const id = (i: number, j: number) => j * nx + i;
295 let Ub = 0;
296 for (let j = 0; j < ny; j++) {
297 for (let i = 0; i < nx; i++) {
298 let a = 0, b = 0;
299 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);
300 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);
301 Ub += nodeWeight(g, i, j) * bendingStiffness(m, pl.h[id(i, j)]) * (a * a + b * b + 2 * nu * a * b);
302 }
303 }
304 for (let cj = 0; cj < ny - 1; cj++) {
305 for (let ci = 0; ci < nx - 1; ci++) {
306 if (!insideCell(g, ci, cj)) continue;
307 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
308 Ub += 2 * (1 - nu) * bendingStiffness(m, cellThickness(pl, ci, cj)) * c * c;
309 }
310 }
311 Ub *= 0.5 * hx * hy;
313 let Um = 0;
314 const q = (hx * hy) / 4;
315 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
316 for (let cj = 0; cj < ny - 1; cj++) {
317 for (let ci = 0; ci < nx - 1; ci++) {
318 if (!insideCell(g, ci, cj)) continue;
319 const hc = cellThickness(pl, ci, cj);
320 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
321 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
322 for (let gp = 0; gp < 4; gp++) {
323 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 324 let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 325 for (let b = 0; b < 4; b++) {
326 const z0 = w0[nn[b]];
327 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
328 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 329 dwx += dNdx[b] * w[nn[b]]; dwy += dNdy[b] * w[nn[b]];
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 330 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
331 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 332 const exx = ux + dwx * (w0x + 0.5 * dwx);
333 const eyy = vy + dwy * (w0y + 0.5 * dwy);
334 const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 335 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
336 }
337 }
338 }
340 let Us = 0;
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 341 if (sr) for (const k of sr.nodes) Us += 0.5 * sr.chain.kInf * (sr.ku[k] * (u[k] * u[k] + v[k] * v[k]) + sr.kw[k] * w[k] * w[k]);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 342 return Ub + Um + Us;
343}
345/**
346 * Largest angular frequency of the linearized discrete system, by power
347 * iteration on a small full grid with the same spacing and a uniform
348 * thickness h (the largest of the plate, which bounds the bending frequency;
349 * the in-plane wave speed does not depend on the thickness). The extreme
350 * mode is a grid-scale mode, so neither the plate size nor its outline
351 * matters.
352 */
353export function maxAngularFrequency(g: Grid, m: Material, h: number, nSmall = 12, iters = 300): number {
354 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
355 const pl = uniformPlate(gs, m, h);
356 const n = gs.n;
357 const x = fields(n), y = fields(n);
358 let seed = 12345;
359 const rnd = () => {
360 seed = (seed * 1664525 + 1013904223) >>> 0;
361 return seed / 4294967296 - 0.5;
362 };
363 for (let k = 0; k < n; k++) {
364 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
365 }
366 const eps = 1e-9;
367 let lambda = 0;
368 for (let it = 0; it < iters; it++) {
369 let norm = 0;
370 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
371 norm = Math.sqrt(norm);
372 for (let k = 0; k < n; k++) {
373 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
374 }
375 accelerations(pl, x, y, false, null);
376 let dot = 0;
377 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];
378 lambda = -dot / (eps * eps);
379 for (let k = 0; k < n; k++) {
380 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
381 }
382 }
383 return Math.sqrt(Math.max(lambda, 0));
384}