96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 1// Web Worker that renders the triangle strike by executing the MATLAB
2// implementation in ../matlab through numbl. The stages run as separate
3// executeCode calls so intermediates can be cached and cheap parameter
4// changes skip the expensive ones:
5// solve tier -> triangle_basis (FEM assembly + eigensolve; the expensive
6// stage, cached on geometry/material/damping params)
7// radiate tier -> triangle_strikes + triangle_radiate (cached on the
8// strike process, the listener, and the synthesis rate)
9// output tier -> triangle_output (always cheap, always re-run)
10// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
12import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
13import { hashMFiles, type MFile } from '../mfiles';
14import { SOLVE_KEYS, RADIATE_KEYS, type TriangleParams } from '../params';
16export interface RenderRequest {
17 type: 'render';
18 id: number;
19 params: TriangleParams;
20 quality: 'preview' | 'full';
21 mfiles: MFile[];
22}
24export interface RenderResponse {
25 type: 'result';
26 id: number;
27 quality: 'preview' | 'full';
28 y: Float32Array; // final signal, [-1, 1]
29 sampleRate: number;
30 strikeT: Float32Array; // strike times, s
31 strikeV: Float32Array; // strike velocities, m/s
32 strikeS: Float32Array; // strike arclength fractions
33 modeF: Float32Array; // mode frequencies, Hz
34 modeZeta: Float32Array; // damping ratio of each mode
35 modeZetaRad: Float32Array; // radiation part of the damping
36 modeFamily: Float32Array; // 0 = in-plane, 1 = out-of-plane
37 nodeX: Float32Array; // rod centerline nodes, m
38 nodeY: Float32Array;
39 nodeS: Float32Array; // arclength fraction of each node
40 shapeUx: Float32Array; // node displacements per mode, nn x n col-major
41 shapeUy: Float32Array;
42 shapeUz: Float32Array;
43 diag: {
44 modesKept: number;
45 modesIn: number;
46 modesOut: number;
47 strikeCount: number;
48 meanV: number;
49 maxV: number;
50 meanTauUs: number;
51 peakForceN: number;
52 durationS: number;
53 peakPa: number;
54 peakSplDb: number;
55 directPeakPa: number;
56 listenerDistM: number;
57 warnings: string[];
58 solveMs: number; // 0 when the cached eigensolve was reused
59 radiateMs: number; // 0 when the cached radiation was reused
60 solveCached: boolean;
61 radiateCached: boolean;
62 oversample: number;
63 };
64}
66export interface ProgressMessage {
67 type: 'progress';
68 id: number;
69 stage: 'solve' | 'radiate' | 'output';
70 frac: number;
71}
73// Raised by a fault in the .m code (the user can edit it), so it is reported
74// rather than left to stall the request.
75export interface ErrorMessage {
76 type: 'error';
77 id: number;
78 message: string;
79 file?: string;
80 line?: number;
81}
83// Basis fields that cross from the solve stage into the radiate stage as
84// numbl RuntimeValues, plus their Float32 copies for the plots.
85interface SolveEntry {
86 vars: Record<string, RuntimeValue>;
87 modeF: Float32Array;
88 modeZeta: Float32Array;
89 modeZetaRad: Float32Array;
90 modeFamily: Float32Array;
91 nodeX: Float32Array;
92 nodeY: Float32Array;
93 nodeS: Float32Array;
94 shapeUx: Float32Array;
95 shapeUy: Float32Array;
96 shapeUz: Float32Array;
97 modesIn: number;
98 modesOut: number;
99}
101interface RadiateEntry {
102 p: RuntimeValue;
103 strikeT: Float32Array;
104 strikeV: Float32Array;
105 strikeS: Float32Array;
106 strikeCount: number;
107 meanV: number;
108 maxV: number;
109 meanTauUs: number;
110 peakForceN: number;
111 durationS: number;
112 directPeakPa: number;
113 listenerDistM: number;
114 listenerInside: boolean;
115 fsHi: number;
116 oversample: number;
117}
119// Small LRU caches. The solve cache does not depend on the oversample rate,
120// so a preview render right after a full one reuses the same eigensolve.
121const solveCache = new Map<string, SolveEntry>();
122const radiateCache = new Map<string, RadiateEntry>();
123const SOLVE_CACHE_MAX = 2;
124const RADIATE_CACHE_MAX = 3;
126function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
127 const v = cache.get(key);
128 if (v !== undefined) {
129 cache.delete(key);
130 cache.set(key, v);
131 }
132 return v;
133}
135function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
136 cache.delete(key);
137 cache.set(key, value);
138 while (cache.size > max) {
139 const oldest = cache.keys().next().value as string;
140 cache.delete(oldest);
141 }
142}
144function keyFor(params: TriangleParams, keys: readonly (keyof TriangleParams)[]): string {
145 return keys.map((k) => String(params[k])).join('|');
146}
148// MATLAB source lines assigning every parameter onto a struct P.
149function pLiteral(params: TriangleParams): string {
150 const lines: string[] = [];
151 for (const [k, v] of Object.entries(params)) {
152 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
153 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
154 }
155 return lines.join('\n');
156}
158function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
159 return executeCode(
160 source,
161 { initialVariableValues: inputs, displayResults: false },
162 mfiles,
163 'stage.m',
164 );
165}
167function asTensor(v: RuntimeValue): RuntimeTensor {
168 const t = v as RuntimeTensor;
169 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
170 return t;
171}
173function asScalar(v: RuntimeValue): number {
174 if (typeof v === 'number') return v;
175 const t = v as RuntimeTensor;
176 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
177 if (typeof v === 'boolean') return v ? 1 : 0;
178 throw new Error('expected scalar result from numbl stage');
179}
181function toF32(v: RuntimeValue): Float32Array {
182 const t = asTensor(v);
183 const out = new Float32Array(t.data.length);
184 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
185 return out;
186}
188// The basis fields handed from the solve stage into the radiate stage, in the
189// order the reassembly script expects.
190const BASIS_FIELDS = [
191 'f', 'omega', 'zeta', 'zetaRad', 'family',
192 'nodeX', 'nodeY', 'sfrac', 'totalLen',
193 'NUx', 'NUy', 'NUz',
194 'elemX', 'elemY', 'elemTanX', 'elemTanY', 'elemLen',
195 'EUpx', 'EUpy', 'EUpz',
196];
198const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
200self.onmessage = (ev: MessageEvent<RenderRequest>) => {
201 const msg = ev.data;
202 if (msg.type !== 'render') return;
204 // numbl reports a file it could not parse through console.warn and then
205 // skips it, so the failure surfaces at the call site as "Undefined
206 // function". Capture the warning to report the actual syntax error instead.
207 const warnings: string[] = [];
208 const origWarn = console.warn;
209 console.warn = (...args: unknown[]) => {
210 warnings.push(args.map(String).join(' '));
211 };
213 try {
214 render(msg);
215 } catch (err) {
216 // Most likely a fault in the (user-editable) .m code. Report it so the
217 // request completes and the message reaches the code panel.
218 const e = err as { message?: string; file?: string; line?: number };
219 let m: ErrorMessage = {
220 type: 'error',
221 id: msg.id,
222 message: e?.message ? String(e.message) : String(err),
223 file: typeof e?.file === 'string' ? e.file : undefined,
224 line: typeof e?.line === 'number' ? e.line : undefined,
225 };
226 for (const w of warnings) {
227 const hit = SKIP_WARNING.exec(w);
228 if (hit) {
229 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
230 break;
231 }
232 }
233 (self as unknown as Worker).postMessage(m);
234 } finally {
235 console.warn = origWarn;
236 }
237};
239function render(msg: RenderRequest): void {
240 const { id, params, mfiles } = msg;
241 let quality = msg.quality;
243 // The .m sources are part of the cache identity: editing the code must
244 // invalidate everything computed from the previous version.
245 const codeKey = hashMFiles(mfiles);
246 const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
248 const fullOversample = params.oversample;
249 const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
250 if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
251 quality = 'full'; // full-rate radiation already cached
252 }
253 const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
254 const fsHi = params.sample_rate * oversample;
256 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
258 // ---- solve stage: FEM assembly and eigensolve ----
259 let solveMs = 0;
260 let solveCached = true;
261 let solveEntry = lruGet(solveCache, solveKey);
262 if (!solveEntry) {
263 solveCached = false;
264 post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
265 const t0 = performance.now();
266 const r = runStage(mfiles, `
267${pLiteral(params)}
268basis = triangle_basis(P);
269${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
270modesIn = sum(basis.family == 0);
271modesOut = sum(basis.family == 1);
272`);
273 solveMs = performance.now() - t0;
275 const vars: Record<string, RuntimeValue> = {};
276 for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
277 solveEntry = {
278 vars,
279 modeF: toF32(r.variableValues.b_f),
280 modeZeta: toF32(r.variableValues.b_zeta),
281 modeZetaRad: toF32(r.variableValues.b_zetaRad),
282 modeFamily: toF32(r.variableValues.b_family),
283 nodeX: toF32(r.variableValues.b_nodeX),
284 nodeY: toF32(r.variableValues.b_nodeY),
285 nodeS: toF32(r.variableValues.b_sfrac),
286 shapeUx: toF32(r.variableValues.b_NUx),
287 shapeUy: toF32(r.variableValues.b_NUy),
288 shapeUz: toF32(r.variableValues.b_NUz),
289 modesIn: asScalar(r.variableValues.modesIn),
290 modesOut: asScalar(r.variableValues.modesOut),
291 };
292 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
293 }
295 // ---- radiate stage: strike process + modal synthesis ----
296 const radiateKey = radiateParamsKey + `|os=${oversample}`;
297 let radiateMs = 0;
298 let radiateCached = true;
299 let radEntry = lruGet(radiateCache, radiateKey);
300 if (!radEntry) {
301 radiateCached = false;
302 post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
303 const t0 = performance.now();
304 const r = runStage(mfiles, `
305${pLiteral(params)}
306${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
307basis.n = length(basis.f);
308[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
309[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, ${fsHi});
310strikeCount = sinfo.count;
311meanV = sinfo.meanV;
312maxV = sinfo.maxV;
313meanTauUs = sinfo.meanTauUs;
314peakForceN = sinfo.peakForceN;
315durationS = sinfo.durationS;
316directPeakPa = radinfo.directPeakPa;
317listenerDistM = radinfo.listenerDistM;
318listenerInside = radinfo.listenerInside;
319`, solveEntry.vars);
320 radiateMs = performance.now() - t0;
321 radEntry = {
322 p: r.variableValues.p,
323 strikeT: toF32(r.variableValues.tstr),
324 strikeV: toF32(r.variableValues.vstr),
325 strikeS: toF32(r.variableValues.sstr),
326 strikeCount: asScalar(r.variableValues.strikeCount),
327 meanV: asScalar(r.variableValues.meanV),
328 maxV: asScalar(r.variableValues.maxV),
329 meanTauUs: asScalar(r.variableValues.meanTauUs),
330 peakForceN: asScalar(r.variableValues.peakForceN),
331 durationS: asScalar(r.variableValues.durationS),
332 directPeakPa: asScalar(r.variableValues.directPeakPa),
333 listenerDistM: asScalar(r.variableValues.listenerDistM),
334 listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
335 fsHi,
336 oversample,
337 };
338 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
339 }
341 // ---- output stage ----
342 post({ type: 'progress', id, stage: 'output', frac: 0.9 });
343 const rOut = runStage(mfiles, `
344${pLiteral(params)}
345[y, peakPa] = triangle_output(P, p, ${radEntry.oversample});
346`, { p: radEntry.p });
348 const y32 = toF32(rOut.variableValues.y);
349 const peakPa = asScalar(rOut.variableValues.peakPa);
351 const warnings: string[] = [];
352 if (radEntry.listenerInside) {
353 warnings.push(
354 'the listener is within 2 cm of the rod; the dipole field is evaluated at a clamped minimum distance',
355 );
356 }
357 if (!params.normalize) {
358 let peak = 0;
359 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
360 if (peak > 1.0) {
361 warnings.push(
362 `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
363 );
364 }
365 }
366 const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
367 if (topMode < 0.4 * params.max_mode_hz) {
368 warnings.push(
369 `the highest retained mode is only ${topMode.toFixed(0)} Hz; the shimmer above it falls entirely on the direct click`,
370 );
371 }
373 const resp: RenderResponse = {
374 type: 'result',
375 id,
376 quality,
377 y: y32,
378 sampleRate: params.sample_rate,
379 strikeT: radEntry.strikeT.slice(),
380 strikeV: radEntry.strikeV.slice(),
381 strikeS: radEntry.strikeS.slice(),
382 modeF: solveEntry.modeF.slice(),
383 modeZeta: solveEntry.modeZeta.slice(),
384 modeZetaRad: solveEntry.modeZetaRad.slice(),
385 modeFamily: solveEntry.modeFamily.slice(),
386 nodeX: solveEntry.nodeX.slice(),
387 nodeY: solveEntry.nodeY.slice(),
388 nodeS: solveEntry.nodeS.slice(),
389 shapeUx: solveEntry.shapeUx.slice(),
390 shapeUy: solveEntry.shapeUy.slice(),
391 shapeUz: solveEntry.shapeUz.slice(),
392 diag: {
393 modesKept: solveEntry.modeF.length,
394 modesIn: solveEntry.modesIn,
395 modesOut: solveEntry.modesOut,
396 strikeCount: radEntry.strikeCount,
397 meanV: radEntry.meanV,
398 maxV: radEntry.maxV,
399 meanTauUs: radEntry.meanTauUs,
400 peakForceN: radEntry.peakForceN,
401 durationS: radEntry.durationS,
402 peakPa,
403 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
404 directPeakPa: radEntry.directPeakPa,
405 listenerDistM: radEntry.listenerDistM,
406 warnings,
407 solveMs,
408 radiateMs,
409 solveCached,
410 radiateCached,
411 oversample: radEntry.oversample,
412 },
413 };
414 (self as unknown as Worker).postMessage(resp, [
415 resp.y.buffer,
416 resp.strikeT.buffer,
417 resp.strikeV.buffer,
418 resp.strikeS.buffer,
419 resp.modeF.buffer,
420 resp.modeZeta.buffer,
421 resp.modeZetaRad.buffer,
422 resp.modeFamily.buffer,
423 resp.nodeX.buffer,
424 resp.nodeY.buffer,
425 resp.nodeS.buffer,
426 resp.shapeUx.buffer,
427 resp.shapeUy.buffer,
428 resp.shapeUz.buffer,
429 ]);
430}