1import { describe, expect, it } from "vitest";
2import { LIBRARY_FILES, buildCompareBundle, isPlainStruct, parseCompareResult } from "../src/index";
3import { maxAbs, runner, solve, testSpec } from "./helpers";
5/** Run a MATLAB snippet against the library and return the struct it leaves in `out`. */
6async function script(source: string) {
7 const run = await runner.run({
8 bundle: {
9 files: [{ name: "c.m", source }, ...LIBRARY_FILES],
10 entry: "c.m",
11 resultVariable: "out",
12 },
13 });
14 if (!run.ok) throw new Error(`${run.error}\n${run.output.join("")}`);
15 if (!isPlainStruct(run.value)) throw new Error("unexpected shape");
16 return run.value.fields as Record<string, number>;
17}
19describe("free edges", () => {
20 it("assembles a symmetric stiffness whose null space is exactly the rigid-body motions and whose interior is the 13-point stencil", async () => {
21 const out = await script(`
22 nx = 19; ny = 13; dx = 0.3/(nx-1); dy = 0.2/(ny-1); nu = 0.3;
23 dom = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), nx, ny);
24 [L, K, M, idx] = plate_matrices(dom, 'free', nu);
25 x = linspace(0, 0.3, nx); y = linspace(0, 0.2, ny); [X, Y] = meshgrid(x, y);
26 out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
27 out.null1 = max(abs(K * ones(nx*ny, 1)));
28 out.nullx = max(abs(K * X(:))) * dx^4;
29 out.nully = max(abs(K * Y(:))) * dx^4;
30 out.lnull = max(abs(L * ones(nx*ny, 1)));
31 % a curved field: the interior rows must be the biharmonic stencil, per unit mass
32 W = sin(3*X) .* cos(5*Y) + X.^2 .* Y;
33 Kw = reshape(K * W(:), ny, nx); Mv = reshape(full(diag(M)), ny, nx);
34 Lw = plate_laplacian(plate_pad(W, 'simply-supported'), dx, dy); % pad only feeds the outer ring, dropped below
35 Bw = (Lw(3:end,2:end-1) - 2*Lw(2:end-1,2:end-1) + Lw(1:end-2,2:end-1))/dy^2 + (Lw(2:end-1,3:end) - 2*Lw(2:end-1,2:end-1) + Lw(2:end-1,1:end-2))/dx^2;
36 inner = Kw(3:end-2, 3:end-2) ./ Mv(3:end-2, 3:end-2);
37 ref = Bw(2:end-1, 2:end-1);
38 out.interior = norm(inner(:) - ref(:)) / norm(ref(:));
39 out.m_inner = Mv(5, 7); out.m_edge = Mv(1, 7); out.m_corner = Mv(1, 1);
40 out.nunk = numel(idx);
41 `);
42 expect(out.symm).toBeLessThan(1e-12);
43 expect(out.null1).toBeLessThan(1e-6);
44 expect(out.nullx).toBeLessThan(1e-9);
45 expect(out.nully).toBeLessThan(1e-9);
46 expect(out.lnull).toBeLessThan(1e-9);
47 expect(out.interior).toBeLessThan(1e-10);
48 expect(out.m_inner).toBe(1);
49 expect(out.m_edge).toBe(0.5);
50 expect(out.m_corner).toBe(0.25);
51 expect(out.nunk).toBe(19 * 13);
52 });
54 it("reproduces the frequencies of a completely free square plate (Leissa, nu = 0.3)", async () => {
55 // omega a^2 sqrt(rho h / D) for the first three elastic modes of the free square plate,
56 // Ritz values tabulated by Leissa (1969): 13.49, 19.79, 24.43. The energy method converges
57 // to them from below at about first order in dx (the edge treatment), so a 17-node grid
58 // should land within a couple of percent, and the three rigid modes must sit at zero.
59 const out = await script(`
60 n = 17; a = 1; h = a/(n-1); nu = 0.3;
61 dom = plate_domain(struct('shape', 'rectangle', 'Lx', a, 'Ly', a), n, n);
62 [~, K, M] = plate_matrices(dom, 'free', nu);
63 Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh;
64 S = (S + S') / 2;
65 ev = sort(eig(S));
66 out.zero3 = max(abs(ev(1:3))) / ev(end);
67 out.f1 = sqrt(ev(4)) * a^2; out.f2 = sqrt(ev(5)) * a^2; out.f3 = sqrt(ev(6)) * a^2;
68 out.lam_ratio = ev(end) / (4/h^2 + 4/h^2)^2;
69 `);
70 expect(out.zero3).toBeLessThan(1e-9);
71 expect(Math.abs(out.f1 - 13.49) / 13.49).toBeLessThan(0.02);
72 // and it lies below the simply supported square's 2 pi^2 = 19.74: free edges are the softest support
73 expect(out.f1).toBeLessThan(2 * Math.PI ** 2);
74 expect(Math.abs(out.f2 - 19.79) / 19.79).toBeLessThan(0.03);
75 expect(Math.abs(out.f3 - 24.43) / 24.43).toBeLessThan(0.03);
76 // the explicit scheme's CFL bound, derived for the interior stencil, remains an upper bound
77 expect(out.lam_ratio).toBeLessThan(1);
78 expect(out.lam_ratio).toBeGreaterThan(0.9);
79 });
81 it("runs the explicit scheme on a free plate: stable below the CFL limit, no rigid-body drift, corners stress-free", async () => {
82 const spec = testSpec();
83 spec.bc = "free";
84 spec.fd.cfl = 0.98;
85 const { result } = await solve(spec, "fd");
86 expect(result.info.stable).toBe(true);
87 expect(result.probe.w.every(Number.isFinite)).toBe(true);
88 expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
89 expect(result.info.notes.some((n) => /rigid-body/.test(n))).toBe(true);
90 // the last frame carries no mean displacement and no mean tilt in the mass inner product
91 // (trapezoidal node weights: 1 inside, 1/2 on an edge, 1/4 at a corner), which is the
92 // inner product the projector uses
93 const { nx, ny, nf, data } = result.frames.w;
94 const k = nf - 1;
95 let sum = 0,
96 sx = 0,
97 sy = 0,
98 peak = 0,
99 mass = 0;
100 for (let j = 0; j < nx; j++)
101 for (let i = 0; i < ny; i++) {
102 const v = data[i + j * ny + k * nx * ny];
103 const m = (j === 0 || j === nx - 1 ? 0.5 : 1) * (i === 0 || i === ny - 1 ? 0.5 : 1);
104 mass += m;
105 sum += m * v;
106 sx += m * v * (j - (nx - 1) / 2);
107 sy += m * v * (i - (ny - 1) / 2);
108 peak = Math.max(peak, Math.abs(v));
109 }
110 expect(peak).toBeGreaterThan(0);
111 expect(Math.abs(sum) / mass).toBeLessThan(1e-9 * peak);
112 expect(Math.abs(sx) / (mass * nx)).toBeLessThan(1e-9 * peak);
113 expect(Math.abs(sy) / (mass * ny)).toBeLessThan(1e-9 * peak);
114 // at a free corner both curvatures and the twist vanish, so every stress there is zero
115 const vm = result.frames.vm;
116 for (const [i, j] of [
117 [0, 0],
118 [0, nx - 1],
119 [ny - 1, 0],
120 [ny - 1, nx - 1],
121 ])
122 expect(vm.data[i + j * ny + k * nx * ny]).toBe(0);
123 // and the edges move, unlike a supported plate
124 let edgePeak = 0;
125 for (let j = 0; j < nx; j++)
126 edgePeak = Math.max(edgePeak, Math.abs(data[0 + j * ny + k * nx * ny]));
127 expect(edgePeak).toBeGreaterThan(0.05 * peak);
128 });
130 it("diverges just above the CFL limit, as the supported plate does", async () => {
131 const spec = testSpec();
132 spec.bc = "free";
133 spec.fd.cfl = 1.04;
134 const { result } = await solve(spec, "fd");
135 expect(result.info.stable).toBe(false);
136 expect(result.probe.w.every(Number.isFinite)).toBe(true);
137 });
139 it("agrees between the implicit and explicit schemes", async () => {
140 const spec = testSpec();
141 spec.bc = "free";
142 spec.implicit = {
143 nx: spec.fd.nx,
144 ny: spec.fd.ny,
145 duration: spec.fd.duration,
146 steps_per_sample: 2,
147 theta: 0.25,
148 };
149 const imp = await solve(spec, "implicit");
150 const exp = await solve(spec, "fd");
151 expect(imp.result.info.stable).toBe(true);
152 expect(imp.result.info.extra.nunknowns).toBe(spec.fd.nx * spec.fd.ny);
153 const cmp = await runner.run({
154 bundle: buildCompareBundle(),
155 inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
156 });
157 expect(cmp.ok).toBe(true);
158 if (!cmp.ok) return;
159 const c = parseCompareResult(cmp.value);
160 expect(c.w_rel_l2).toBeLessThan(0.03);
161 });
163 it("is refused by the modal method with a legible message", async () => {
164 const spec = testSpec();
165 spec.bc = "free";
166 await expect(solve(spec, "modal")).rejects.toThrow(/simply supported/);
167 });
168});
170describe("pinned centre", () => {
171 it("removes the centre node from the system and leaves exactly the two rotations about it at zero frequency", async () => {
172 const out = await script(`
173 n = 15; a = 1; h = a/(n-1); nu = 0.3;
174 dom = plate_domain(struct('shape', 'rectangle', 'Lx', a, 'Ly', a), n, n);
175 [L, K, M, idx, fixed] = plate_matrices(dom, 'pinned-center', nu);
176 ic = (floor(n/2) + 1) + floor(n/2) * n;
177 out.fixed = fixed; out.ic = ic;
178 out.row = max(abs(K(ic, :))); out.col = max(abs(K(:, ic))); out.lrow = max(abs(L(ic, :)));
179 out.mdiag = full(M(ic, ic));
180 x = linspace(0, a, n); [X, Y] = meshgrid(x, x);
181 out.rotx = max(abs(K * (X(:) - X(ic)))) * h^4; out.roty = max(abs(K * (Y(:) - Y(ic)))) * h^4;
182 out.const = max(abs(K * ones(n*n, 1))) * h^4;
183 % the pinned node's own row and column are zero, so drop it before looking at the spectrum
184 keep = true(n*n, 1); keep(ic) = false;
185 Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh; S = (S + S') / 2; S = S(keep, keep);
186 ev = sort(eig(S));
187 out.zero2 = max(abs(ev(1:2))) / ev(end);
188 out.third = ev(3) / ev(end);
189 out.f1 = sqrt(ev(3)) * a^2; out.f2 = sqrt(ev(4)) * a^2;
190 [~, Kf, Mf] = plate_matrices(dom, 'free', nu);
191 Dhf = diag(1 ./ sqrt(full(diag(Mf)))); Sf = Dhf * full(Kf) * Dhf; Sf = (Sf + Sf') / 2;
192 evf = sort(eig(Sf));
193 out.free_f1 = sqrt(evf(4)) * a^2;
194 `);
195 expect(out.fixed).toBe(out.ic);
196 expect(out.row).toBe(0);
197 expect(out.col).toBe(0);
198 expect(out.lrow).toBe(0);
199 expect(out.mdiag).toBe(1);
200 expect(out.rotx).toBeLessThan(1e-9);
201 expect(out.roty).toBeLessThan(1e-9);
202 expect(out.const).toBeGreaterThan(1e-3);
203 expect(out.zero2).toBeLessThan(1e-9);
204 // the pin removes the translation, so the third eigenvalue is an elastic mode, not a zero
205 expect(out.third).toBeGreaterThan(1e-6);
206 // Rayleigh's interlacing theorem: one constraint lifts each eigenvalue to at most the next one,
207 // so the first elastic mode of the pinned plate (the rim moving in phase, centre held) lies
208 // between the free plate's third zero and its fundamental
209 expect(out.f1).toBeGreaterThan(0);
210 expect(out.f1).toBeLessThan(out.free_f1);
211 // the free plate's first elastic mode twists about both diagonals and is zero at the centre
212 // node, so it is an eigenvector of the pinned plate as well, at the same frequency
213 expect(Math.abs(out.f2 - out.free_f1) / out.free_f1).toBeLessThan(1e-8);
214 });
216 it("holds the centre at zero through both schemes while the edges move, and reports the support", async () => {
217 const spec = testSpec();
218 spec.bc = "pinned-center";
219 spec.implicit = {
220 nx: spec.fd.nx,
221 ny: spec.fd.ny,
222 duration: spec.fd.duration,
223 steps_per_sample: 2,
224 theta: 0.25,
225 };
226 for (const id of ["fd", "implicit"]) {
227 const { result, run } = await solve(spec, id);
228 expect(result.info.stable).toBe(true);
229 expect(result.info.notes.some((n) => /pinned/.test(n))).toBe(true);
230 const { nx, ny, nf, data } = result.frames.w;
231 const jc = Math.floor(nx / 2);
232 const ic = Math.floor(ny / 2);
233 expect(result.info.extra.support_x).toBeCloseTo(result.grid.x[jc], 12);
234 expect(result.info.extra.support_y).toBeCloseTo(result.grid.y[ic], 12);
235 let peak = 0;
236 let edgePeak = 0;
237 for (let k = 0; k < nf; k++) {
238 expect(data[ic + jc * ny + k * nx * ny]).toBe(0);
239 for (let j = 0; j < nx; j++) {
240 edgePeak = Math.max(edgePeak, Math.abs(data[0 + j * ny + k * nx * ny]));
241 for (let i = 0; i < ny; i++)
242 peak = Math.max(peak, Math.abs(data[i + j * ny + k * nx * ny]));
243 }
244 }
245 expect(peak).toBeGreaterThan(0);
246 expect(edgePeak).toBeGreaterThan(0.05 * peak);
247 if (id === "fd") (spec as unknown as { _fd: typeof run })._fd = run;
248 }
249 });
251 it("agrees between the implicit and explicit schemes", async () => {
252 const spec = testSpec();
253 spec.bc = "pinned-center";
254 spec.implicit = {
255 nx: spec.fd.nx,
256 ny: spec.fd.ny,
257 duration: spec.fd.duration,
258 steps_per_sample: 2,
259 theta: 0.25,
260 };
261 const imp = await solve(spec, "implicit");
262 const exp = await solve(spec, "fd");
263 const cmp = await runner.run({
264 bundle: buildCompareBundle(),
265 inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
266 });
267 expect(cmp.ok).toBe(true);
268 if (!cmp.ok) return;
269 expect(parseCompareResult(cmp.value).w_rel_l2).toBeLessThan(0.03);
270 });
272 it("is refused by the modal method", async () => {
273 const spec = testSpec();
274 spec.bc = "pinned-center";
275 await expect(solve(spec, "modal")).rejects.toThrow(/simply supported/);
276 });
277});