1// The design: a MATLAB-syntax script, run by numbl, that fills in a struct
2// named `design` with the material, the bounding box, the outline, the
3// thickness field, the stress-free profile, and the loss constants of the
4// material (from which the decay is derived, loss.ts). The stand is not part
5// of the design: it is a parameter of the setup (params.ts).
6// The outline, the thickness, and the profile are function handles
7// @(x, y, r, theta), vectorized over grid matrices (or plain numbers).
8// compileDesign runs the script and checks the struct; evaluatePlate runs it
9// again with a short trailer that samples the fields on meshgrid matrices of
10// the requested spacing, inside numbl, and reads the matrices back.
11//
12// Lengths are meters throughout. A script is self-contained: nothing is in
13// scope but numbl's own functions and whatever the script defines itself,
14// including ordinary MATLAB functions written at its bottom (numbl hoists
15// them, so the trailer may still follow them).
17import { executeCode, type RuntimeValue, type RuntimeTensor, type RuntimeStruct } from "numbl";
18import { RUNTIME_FILES } from "./mruntime.ts";
19import {
20 type Grid,
21 type Material,
22 type Plate,
23 type Stand,
24 type Thermal,
25 ROOM_TEMPERATURE,
26 insideNode,
27 makeGrid,
28 nodeWeight,
29 nodeX,
30 nodeY,
31 withCellMask,
32} from "./params.ts";
34/** A compiled design: the checked contents of the `design` struct, plus the script that produced it. */
35export interface Design {
36 name: string;
37 material: Material;
38 Lx: number;
39 Ly: number;
40 /** Whether the script gave an outline; without one the plate is the whole bounding box. */
41 hasInside: boolean;
42 /** The thickness when it is a plain number; null for a function. */
43 thicknessConst: number | null;
44 /** The profile when it is a plain number (or absent); null for a function. */
45 profileConst: number | null;
46 /** Temperature of the plate [K], for the thermoelastic damping. */
47 temperature: number;
48 source: string;
49}
51/** Message of an error thrown by the script or by the checks. */
52export class DesignError extends Error {}
54const FILE_NAME = "design.m";
56/**
57 * Runs a design script and returns its workspace variables. `userLines` is
58 * the length of the part of the source the user wrote, so that an error in
59 * the appended trailer is not reported at a line the user cannot see.
60 */
61export function runDesignScript(source: string, userLines = Infinity): Record<string, RuntimeValue> {
62 try {
63 return executeCode(source, { displayResults: false, onOutput: () => {} }, RUNTIME_FILES, FILE_NAME).variableValues;
64 } catch (e) {
65 throw new DesignError(describeError(e, Math.min(source.split("\n").length, userLines)));
66 }
67}
69function describeError(e: unknown, userLines: number): string {
70 if (!(e instanceof Error)) return String(e);
71 const { line, file } = e as { line?: unknown; file?: unknown };
72 const inScript = file === undefined || file === null || file === FILE_NAME;
73 const where = inScript && typeof line === "number" && line >= 1 && line <= userLines ? `line ${line}: ` : "";
74 return `${where}${e.message}`;
75}
77// ---------- Reading runtime values ----------
79function isStruct(v: RuntimeValue | undefined): v is RuntimeStruct {
80 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "struct";
81}
82function isTensor(v: RuntimeValue | undefined): v is RuntimeTensor {
83 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor";
84}
85function isFunction(v: RuntimeValue | undefined): boolean {
86 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "function";
87}
88function isEmpty(v: RuntimeValue | undefined): boolean {
89 return v === undefined || (isTensor(v) && v.data.length === 0);
90}
92/** A scalar number from a runtime value, or null when it is not one. */
93function scalar(v: RuntimeValue | undefined): number | null {
94 if (typeof v === "number") return v;
95 if (typeof v === "boolean") return v ? 1 : 0;
96 if (isTensor(v) && v.data.length === 1 && (!v.imag || v.imag[0] === 0)) return v.data[0];
97 return null;
98}
100function num(v: RuntimeValue | undefined, what: string, min = -Infinity): number {
101 const s = scalar(v);
102 if (s === null || !isFinite(s)) throw new DesignError(`${what} must be a finite number`);
103 if (s < min) throw new DesignError(`${what} must be at least ${min}`);
104 return s;
105}
107function field(s: RuntimeStruct, name: string): RuntimeValue | undefined {
108 return s.fields.get(name);
109}
111function text(v: RuntimeValue | undefined): string | null {
112 if (typeof v === "string") return v;
113 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") return (v as { value: string }).value;
114 return null;
115}
117// ---------- Compilation ----------
119/**
120 * Runs the script and checks the `design` struct it leaves in the
121 * workspace. The field functions are not called here; evaluatePlate samples
122 * them on the grid.
123 */
124export function compileDesign(source: string): Design {
125 const vars = runDesignScript(source);
126 const d = vars.design;
127 if (d === undefined) throw new DesignError("the script must fill in a struct named design (design.material, design.bounds, design.thickness, ...)");
128 if (!isStruct(d)) throw new DesignError("design must be a struct");
130 const nameV = field(d, "name");
131 const name = (nameV !== undefined ? text(nameV) : null)?.trim() || "Untitled design";
133 const m = field(d, "material");
134 if (!isStruct(m)) throw new DesignError("design.material must be a struct with fields E, nu, rho, lossFactor, and the thermal constants alpha, cp, k");
135 const has = (s: RuntimeStruct, k: string) => field(s, k) !== undefined && !isEmpty(field(s, k));
136 const thermalKeys = ["alpha", "cp", "k"];
137 const given = thermalKeys.filter((k) => has(m, k));
138 let thermal: Thermal | null = null;
139 if (given.length === 3) {
140 thermal = { alpha: num(field(m, "alpha"), "design.material.alpha (thermal expansion, 1/K)", 0), cp: num(field(m, "cp"), "design.material.cp (specific heat, J/(kg K))", 1), k: num(field(m, "k"), "design.material.k (thermal conductivity, W/(m K))", 1e-3) };
141 } else if (given.length > 0) {
142 throw new DesignError("design.material needs all three thermal constants for the thermoelastic damping, alpha (1/K), cp (J/(kg K)) and k (W/(m K)), or none of them");
143 }
144 const material: Material = {
145 E: num(field(m, "E"), "design.material.E", 1e6),
146 nu: num(field(m, "nu"), "design.material.nu", 0),
147 rho: num(field(m, "rho"), "design.material.rho", 1),
148 eta: has(m, "lossFactor") ? num(field(m, "lossFactor"), "design.material.lossFactor", 0) : 0,
149 thermal,
150 };
151 if (material.nu >= 0.5) throw new DesignError("design.material.nu must be below 0.5");
152 if (material.eta >= 0.1) throw new DesignError("design.material.lossFactor must be small (metals have loss factors of 1e-4 to 1e-3)");
153 if (field(d, "loss") !== undefined) throw new DesignError("design.loss is no longer used: the decay of each mode follows from design.material.lossFactor, the thermal constants alpha, cp and k, the radiation of the plate, and the felt of the stand");
154 const temperature = has(d, "temperature") ? num(field(d, "temperature"), "design.temperature (K)", 1) : ROOM_TEMPERATURE;
156 const b = field(d, "bounds");
157 if (!isStruct(b)) throw new DesignError("design.bounds must be a struct with fields width, height");
158 const Lx = num(field(b, "width"), "design.bounds.width", 1e-3), Ly = num(field(b, "height"), "design.bounds.height", 1e-3);
160 const inside = field(d, "inside");
161 let hasInside = false;
162 if (!isEmpty(inside)) {
163 if (!isFunction(inside)) throw new DesignError("design.inside must be a function handle @(x, y, r, theta) returning true inside the plate, for example @(x, y, r, theta) r <= D/2");
164 hasInside = true;
165 }
167 const th = field(d, "thickness");
168 if (th === undefined) throw new DesignError("design.thickness is required (a number in meters, or a function handle @(x, y, r, theta))");
169 let thicknessConst: number | null = null;
170 if (!isFunction(th)) {
171 thicknessConst = scalar(th);
172 if (thicknessConst === null || !isFinite(thicknessConst)) throw new DesignError("design.thickness must be a number or a function handle @(x, y, r, theta)");
173 if (thicknessConst <= 0) throw new DesignError("design.thickness must be positive");
174 }
176 const pr = field(d, "profile");
177 let profileConst: number | null = 0;
178 if (pr !== undefined && !isEmpty(pr)) {
179 if (isFunction(pr)) profileConst = null;
180 else {
181 profileConst = scalar(pr);
182 if (profileConst === null || !isFinite(profileConst)) throw new DesignError("design.profile must be a number or a function handle @(x, y, r, theta)");
183 }
184 }
186 if (field(d, "stand") !== undefined) throw new DesignError("design.stand is no longer part of the script: the stand (the felt pads, the preload, the position of the rod) is set in the Stand parameters beside the script and applies to any design");
188 return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, temperature, source };
189}
191// ---------- Evaluation on a grid ----------
193/** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
194export function removePlane(f: Float64Array, g: Grid): void {
195 let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
196 for (let j = 0; j < g.ny; j++) {
197 for (let i = 0; i < g.nx; i++) {
198 const wgt = nodeWeight(g, i, j);
199 if (wgt === 0) continue;
200 const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
201 sw += wgt; sx += wgt * x; sy += wgt * y;
202 sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
203 sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
204 }
205 }
206 if (sw === 0) return;
207 const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
208 if (Math.abs(det) < 1e-300) return;
209 const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
210 const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
211 const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
212 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) f[j * g.nx + i] -= a + b * nodeX(g, i) + c * nodeY(g, j);
213}
215/** The trailer appended to the script to sample the fields on the grid (all names start with cymbalsim_). */
216function trailer(d: Design, g: Grid): string {
217 const fmt = (v: number) => v.toPrecision(17);
218 return [
219 "",
220 "% ---- appended by the simulator: sample the fields on the grid ----",
221 `cymbalsim_nx = ${g.nx}; cymbalsim_ny = ${g.ny};`,
222 `cymbalsim_hx = ${fmt(g.hx)}; cymbalsim_hy = ${fmt(g.hy)};`,
223 `[cymbalsim_X, cymbalsim_Y] = meshgrid(((0:cymbalsim_nx-1) - (cymbalsim_nx-1)/2) * cymbalsim_hx, ((0:cymbalsim_ny-1) - (cymbalsim_ny-1)/2) * cymbalsim_hy);`,
224 "cymbalsim_R = hypot(cymbalsim_X, cymbalsim_Y); cymbalsim_T = atan2(cymbalsim_Y, cymbalsim_X);",
225 "cymbalsim_H = cymbalsim_field(design.thickness, 'thickness', cymbalsim_X, cymbalsim_Y, cymbalsim_R, cymbalsim_T);",
226 d.profileConst === null
227 ? "cymbalsim_Z = cymbalsim_field(design.profile, 'profile', cymbalsim_X, cymbalsim_Y, cymbalsim_R, cymbalsim_T);"
228 : `cymbalsim_Z = ${fmt(d.profileConst)} * ones(cymbalsim_ny, cymbalsim_nx);`,
229 d.hasInside
230 ? [
231 `[cymbalsim_XC, cymbalsim_YC] = meshgrid(((0:cymbalsim_nx-2) - (cymbalsim_nx-2)/2) * cymbalsim_hx, ((0:cymbalsim_ny-2) - (cymbalsim_ny-2)/2) * cymbalsim_hy);`,
232 "cymbalsim_IN = cymbalsim_field(design.inside, 'inside', cymbalsim_XC, cymbalsim_YC, hypot(cymbalsim_XC, cymbalsim_YC), atan2(cymbalsim_YC, cymbalsim_XC));",
233 ].join("\n")
234 : "",
235 "",
236 ].join("\n");
237}
239/** A real [rows, cols] matrix from the workspace as a Float64Array in column-major order. */
240function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number, cols: number, what: string): Float64Array {
241 const v = vars[name];
242 if (!isTensor(v) || v.shape.length !== 2 || v.shape[0] !== rows || v.shape[1] !== cols) throw new DesignError(`${what} did not produce a ${rows} x ${cols} matrix`);
243 if (v.imag && v.imag.some((x) => x !== 0)) throw new DesignError(`${what} returned complex values`);
244 return v.data;
245}
247/**
248 * Samples the design on a grid of about the requested spacing, with the given
249 * stand (from the parameters, not the script): the mask from
250 * the outline at the cell centers, the thickness at the nodes (which must be
251 * positive on the plate), the profile at the nodes with its mean plane
252 * removed. Only the slopes of
253 * the profile enter the physics, so its mean is immaterial, but its tilt is
254 * not: the von Karman strains are not invariant under rotation, and a tilted
255 * stress-free shape would add a spurious membrane coupling.
256 */
257export function evaluatePlate(d: Design, spacing: number, stand: Stand | null = null): Plate {
258 const full = makeGrid(d.Lx, d.Ly, spacing);
259 const { nx, ny } = full;
260 const vars = runDesignScript(d.source + trailer(d, full), d.source.split("\n").length);
261 // meshgrid matrices are [ny, nx] in column-major order: entry (i, j) sits at i * ny + j.
262 const H = matrix(vars, "cymbalsim_H", ny, nx, "thickness");
263 const Z = matrix(vars, "cymbalsim_Z", ny, nx, "profile");
264 let g = full;
265 if (d.hasInside) {
266 const IN = matrix(vars, "cymbalsim_IN", ny - 1, nx - 1, "inside");
267 g = withCellMask(full, (ci, cj) => IN[ci * (ny - 1) + cj] !== 0);
268 }
269 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
270 let any = false;
271 for (let j = 0; j < ny; j++) {
272 for (let i = 0; i < nx; i++) {
273 if (!insideNode(g, i, j)) continue;
274 any = true;
275 const k = j * nx + i, m = i * ny + j;
276 const x = nodeX(g, i), y = nodeY(g, j);
277 const at = `(${x.toFixed(4)}, ${y.toFixed(4)})`;
278 if (!isFinite(H[m])) throw new DesignError(`thickness is ${H[m]} at ${at}`);
279 if (H[m] <= 0) throw new DesignError(`thickness must be positive; it is ${H[m]} at ${at}`);
280 if (!isFinite(Z[m])) throw new DesignError(`profile is ${Z[m]} at ${at}`);
281 h[k] = H[m];
282 w0[k] = Z[m];
283 }
284 }
285 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
286 removePlane(w0, g);
287 for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) if (!insideNode(g, i, j)) w0[j * nx + i] = 0;
288 return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, stand, temperature: d.temperature };
289}
291/** Whether the plate is flat (no profile), in which case bending and in-plane motion decouple in the linear model. */
292export function isFlat(pl: Plate): boolean {
293 for (let k = 0; k < pl.w0.length; k++) if (pl.w0[k] !== 0) return false;
294 return true;
295}
297/** Peak-to-peak height and the largest slope of the profile over the plate nodes. */
298export function profileStats(pl: Plate): { peakToPeak: number; maxSlope: number } {
299 const g = pl.grid, w0 = pl.w0;
300 let lo = Infinity, hi = -Infinity, maxSlope = 0;
301 const at = (i: number, j: number) => w0[j * g.nx + i];
302 const slope = (i: number, j: number, di: number, dj: number, hg: number) => {
303 const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
304 if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * hg);
305 if (hasP) return (at(i + di, j + dj) - at(i, j)) / hg;
306 if (hasM) return (at(i, j) - at(i - di, j - dj)) / hg;
307 return 0;
308 };
309 for (let j = 0; j < g.ny; j++) {
310 for (let i = 0; i < g.nx; i++) {
311 if (!insideNode(g, i, j)) continue;
312 const v = at(i, j);
313 lo = Math.min(lo, v); hi = Math.max(hi, v);
314 maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
315 }
316 }
317 return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
318}