1// The float64 reference discretization: symmetry and rigid modes of the
2// linear stiffness, forces as the gradient of the energy, dissipativity,
3// thickness scaling, and the free-plate frequencies.
4import { test } from "node:test";
5import assert from "node:assert/strict";
6import { makeGrid, withMask, type Grid, type Plate, nodeMasses, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
7import { uniformPlate, accelerations, addDamping, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
8import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
9import { fullLoss } from "../src/model/reference.ts";
11const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
13function rnd(seed: number): () => number {
14 let s = seed >>> 0;
15 return () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
16}
17function randomFields(pl: Plate, amp: number, seed: number): Fields {
18 const g = pl.grid, r = rnd(seed), f = fields(g.n);
19 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
20 if (!insideNode(g, i, j)) continue;
21 const k = j * g.nx + i;
22 f.u[k] = 0.02 * amp * r(); f.v[k] = 0.02 * amp * r(); f.w[k] = amp * r();
23 }
24 return f;
25}
26function bowed(g: Grid, sag: number): Float64Array {
27 const w0 = new Float64Array(g.n);
28 const Lx = (g.nx - 1) * g.hx;
29 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = sag * (1 - (2 * nodeX(g, i) / Lx) ** 2) + 0.3 * sag * Math.sin(3 * nodeY(g, j) / Lx);
30 return w0;
31}
32test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
33 const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
34 const g1 = makeGrid(0.12, 0.09, 0.015);
35 cases.push({ name: "uniform rectangle", pl: uniformPlate(g1, steel, 1e-3), fields: 3 });
36 const tapered = uniformPlate(g1, steel, 1e-3);
37 for (let k = 0; k < g1.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.5 * (k % g1.nx) / g1.nx);
38 cases.push({ name: "tapered rectangle", pl: tapered, fields: 3 });
39 const gd = withMask(makeGrid(0.16, 0.16, 0.01), (x, y) => x * x + y * y <= 0.08 ** 2 && x * x + y * y >= 0.015 ** 2);
40 cases.push({ name: "disc with hole", pl: uniformPlate(gd, steel, 1e-3), fields: 1 });
41 for (const c of cases) {
42 const sys = assemble(c.pl, c.fields);
43 // Symmetry is built into the banded storage; check the rigid modes: K x = 0 for translations and rotations.
44 const g = c.pl.grid;
45 const nullVectors: Float64Array[] = [];
46 const mk = (fn: (i: number, j: number, comp: number) => number) => {
47 const x = new Float64Array(sys.n);
48 for (let d = 0; d < sys.n; d++) {
49 const k = sys.dofNode[d];
50 const comp = sys.fields === 1 ? 2 : d % 3;
51 x[d] = fn(k % g.nx, Math.floor(k / g.nx), comp);
52 }
53 return x;
54 };
55 nullVectors.push(mk((_i, _j, c) => (c === 2 ? 1 : 0)));
56 nullVectors.push(mk((i, _j, c) => (c === 2 ? nodeX(g, i) : 0)));
57 nullVectors.push(mk((_i, j, c) => (c === 2 ? nodeY(g, j) : 0)));
58 if (c.fields === 3) {
59 nullVectors.push(mk((_i, _j, c) => (c === 0 ? 1 : 0)));
60 nullVectors.push(mk((_i, _j, c) => (c === 1 ? 1 : 0)));
61 nullVectors.push(mk((i, j, c) => (c === 0 ? -nodeY(g, j) : c === 1 ? nodeX(g, i) : 0)));
62 }
63 const y = new Float64Array(sys.n);
64 for (const x of nullVectors) {
65 const { stiffnessMatvec } = require_assemble();
66 stiffnessMatvec(sys, x, y);
67 let ny = 0, nx = 0;
68 for (let d = 0; d < sys.n; d++) { ny += y[d] * y[d]; nx += x[d] * x[d]; }
69 const scale = Math.abs(stiffnessEntry(sys, 0, 0)) * Math.sqrt(nx);
70 assert.ok(Math.sqrt(ny) < 1e-9 * scale, `${c.name}: rigid mode residual ${Math.sqrt(ny) / scale}`);
71 }
72 // A non-rigid vector must not be in the null space.
73 const z = mk((i, j, c) => (c === 2 ? Math.sin(i * 0.9) * Math.cos(j * 1.1) : 0));
74 const { stiffnessMatvec } = require_assemble();
75 stiffnessMatvec(sys, z, y);
76 assert.ok(y.some((v) => Math.abs(v) > 1e-3), `${c.name}: bending is not trivial`);
77 }
78});
80import * as assembleModule from "../src/model/assemble.ts";
81function require_assemble(): typeof assembleModule { return assembleModule; }
83test("forces are the exact gradient of the discrete energy (flat and bowed, tapered)", () => {
84 const g = makeGrid(0.1, 0.08, 0.01);
85 const variants: { name: string; pl: Plate }[] = [];
86 variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
87 variants.push({ name: "bowed free", pl: uniformPlate(g, steel, 1e-3, bowed(g, 3e-3)) });
88 const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
89 for (let k = 0; k < g.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.4 * (k % g.nx) / g.nx);
90 variants.push({ name: "bowed tapered", pl: tapered });
91 for (const v of variants) {
92 const pl = v.pl;
93 const f = randomFields(pl, 2.5e-3, 7);
94 const masses = nodeMasses(pl);
95 const acc = fields(g.n);
96 accelerations(pl, f, acc);
97 const dir = randomFields(pl, 1, 11);
98 // Directional derivative of U along dir versus -F . dir.
99 const eps = 1e-7;
100 const shift = (s: number): Fields => ({ u: f.u.map((x, k) => x + s * dir.u[k]), v: f.v.map((x, k) => x + s * dir.v[k]), w: f.w.map((x, k) => x + s * dir.w[k]) });
101 const dU = (potentialEnergy(pl, shift(eps)) - potentialEnergy(pl, shift(-eps))) / (2 * eps);
102 let fd = 0;
103 for (let k = 0; k < g.n; k++) fd -= masses[k] * (acc.u[k] * dir.u[k] + acc.v[k] * dir.v[k] + acc.w[k] * dir.w[k]);
104 assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
105 }
106});
108test("the loss terms only remove energy", () => {
109 const g = makeGrid(0.1, 0.08, 0.01);
110 const pl = uniformPlate(g, steel, 1e-3);
111 const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, 20);
112 const masses = nodeMasses(pl);
113 for (const seed of [1, 2, 3]) {
114 const vel = randomFields(pl, 1, seed);
115 const acc = fields(g.n);
116 addDamping(pl, loss, vel, acc);
117 let power = 0;
118 for (let k = 0; k < g.n; k++) power += masses[k] * (acc.u[k] * vel.u[k] + acc.v[k] * vel.v[k] + acc.w[k] * vel.w[k]);
119 assert.ok(power < 0, `power ${power}`);
120 }
121});
123test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
124 const g = makeGrid(0.1, 0.1, 0.01);
125 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
126 // The largest frequency is set by the in-plane waves, independent of h.
127 assert.ok(Math.abs(w1 / w2 - 1) < 1e-6, `in-plane omega_max ${w1} vs ${w2}`);
128 const { stiffnessEntry: ke } = assembleModule;
129 const s1 = assemble(uniformPlate(g, steel, 1e-3)), s2 = assemble(uniformPlate(g, steel, 2e-3));
130 const c = Math.floor(s1.n / 2);
131 // Bending stiffness scales as h^3 and the mass as h: omega^2 scales as h^2.
132 assert.ok(Math.abs(ke(s2, c, c) / ke(s1, c, c) - 8) < 1e-9);
133 assert.ok(Math.abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12);
134 assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
135});