1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1import test from "node:test";
2import assert from "node:assert/strict";
3import { dktStiffness } from "../src/model/dkt.ts";
4import { symmetricEigenvalues, planeResiduals } from "../src/model/linalg.ts";
5import { fixture, maxAbs, maxAbsDiff } from "./helpers.ts";
7test("DKT stiffness matches the Python (autodiff) element", () => {
8 const f = fixture.dkt;
9 const Ke = dktStiffness(f.x, f.y, f.D, f.nu);
10 assert.ok(maxAbsDiff(Ke, f.Ke) < 1e-9 * maxAbs(f.Ke));
11});
13test("DKT rigid modes, constant curvature and twist", () => {
14 const x = [0, 0.03, 0.01], y = [0, 0.005, 0.025], D = 10, nu = 0.3;
15 const Ke = dktStiffness(x, y, D, nu);
16 const apply = (U: number[]) => Array.from({ length: 9 }, (_, i) => U.reduce((s, u, j) => s + Ke[9 * i + j] * u, 0));
17 const energy = (U: number[]) => 0.5 * U.reduce((s, u, i) => s + u * apply(U)[i], 0);
18 const A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
19 const scale = maxAbs(Ke);
20 for (const U of [[1, 0, 0, 1, 0, 0, 1, 0, 0], [y[0], 1, 0, y[1], 1, 0, y[2], 1, 0], [-x[0], 0, 1, -x[1], 0, 1, -x[2], 0, 1]])
21 assert.ok(maxAbs(apply(U)) < 1e-12 * scale);
22 const Ucur = x.flatMap((xi) => [0.5 * xi * xi, 0, -xi]);
23 assert.ok(Math.abs(energy(Ucur) / (0.5 * D * A) - 1) < 1e-9);
24 const Utw = x.map((xi, i) => [xi * y[i], xi, -y[i]]).flat();
25 assert.ok(Math.abs(energy(Utw) / (D * (1 - nu) * A) - 1) < 1e-9);
26 const ev = symmetricEigenvalues(Ke, 9);
27 assert.ok(ev[2] < 1e-10 * ev[8] && ev[3] > 1e-6 * ev[8]);
28});
30test("Jacobi eigenvalues and plane fit", () => {
31 const ev = symmetricEigenvalues([2, 1, 0, 1, 2, 0, 0, 0, 5], 3);
32 assert.ok(maxAbsDiff(ev, [1, 3, 5]) < 1e-12);
33 const res = planeResiduals([0, 1, 0, 1], [0, 0, 1, 1], [1, 3, 4, 6]);
34 assert.ok(maxAbs(res) < 1e-12);
6d7f34aRank-safe plane fit for contact patchesJeremy Magland 35 // collinear points (a stick's line): the fit degenerates to a line, as numpy's lstsq gives
36 const line = planeResiduals([-1, 0, 1], [-2, 0, 2], [1, 2, 3.5]);
37 assert.ok(maxAbsDiff(line, [1 / 12, -1 / 6, 1 / 12]) < 1e-12);
38 // a single point: zero residual
39 assert.ok(maxAbs(planeResiduals([0.3], [0.1], [0.02])) < 1e-15);