1// The cymbal design: a MATLAB-syntax script, run by numbl, that fills in a
2// struct named `design` with the material, the bounding box, the outline,
3// the thickness field, the stress-free profile, the mount, and the losses.
4// The outline, the thickness, and the profile are function handles
5// @(x, y, r, theta), vectorized over grid matrices (or plain numbers).
6// compileDesign runs the script and checks the struct; evaluatePlate runs it
7// again with a short trailer that samples the fields on meshgrid matrices of
8// the requested spacing, inside numbl, and reads the matrices back.
9//
10// Lengths are meters throughout; `inch` and `mm` are helper functions. The
11// helpers a script may call are listed in mhelpers.ts.
13import { executeCode, type RuntimeValue, type RuntimeTensor, type RuntimeStruct } from "numbl";
14import { HELPER_FILES } from "./mhelpers.ts";
15import {
16 type Grid,
17 type Loss,
18 type Material,
19 type Mount,
20 type Plate,
21 insideNode,
22 makeGrid,
23 nodeWeight,
24 nodeX,
25 nodeY,
26 withCellMask,
27} from "./params.ts";
29/** A compiled design: the checked contents of the `design` struct, plus the script that produced it. */
30export interface Design {
31 name: string;
32 material: Material;
33 Lx: number;
34 Ly: number;
35 /** Whether the script gave an outline; without one the plate is the whole bounding box. */
36 hasInside: boolean;
37 /** The thickness when it is a plain number; null for a function. */
38 thicknessConst: number | null;
39 /** The profile when it is a plain number (or absent); null for a function. */
40 profileConst: number | null;
41 mount: Mount | null;
42 loss: Loss;
43 source: string;
44}
46export const DEFAULT_LOSS: Loss = { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 };
48/** Message of an error thrown by the script or by the checks. */
49export class DesignError extends Error {}
51const FILE_NAME = "design.m";
53/** Runs a design script (with the helpers in scope) and returns its workspace variables. */
54export function runDesignScript(source: string): Record<string, RuntimeValue> {
55 try {
56 return executeCode(source, { displayResults: false, onOutput: () => {} }, HELPER_FILES, FILE_NAME).variableValues;
57 } catch (e) {
58 throw new DesignError(describeError(e, source));
59 }
60}
62/** The message of a numbl error (a RuntimeError or a parser SyntaxError), with its line when it lies in the user's script. */
63function describeError(e: unknown, source: string): string {
64 if (!(e instanceof Error)) return String(e);
65 const { line, file } = e as { line?: unknown; file?: unknown };
66 const lines = source.split("\n").length;
67 // A line is only meaningful when the error arose in the script itself, not in a helper .m file.
68 const inScript = file === undefined || file === null || file === FILE_NAME;
69 const where = inScript && typeof line === "number" && line >= 1 && line <= lines ? `line ${line}: ` : "";
70 return `${where}${e.message}`;
71}
73// ---------- Reading runtime values ----------
75function isStruct(v: RuntimeValue | undefined): v is RuntimeStruct {
76 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "struct";
77}
78function isTensor(v: RuntimeValue | undefined): v is RuntimeTensor {
79 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor";
80}
81function isFunction(v: RuntimeValue | undefined): boolean {
82 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "function";
83}
84function isEmpty(v: RuntimeValue | undefined): boolean {
85 return v === undefined || (isTensor(v) && v.data.length === 0);
86}
88/** A scalar number from a runtime value, or null when it is not one. */
89function scalar(v: RuntimeValue | undefined): number | null {
90 if (typeof v === "number") return v;
91 if (typeof v === "boolean") return v ? 1 : 0;
92 if (isTensor(v) && v.data.length === 1 && (!v.imag || v.imag[0] === 0)) return v.data[0];
93 return null;
94}
96function num(v: RuntimeValue | undefined, what: string, min = -Infinity): number {
97 const s = scalar(v);
98 if (s === null || !isFinite(s)) throw new DesignError(`${what} must be a finite number`);
99 if (s < min) throw new DesignError(`${what} must be at least ${min}`);
100 return s;
101}
103function field(s: RuntimeStruct, name: string): RuntimeValue | undefined {
104 return s.fields.get(name);
105}
107function text(v: RuntimeValue | undefined): string | null {
108 if (typeof v === "string") return v;
109 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") return (v as { value: string }).value;
110 return null;
111}
113// ---------- Compilation ----------
115/**
116 * Runs the script and checks the `design` struct it leaves in the
117 * workspace. The field functions are not called here; evaluatePlate samples
118 * them on the grid.
119 */
120export function compileDesign(source: string): Design {
121 const vars = runDesignScript(source);
122 const d = vars.design;
123 if (d === undefined) throw new DesignError("the script must fill in a struct named design (design.material, design.bounds, design.thickness, ...)");
124 if (!isStruct(d)) throw new DesignError("design must be a struct");
126 const nameV = field(d, "name");
127 const name = (nameV !== undefined ? text(nameV) : null)?.trim() || "Untitled design";
129 const m = field(d, "material");
130 if (!isStruct(m)) throw new DesignError("design.material must be a struct with fields E, nu, rho");
131 const material: Material = { E: num(field(m, "E"), "design.material.E", 1e6), nu: num(field(m, "nu"), "design.material.nu", 0), rho: num(field(m, "rho"), "design.material.rho", 1) };
132 if (material.nu >= 0.5) throw new DesignError("design.material.nu must be below 0.5");
134 const b = field(d, "bounds");
135 if (!isStruct(b)) throw new DesignError("design.bounds must be a struct with fields width, height");
136 const Lx = num(field(b, "width"), "design.bounds.width", 1e-3), Ly = num(field(b, "height"), "design.bounds.height", 1e-3);
138 const inside = field(d, "inside");
139 let hasInside = false;
140 if (!isEmpty(inside)) {
141 if (!isFunction(inside)) throw new DesignError("design.inside must be a function handle @(x, y, r, theta) returning true inside the plate, for example circle(D, hole)");
142 hasInside = true;
143 }
145 const th = field(d, "thickness");
146 if (th === undefined) throw new DesignError("design.thickness is required (a number in meters, or a function handle @(x, y, r, theta))");
147 let thicknessConst: number | null = null;
148 if (!isFunction(th)) {
149 thicknessConst = scalar(th);
150 if (thicknessConst === null || !isFinite(thicknessConst)) throw new DesignError("design.thickness must be a number or a function handle @(x, y, r, theta)");
151 if (thicknessConst <= 0) throw new DesignError("design.thickness must be positive");
152 }
154 const pr = field(d, "profile");
155 let profileConst: number | null = 0;
156 if (pr !== undefined && !isEmpty(pr)) {
157 if (isFunction(pr)) profileConst = null;
158 else {
159 profileConst = scalar(pr);
160 if (profileConst === null || !isFinite(profileConst)) throw new DesignError("design.profile must be a number or a function handle @(x, y, r, theta)");
161 }
162 }
164 let mount: Mount | null = null;
165 const mo = field(d, "mount");
166 if (mo !== undefined && !isEmpty(mo) && mo !== false) {
167 if (!isStruct(mo)) throw new DesignError("design.mount must be a struct with fields radius, stiffness, damping (and optionally x, y), or [] for a free plate");
168 const opt = (k: string, what: string, dflt: number, min = -Infinity) => (field(mo, k) === undefined ? dflt : num(field(mo, k), what, min));
169 mount = {
170 x: opt("x", "design.mount.x", 0),
171 y: opt("y", "design.mount.y", 0),
172 radius: num(field(mo, "radius"), "design.mount.radius", 0),
173 stiffness: num(field(mo, "stiffness"), "design.mount.stiffness", 0),
174 damping: opt("damping", "design.mount.damping", 0, 0),
175 };
176 }
178 const lo = field(d, "loss");
179 const loss: Loss = { ...DEFAULT_LOSS };
180 if (lo !== undefined && !isEmpty(lo)) {
181 if (!isStruct(lo)) throw new DesignError("design.loss must be a struct with fields sigma0, sigma1, sigmaIn");
182 for (const k of ["sigma0", "sigma1", "sigmaIn"] as const) if (field(lo, k) !== undefined) loss[k] = num(field(lo, k), `design.loss.${k}`, 0);
183 }
185 return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, mount, loss, source };
186}
188// ---------- Evaluation on a grid ----------
190/** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
191export function removePlane(f: Float64Array, g: Grid): void {
192 let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
193 for (let j = 0; j < g.ny; j++) {
194 for (let i = 0; i < g.nx; i++) {
195 const wgt = nodeWeight(g, i, j);
196 if (wgt === 0) continue;
197 const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
198 sw += wgt; sx += wgt * x; sy += wgt * y;
199 sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
200 sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
201 }
202 }
203 if (sw === 0) return;
204 // Normal equations of the 3-parameter fit, solved by Cramer's rule.
205 const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
206 if (Math.abs(det) < 1e-300) return;
207 const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
208 const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
209 const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
210 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);
211}
213/** The trailer appended to the script to sample the fields on the grid (all names start with steelsim_). */
214function trailer(d: Design, g: Grid): string {
215 const fmt = (v: number) => v.toPrecision(17);
216 return [
217 "",
218 "% ---- appended by the simulator: sample the fields on the grid ----",
219 `steelsim_nx = ${g.nx}; steelsim_ny = ${g.ny};`,
220 `steelsim_hx = ${fmt(g.hx)}; steelsim_hy = ${fmt(g.hy)};`,
221 `[steelsim_X, steelsim_Y] = meshgrid(((0:steelsim_nx-1) - (steelsim_nx-1)/2) * steelsim_hx, ((0:steelsim_ny-1) - (steelsim_ny-1)/2) * steelsim_hy);`,
222 "steelsim_R = hypot(steelsim_X, steelsim_Y); steelsim_T = atan2(steelsim_Y, steelsim_X);",
223 "steelsim_H = steelsim_field(design.thickness, 'thickness', steelsim_X, steelsim_Y, steelsim_R, steelsim_T);",
224 d.profileConst === null
225 ? "steelsim_Z = steelsim_field(design.profile, 'profile', steelsim_X, steelsim_Y, steelsim_R, steelsim_T);"
226 : `steelsim_Z = ${fmt(d.profileConst)} * ones(steelsim_ny, steelsim_nx);`,
227 d.hasInside
228 ? [
229 `[steelsim_XC, steelsim_YC] = meshgrid(((0:steelsim_nx-2) - (steelsim_nx-2)/2) * steelsim_hx, ((0:steelsim_ny-2) - (steelsim_ny-2)/2) * steelsim_hy);`,
230 "steelsim_IN = steelsim_field(design.inside, 'inside', steelsim_XC, steelsim_YC, hypot(steelsim_XC, steelsim_YC), atan2(steelsim_YC, steelsim_XC));",
231 ].join("\n")
232 : "",
233 "",
234 ].join("\n");
235}
237/** A real [rows, cols] matrix from the workspace as a Float64Array in column-major order. */
238function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number, cols: number, what: string): Float64Array {
239 const v = vars[name];
240 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`);
241 if (v.imag && v.imag.some((x) => x !== 0)) throw new DesignError(`${what} returned complex values`);
242 return v.data;
243}
245/**
246 * Samples the design on a grid of about the requested spacing: the mask from
247 * the outline at the cell centers, the thickness at the nodes (which must be
248 * positive on the plate), and the profile at the nodes with its mean plane
249 * removed. Only the slopes of the profile enter the physics, so its mean is
250 * immaterial, but its tilt is not: the von Karman strains are not invariant
251 * under rotation, and a tilted stress-free shape would add a spurious
252 * membrane coupling.
253 */
254export function evaluatePlate(d: Design, spacing: number): Plate {
255 const full = makeGrid(d.Lx, d.Ly, spacing);
256 const { nx, ny } = full;
257 const vars = runDesignScript(d.source + trailer(d, full));
258 // meshgrid matrices are [ny, nx] in column-major order: entry (i, j) sits at i * ny + j.
259 const H = matrix(vars, "steelsim_H", ny, nx, "thickness");
260 const Z = matrix(vars, "steelsim_Z", ny, nx, "profile");
261 let g = full;
262 if (d.hasInside) {
263 const IN = matrix(vars, "steelsim_IN", ny - 1, nx - 1, "inside");
264 g = withCellMask(full, (ci, cj) => IN[ci * (ny - 1) + cj] !== 0);
265 }
266 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
267 let any = false;
268 for (let j = 0; j < ny; j++) {
269 for (let i = 0; i < nx; i++) {
270 if (!insideNode(g, i, j)) continue;
271 any = true;
272 const k = j * nx + i, m = i * ny + j;
273 const x = nodeX(g, i), y = nodeY(g, j);
274 const at = `(${x.toFixed(4)}, ${y.toFixed(4)})`;
275 if (!isFinite(H[m])) throw new DesignError(`thickness is ${H[m]} at ${at}`);
276 if (H[m] <= 0) throw new DesignError(`thickness must be positive; it is ${H[m]} at ${at}`);
277 if (!isFinite(Z[m])) throw new DesignError(`profile is ${Z[m]} at ${at}`);
278 h[k] = H[m];
279 w0[k] = Z[m];
280 }
281 }
282 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
283 removePlane(w0, g);
284 for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) if (!insideNode(g, i, j)) w0[j * nx + i] = 0;
285 return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, mount: d.mount, loss: d.loss };
286}
288/** Peak-to-peak height and the largest slope of the profile over the plate nodes. */
289export function profileStats(pl: Plate): { peakToPeak: number; maxSlope: number } {
290 const g = pl.grid, w0 = pl.w0;
291 let lo = Infinity, hi = -Infinity, maxSlope = 0;
292 const at = (i: number, j: number) => w0[j * g.nx + i];
293 const slope = (i: number, j: number, di: number, dj: number, hg: number) => {
294 const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
295 if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * hg);
296 if (hasP) return (at(i + di, j + dj) - at(i, j)) / hg;
297 if (hasM) return (at(i, j) - at(i - di, j - dj)) / hg;
298 return 0;
299 };
300 for (let j = 0; j < g.ny; j++) {
301 for (let i = 0; i < g.nx; i++) {
302 if (!insideNode(g, i, j)) continue;
303 const v = at(i, j);
304 lo = Math.min(lo, v); hi = Math.max(hi, v);
305 maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
306 }
307 }
308 return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
309}