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