2004981Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1import { readFileSync } from "node:fs";
2import { fileURLToPath } from "node:url";
3import path from "node:path";
4import { NodeRunner } from "../src/node";
5import {
6 LIBRARY_FILES,
7 buildModesBundle,
8 buildRunBundle,
9 buildStrikeBundle,
10 cloneSpec,
11 generateParamsBlock,
12 isPlainStruct,
13 methodById,
14 parsePlateResult,
15 scenarioPreset,
16 type MethodId,
17 type PlainStruct,
18 type PlainValue,
19 type PlateResult,
20 type ProblemSpec,
21 type RunResult,
22} from "../src/index";
24export const runner = new NodeRunner();
26/**
27 * The reference problem: the centre hit of the two Python simulators, with a
28 * short render (the reference WAV excerpts are 50 ms) and the finite-element
29 * mesh set per test. The strike point and the mesh are what the tests vary.
30 */
31export function testSpec(scenario = "center-hit"): ProblemSpec {
32 const s = cloneSpec(scenarioPreset(scenario).spec);
33 s.render.duration = 0.1;
34 s.render.nframes = 6;
35 s.render.frame_span = 0.001;
36 s.fem.mesh_x = 16;
37 s.fem.mesh_y = 16;
38 return s;
39}
41/** Run a method as its two stages and return the typed result with the raw modes. */
42export async function solve(
43 spec: ProblemSpec,
44 method: MethodId,
45): Promise<{ modes: PlainValue; result: PlateResult; modesRun: RunResult; strikeRun: RunResult }> {
46 const m = methodById(method);
47 const modesRun = await runner.run({ bundle: buildModesBundle(spec, m) });
48 if (!modesRun.ok)
49 throw new Error(`${method} modes failed: ${modesRun.error}\n${modesRun.output.join("")}`);
50 const strikeRun = await runner.run({
51 bundle: buildStrikeBundle(spec, m),
52 inputs: { modes: modesRun.value },
53 });
54 if (!strikeRun.ok)
55 throw new Error(`${method} strike failed: ${strikeRun.error}\n${strikeRun.output.join("")}`);
56 return { modes: modesRun.value, result: parsePlateResult(strikeRun.value), modesRun, strikeRun };
57}
59/** Run a method as one script (the generated main.m). */
60export async function solveSingle(spec: ProblemSpec, method: MethodId): Promise<PlateResult> {
61 const run = await runner.run({ bundle: buildRunBundle(spec, methodById(method)) });
62 if (!run.ok) throw new Error(`${method} failed: ${run.error}\n${run.output.join("")}`);
63 return parsePlateResult(run.value);
64}
66/** Run a custom script after the params block and return the struct it leaves in `out`. */
67export async function runScript(spec: ProblemSpec, body: string): Promise<PlainStruct> {
68 const main = [generateParamsBlock(spec), "", body, ""].join("\n");
69 const run = await runner.run({
70 bundle: {
71 files: [{ name: "script.m", source: main }, ...LIBRARY_FILES],
72 entry: "script.m",
73 resultVariable: "out",
74 },
75 });
76 if (!run.ok) throw new Error(`script failed: ${run.error}\n${run.output.join("")}`);
77 if (!isPlainStruct(run.value)) throw new Error("script did not leave a struct in out");
78 return run.value;
79}
81export function scalar(s: PlainStruct, name: string): number {
82 const v: PlainValue | undefined = s.fields[name];
83 if (typeof v === "number") return v;
84 if (typeof v === "boolean") return v ? 1 : 0;
85 if (v && typeof v === "object" && v.kind === "tensor" && v.data.length >= 1) return v.data[0];
86 throw new Error(`${name} is not a scalar`);
87}
89export function vector(s: PlainStruct, name: string): Float64Array {
90 const v: PlainValue | undefined = s.fields[name];
91 if (typeof v === "number") return Float64Array.of(v);
92 if (v && typeof v === "object" && v.kind === "tensor") return v.data;
93 throw new Error(`${name} is not a vector`);
94}
96export function relL2(a: ArrayLike<number>, b: ArrayLike<number>): number {
97 let num = 0;
98 let den = 0;
99 for (let i = 0; i < a.length; i++) {
100 num += (a[i] - b[i]) ** 2;
101 den += a[i] ** 2;
102 }
103 return Math.sqrt(num / Math.max(den, 1e-300));
104}
106export function maxAbs(a: ArrayLike<number>): number {
107 let m = 0;
108 for (let i = 0; i < a.length; i++) m = Math.max(m, Math.abs(a[i]));
109 return m;
110}
112// ---- the reference diagnostics of the Python simulators (sim-docs)
114export interface ReferenceMode {
115 mode: number;
116 frequency_hz: number;
117 q_total: number;
118 hit_participation: number;
119 radiation_participation: number;
120 support_participation?: number;
121}
123export interface Reference {
124 model: string;
125 derived: Record<string, number>;
126 modes: ReferenceMode[];
127}
129const fixtures = path.join(path.dirname(fileURLToPath(import.meta.url)), "fixtures");
131export function reference(name: string): Reference {
132 return JSON.parse(readFileSync(path.join(fixtures, `${name}.json`), "utf8")) as Reference;
133}
135/** The first 50 ms of the reference WAV, 24-bit samples as fractions of full scale. */
136export function referenceWav(name: string): { fs: number; samples: number[] } {
137 return JSON.parse(readFileSync(path.join(fixtures, `${name}_wav.json`), "utf8")) as {
138 fs: number;
139 samples: number[];
140 };
141}
143/** Frequencies of the audible modes of a result, in order. */
144export function audibleFreqs(r: PlateResult): number[] {
145 return Array.from(r.modes.freq).filter((_, k) => r.modes.audible[k]);
146}
148/** Indices of the audible modes of a result whose neighbours are at least `gap` Hz away (no degeneracy). */
149export function isolatedAudible(r: PlateResult, gap = 5): number[] {
150 const f = r.modes.freq;
151 const out: number[] = [];
152 for (let k = 0; k < f.length; k++) {
153 if (!r.modes.audible[k]) continue;
154 const near =
155 (k > 0 && Math.abs(f[k] - f[k - 1]) < gap) ||
156 (k + 1 < f.length && Math.abs(f[k + 1] - f[k]) < gap);
157 if (!near) out.push(k);
158 }
159 return out;
160}
162export { cloneSpec };