1import { test } from "node:test";
2import assert from "node:assert/strict";
3import { accelerations, potentialEnergy, maxAngularFrequency, type Fields } from "../src/sim/cpuRef.ts";
4import { type Grid, type Material, countInside, nodeWeight, withEllipseMask } from "../src/params.ts";
6const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
8function fields(n: number): Fields {
9 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
10}
12/** Jacobi eigenvalues of a symmetric matrix (row-major, n x n). */
13function symEig(A: Float64Array, n: number): number[] {
14 const a = Float64Array.from(A);
15 for (let sweep = 0; sweep < 60; sweep++) {
16 let off = 0;
17 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;
18 if (off < 1e-22) break;
19 for (let p = 0; p < n; p++) {
20 for (let q = p + 1; q < n; q++) {
21 const apq = a[p * n + q];
22 if (Math.abs(apq) < 1e-300) continue;
23 const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);
24 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
25 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
26 for (let k = 0; k < n; k++) {
27 const akp = a[k * n + p], akq = a[k * n + q];
28 a[k * n + p] = c * akp - s * akq;
29 a[k * n + q] = s * akp + c * akq;
30 }
31 for (let k = 0; k < n; k++) {
32 const apk = a[p * n + k], aqk = a[q * n + k];
33 a[p * n + k] = c * apk - s * aqk;
34 a[q * n + k] = s * apk + c * aqk;
35 }
36 }
37 }
38 }
39 const ev: number[] = [];
40 for (let k = 0; k < n; k++) ev.push(a[k * n + k]);
41 return ev.sort((x, y) => x - y);
42}
44/** Mass-symmetrized linear stiffness of one block ("w" or "uv") over the nodes of the sheet. */
45function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
46 const nodes: number[] = [];
47 for (let k = 0; k < g.n; k++) if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) > 0) nodes.push(k);
48 const N = nodes.length;
49 const n = block === "w" ? N : 2 * N;
50 const K = new Float64Array(n * n);
51 const x = fields(g.n), y = fields(g.n);
52 const eps = 1e-9;
53 const mass = (d: number) => mat.rhoh * g.hx * g.hy * nodeWeight(g, nodes[d] % g.nx, Math.floor(nodes[d] / g.nx));
54 for (let col = 0; col < n; col++) {
55 x.u.fill(0); x.v.fill(0); x.w.fill(0);
56 if (block === "w") x.w[nodes[col]] = eps;
57 else if (col < N) x.u[nodes[col]] = eps;
58 else x.v[nodes[col - N]] = eps;
59 accelerations(g, mat, x, y);
60 const mc = mass(col % N);
61 for (let row = 0; row < n; row++) {
62 const acc = block === "w" ? y.w[nodes[row]] : row < N ? y.u[nodes[row]] : y.v[nodes[row - N]];
63 const mr = mass(row % N);
64 // K_sym = M^{1/2} (M^{-1} K) M^{-1/2}
65 K[row * n + col] = (-acc / eps) * Math.sqrt(mr / mc);
66 }
67 }
68 return { K, n };
69}
71function checkSymmetricWithThreeRigidModes(K: Float64Array, n: number): number[] {
72 let asym = 0, scale = 0;
73 for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
74 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
75 scale = Math.max(scale, Math.abs(K[r * n + c]));
76 }
77 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
78 const ev = symEig(K, n);
79 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
80 assert.ok(ev[3] > 1e-4 * ev[n - 1], `fourth mode is not rigid (${ev[3]} vs ${ev[n - 1]})`);
81 return ev;
82}
84/** Relative gradient mismatch between the forces and the finite-difference gradient of the energy. */
85function gradientMismatch(g: Grid, f: Fields, w0?: Float64Array): number {
86 const acc = fields(g.n);
87 accelerations(g, mat, f, acc, w0);
88 const delta = 1e-7;
89 let worst = 0;
90 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
91 for (let k = 0; k < g.n; k += 5) {
92 const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
93 if (mass === 0) continue;
94 const orig = arr[k];
95 arr[k] = orig + delta;
96 const Ep = potentialEnergy(g, mat, f, w0);
97 arr[k] = orig - delta;
98 const Em = potentialEnergy(g, mat, f, w0);
99 arr[k] = orig;
100 const Ffd = -(Ep - Em) / (2 * delta);
101 const Fan = accArr[k] * mass;
102 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
103 }
104 }
105 return worst;
106}
108test("bending stiffness is symmetric and has exactly three rigid modes", () => {
109 const nx = 9, ny = 7;
110 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
111 const { K, n } = stiffness(g, "w");
112 let asym = 0, scale = 0;
113 for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
114 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
115 scale = Math.max(scale, Math.abs(K[r * n + c]));
116 }
117 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
118 const ev = symEig(K, n);
119 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
120 assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
121});
123test("membrane stiffness is symmetric and has exactly three rigid modes", () => {
124 const nx = 7, ny = 6;
125 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
126 const { K, n } = stiffness(g, "uv");
127 let asym = 0, scale = 0;
128 for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
129 asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
130 scale = Math.max(scale, Math.abs(K[r * n + c]));
131 }
132 assert.ok(asym < 1e-6 * scale, `asymmetry ${asym}`);
133 const ev = symEig(K, n);
134 assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
135 assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
136});
138test("free square plate frequencies match the Leissa values", () => {
139 const nx = 21, ny = 21, a = 0.3;
140 const g: Grid = { nx, ny, hx: a / (nx - 1), hy: a / (ny - 1), n: nx * ny };
141 const { K, n } = stiffness(g, "w");
142 const ev = symEig(K, n);
143 const lam2 = ev.slice(3, 8).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
144 // Completely free square plate, nu = 0.3 (Leissa / Gorman): first five nonzero modes.
145 const ref = [13.468, 19.596, 24.27, 34.801, 34.801];
146 for (let k = 0; k < ref.length; k++) {
147 const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
148 assert.ok(relErr < 0.03, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
149 }
150});
152test("forces are the negative gradient of the discrete energy (nonlinear state)", () => {
153 const nx = 8, ny = 7;
154 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
155 const f = fields(g.n), acc = fields(g.n);
156 let seed = 7;
157 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
158 for (let k = 0; k < g.n; k++) {
159 f.w[k] = 0.02 * rnd(); // ~ 25 thicknesses, strongly nonlinear
160 f.u[k] = 1e-4 * rnd();
161 f.v[k] = 1e-4 * rnd();
162 }
163 accelerations(g, mat, f, acc);
164 const delta = 1e-7;
165 let worst = 0;
166 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
167 for (let k = 0; k < g.n; k += 5) {
168 const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % nx, Math.floor(k / nx));
169 const orig = arr[k];
170 arr[k] = orig + delta;
171 const Ep = potentialEnergy(g, mat, f);
172 arr[k] = orig - delta;
173 const Em = potentialEnergy(g, mat, f);
174 arr[k] = orig;
175 const Ffd = -(Ep - Em) / (2 * delta);
176 const Fan = accArr[k] * mass;
177 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
178 }
179 }
180 assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);
181});
183test("maximum frequency estimate is close to the in-plane CFL estimate", () => {
184 const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };
185 const wmax = maxAngularFrequency(g, mat);
186 const cp = Math.sqrt(mat.C / mat.rhoh);
187 const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)
188 assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);
189});
191test("forces are the energy gradient with an initial-shape imperfection", () => {
192 const nx = 8, ny = 7;
193 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
194 const f = fields(g.n), acc = fields(g.n);
195 const w0 = new Float64Array(g.n);
196 let seed = 11;
197 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
198 for (let k = 0; k < g.n; k++) {
199 const x = (k % nx) / (nx - 1) - 0.5;
200 w0[k] = 0.003 * (1 - 4 * x * x); // cylindrical bow, 3 mm
201 f.w[k] = 0.01 * rnd();
202 f.u[k] = 1e-4 * rnd();
203 f.v[k] = 1e-4 * rnd();
204 }
205 accelerations(g, mat, f, acc, w0);
206 const delta = 1e-7;
207 let worst = 0;
208 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
209 for (let k = 0; k < g.n; k += 5) {
210 const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % nx, Math.floor(k / nx));
211 const orig = arr[k];
212 arr[k] = orig + delta;
213 const Ep = potentialEnergy(g, mat, f, w0);
214 arr[k] = orig - delta;
215 const Em = potentialEnergy(g, mat, f, w0);
216 arr[k] = orig;
217 const Ffd = -(Ep - Em) / (2 * delta);
218 const Fan = accArr[k] * mass;
219 worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
220 }
221 }
222 assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);
223 // The stress-free shape itself carries no membrane energy.
224 const rest = fields(g.n);
225 assert.ok(potentialEnergy(g, mat, rest, w0) < 1e-12, "w0 is stress free");
226});
228// ---------- Masked (elliptical) grids ----------
230test("elliptical sheet: bending and membrane stiffness are symmetric with exactly three rigid modes", () => {
231 const nx = 13, ny = 10;
232 const g = withEllipseMask({ nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny });
233 assert.ok(countInside(g) < g.n && countInside(g) > g.n / 2, "mask removes some nodes");
234 for (let k = 0; k < g.n; k++) {
235 // Every node of the sheet carries mass (no whiskers left after the mask is trimmed).
236 if (g.mask![k]) assert.ok(nodeWeight(g, k % nx, Math.floor(k / nx)) > 0, `node ${k} has no mass`);
237 }
238 const b = stiffness(g, "w");
239 checkSymmetricWithThreeRigidModes(b.K, b.n);
240 const m = stiffness(g, "uv");
241 checkSymmetricWithThreeRigidModes(m.K, m.n);
242});
244test("elliptical sheet: forces are the energy gradient, with and without an initial shape", () => {
245 const nx = 12, ny = 9;
246 const g = withEllipseMask({ nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny });
247 const f = fields(g.n);
248 const w0 = new Float64Array(g.n);
249 let seed = 5;
250 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
251 for (let k = 0; k < g.n; k++) {
252 const x = (k % nx) / (nx - 1) - 0.5;
253 w0[k] = 0.003 * (1 - 4 * x * x);
254 f.w[k] = 0.01 * rnd();
255 f.u[k] = 1e-4 * rnd();
256 f.v[k] = 1e-4 * rnd();
257 }
258 assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
259 assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
260});
262test("free circular plate frequencies approach the Leissa values", () => {
263 const nx = 25, a = 0.15; // radius; 24 cells across the diameter
264 const g = withEllipseMask({ nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx });
265 const { K, n } = stiffness(g, "w");
266 const ev = symEig(K, n);
267 const lam2 = ev.slice(3, 12).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
268 // Completely free circular plate, nu = 0.3 (Leissa, Table 2.4): modes (2,0) x2, (0,1),
269 // (3,0) x2, (1,1) x2, (4,0) x2 by nodal diameters and circles. The staircase edge is
270 // first-order accurate: about 2 percent here, under 1 percent at 60 cells across.
271 const ref = [5.358, 5.358, 9.003, 12.439, 12.439, 20.475, 20.475, 21.835, 21.835];
272 for (let k = 0; k < ref.length; k++) {
273 const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
274 assert.ok(relErr < 0.04, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
275 }
276});