Design scripts in MATLAB syntax, run by numbl - proof-of-concept/steelsim41

proof-of-concept / steelsim41

Design scripts in MATLAB syntax, run by numbl

The design script now fills in a struct named design in MATLAB syntax; numbl (an npm dependency) runs it in the page and in the tests. The outline, thickness, and profile are function handles @(x, y, r, theta) that design.ts samples on meshgrid matrices through an appended trailer, so the fields are evaluated inside numbl, vectorized. The helpers (inch, mm, circle, ellipse, rectangle, clamp, lerp, smoothstep, bump, dome, bow, dents, seededrand, noise2d) are .m workspace files in src/mhelpers.ts. All presets are rewritten in MATLAB style. Errors from the script are reported with their line when they arise in the script itself.

Jeremy Magland <jeremy.magland@gmail.com> committed 1 month ago commit 8db04e05b845 parent 19b58d0 Browse files

10 changed files+2490−466

README.mdmodified+23−18View file

@@ -1,6 +1,6 @@

11 # steelsim41

22

3-A browser application for designing a virtual cymbal and simulating a drumstick strike on it. The cymbal is described by a short JavaScript script (material, bounding box, outline, thickness field, stress-free profile, mount, and losses) that the user edits in the page, or that an AI writes for them. The sound is computed by solving the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly on the GPU through WebGPU. The model and the numerical method are described in [THEORY.md](THEORY.md).

3+A browser application for designing a virtual cymbal and simulating a drumstick strike on it. The cymbal is described by a short script in MATLAB syntax (material, bounding box, outline, thickness field, stress-free profile, mount, and losses), run in the page by [numbl](https://numbl.org), that the user edits in the page or that an AI writes for them. The sound is computed by solving the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly on the GPU through WebGPU. The model and the numerical method are described in [THEORY.md](THEORY.md).

44

55 The physics and most of the interface come from steelsim31 (the hanging steel sheet); the changes here are the scripted design, the thickness field, the mount, an arbitrary outline (including a center hole), and the removal of the push excitation, the hit series, and the feature-list shape editor.

66

@@ -21,23 +21,27 @@ The sidebar has three parts: the design script, the strike and simulation parame

2121

2222 **Design script.** Pick a preset from the dropdown (a flat steel rectangle by default; a bowed rectangle, a flat disc on a stand, a 16 inch crash, the same crash with a dozen dents, a 20 inch ride, an ellipse with a twist, an 18 inch china) and edit it freely; the script is compiled as you type (or with Ctrl+Enter) and the view shows the resulting stress-free shape, colored by height or by thickness. Errors are reported under the script, and a summary lists the bounding box, the material, the thickness range, the mass, the grid, the profile's peak-to-peak height and largest slope (with a warning above the shallow-shell range), the mount, and the losses. Revert restores the selected preset.

2323

24-The script is a plain JavaScript function body that returns the design object. All lengths are meters; `inch` and `mm` are conversion constants. Every field function receives `(x, y, r, theta)`: the point relative to the center of the bounding box, its distance from the center, and its angle from the x axis.

25-

26-```js

27-const D = 16 * inch, R = D / 2;

28-return {

29- name: "16 inch crash",

30- material: { E: 110e9, nu: 0.34, rho: 8800 }, // Pa, -, kg/m^3

31- bounds: { width: D, height: D }, // bounding box [m]

32- inside: circle(D, 12 * mm), // outline, or (x, y, r, theta) => boolean; omit for the whole box

33- thickness: (x, y, r) => lerp(1.4 * mm, 0.8 * mm, clamp((r - 60 * mm) / (R - 60 * mm), 0, 1)),

34- profile: (x, y, r) => 16 * mm * (1 - (r / R) ** 2) + dome(r, 60 * mm, 12 * mm), // stress-free height [m]

35- mount: { radius: 20 * mm, stiffness: 2e4, damping: 5 }, // or null for a free plate

36- loss: { sigma0: 0.6, sigma1: 2.5e-4, sigmaIn: 20 },

37-};

24+The script is written in MATLAB syntax and run by [numbl](https://numbl.org). It fills in a struct named `design`. All lengths are meters; `inch` and `mm` are conversion constants. The outline, the thickness, and the profile are numbers or function handles `@(x, y, r, theta)` of the point relative to the center of the bounding box, its distance from the center, and its angle from the x axis; the handles are applied to whole grid matrices, so they use elementwise operators.

25+

26+```matlab

27+D = 16*inch; R = D/2;

28+bellR = 60*mm; bellH = 12*mm; bowH = 16*mm;

29+design.name = '16 inch crash';

30+design.material.E = 110e9; % Pa

31+design.material.nu = 0.34;

32+design.material.rho = 8800; % kg/m^3

33+design.bounds.width = D; % bounding box [m]

34+design.bounds.height = D;

35+design.inside = circle(D, 12*mm); % outline with a 12 mm hole; omit for the whole box

36+design.thickness = @(x, y, r, theta) lerp(1.4*mm, 0.8*mm, clamp((r - bellR) ./ (R - bellR), 0, 1));

37+design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH); % stress-free height [m]

38+design.mount.radius = 20*mm; % or design.mount = [] for a free plate

39+design.mount.stiffness = 2e4; % N/m

40+design.mount.damping = 5; % N s/m

41+design.loss.sigma0 = 0.6; design.loss.sigma1 = 2.5e-4; design.loss.sigmaIn = 20;

3842 ```

3943

40-`thickness` and `profile` may be numbers or functions; `inside`, `mount`, and `loss` are optional (`mount` may also give `x` and `y`). Helpers in scope: `circle(diameter, hole)`, `ellipse(width, height, hole)`, `rectangle(width, height, hole)` (outlines), `clamp`, `lerp`, `smoothstep(e0, e1, x)`, `bump(d, radius, height)` (Gaussian), `dome(d, radius, height)` (cosine dome with zero slope at the rim), `bow(s, span, sag, offset)` (parabolic bow), and `noise2d(seed, length)` (a smooth random field of unit RMS). The profile's mean plane is removed before use, since only its slopes enter the physics and a tilt would add a spurious membrane coupling.

44+`inside`, `mount`, `profile`, and `loss` are optional (`mount` may also give `x` and `y`). Helper functions in scope: `circle(diameter, hole)`, `ellipse(width, height, hole)`, `rectangle(width, height, hole)` (outlines), `clamp`, `lerp`, `smoothstep(e0, e1, x)`, `bump(d, radius, height)` (Gaussian), `dome(d, radius, height)` (cosine dome with zero slope at the rim), `bow(s, span, sag, offset)` (parabolic bow), `dents(x, y, cx, cy, radius, height)` (a sum of bumps), `seededrand(seed, n)`, and `noise2d(seed, length)` (a smooth random field of unit RMS, returned as a handle `f(x, y)`). Their sources are in `src/mhelpers.ts`. The profile's mean plane is removed before use, since only its slopes enter the physics and a tilt would add a spurious membrane coupling.

4145

4246 **Strike and simulation.** The strike point is given as fractions of the bounding box (the yellow ring in the view), with the stick speed, its effective mass, the Hertz contact stiffness and exponent, and the tip radius. The listener distance and direction can be changed while a simulation runs. The grid spacing (8 mm by default), the in-plane mass factor, the sample rate, and the frame rate take effect at the next Start; the panel reports the grid, the frequencies at which grid dispersion lowers the partials by 3 and 10 percent, the time step, and the 60 dB decay times of the design's losses. The rigid-body projection applies to free plates only; a mounted plate takes the reaction of the strike through the mount.

4347

@@ -52,7 +56,7 @@ The default grid spacing of 8 mm keeps the simulation interactive but lowers the

5256 `npm test` runs three groups of tests with Node's built-in runner.

5357

5458 - `tests/cpuRef.test.ts` checks the float64 reference of the spatial discretization without a GPU: symmetry of the stiffness with exactly three rigid modes per block on uniform, tapered, elliptical, and holed plates, forces as the exact gradient of the discrete energy (flat and bowed, uniform and tapered, free and mounted), dissipativity of the loss terms and the mount's dashpot, the scaling of the bending frequencies with the thickness, the mount's stiffness against rigid translation and tilt, and convergence of the free square and circular plate frequencies to the reference values of Gorman and Leissa.

55-- `tests/design.test.ts` checks the design script: compilation, the checks and their messages, the helpers, the mask of a disc with a hole, the removal of the mean plane, the mount patch, and that every preset evaluates on the default grid.

59+- `tests/design.test.ts` checks the design script: the helpers, compilation and the checks with their messages, the evaluation on the grid (including the orientation of the sampled matrices), the mask of a disc with a hole, the removal of the mean plane, the mount patch, and that every preset evaluates on the default grid.

5660 - `tests/gpu.test.ts` serves the app with Vite and drives a headless Chromium through Playwright. It compares one step of the WGSL kernel with the float64 reference on flat, bowed, tapered, mounted, and elliptical presets (`?selftest=1` in the browser shows the same numbers), runs a short strike on a free plate and on a mounted cymbal to completion, and checks that a script with an error is reported and cannot start. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first two groups.

5761

5862 ## Benchmark

@@ -62,7 +66,8 @@ The default grid spacing of 8 mm keeps the simulation interactive but lowers the

6266 ## Layout

6367

6468 - `src/params.ts`: simulation parameters, grid and mask, the plate as evaluated on the grid, material functions, patches of the strike and the mount, grid dispersion estimate.

65-- `src/design.ts`: the script API, the helpers, compilation and checks, evaluation on the grid.

69+- `src/design.ts`: running the script with numbl, the checks on the design struct, evaluation on the grid.

70+- `src/mhelpers.ts`: the MATLAB helper functions available to a script, as numbl workspace files.

6671 - `src/presets.ts`: the preset scripts.

6772 - `src/sim/cpuRef.ts`: float64 reference of the spatial discretization; used for the time-step estimate and the tests.

6873 - `src/sim/shaders/`: WGSL kernels (time step, rigid-body sums and projection, listener sum, display extraction).

package-lock.jsonmodified+1846−85View file

This diff is 2,049 lines long and is not shown.

package.jsonmodified+3−0View file

@@ -19,5 +19,8 @@

1919 "playwright": "^1.50.0",

2020 "typescript": "^5.6.0",

2121 "vite": "^6.0.0"

22+ },

23+ "dependencies": {

24+ "numbl": "^0.4.21"

2225 }

2326 }

src/design.tsmodified+188−163View file

@@ -1,18 +1,19 @@

1-// The cymbal design: a JavaScript script that returns an object describing

2-// the material, the bounding box, the outline, the thickness field, the

3-// stress-free profile, the mount, and the losses. The script runs in the page

4-// with a set of helper functions in scope (see HELPERS) and returns a plain

5-// object; compileDesign checks it and evaluatePlate samples the fields on a

6-// grid of the requested spacing.

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.

79 //

8-// Lengths are meters throughout; `inch` and `mm` are available as constants.

9-// Every field function receives (x, y, r, theta): the point relative to the

10-// center of the bounding box, its distance from the center, and its angle

11-// from the x axis in radians.

10+// Lengths are meters throughout; `inch` and `mm` are helper functions. The

11+// helpers a script may call are listed in mhelpers.ts.

1212

13+import { executeCode, type RuntimeValue, type RuntimeTensor, type RuntimeStruct } from "numbl";

14+import { HELPER_FILES } from "./mhelpers.ts";

1315 import {

1416 type Grid,

15- type InsideFn,

1617 type Loss,

1718 type Material,

1819 type Mount,

@@ -22,167 +23,166 @@ import {

2223 nodeWeight,

2324 nodeX,

2425 nodeY,

26+ withCellMask,

2527 } from "./params.ts";

2628

27-export type PointFn = (x: number, y: number, r: number, theta: number) => number;

28-export type OutlineFn = (x: number, y: number, r: number, theta: number) => boolean;

29-

30-/** A compiled design: the script's object, checked and with defaults filled in. */

29+/** A compiled design: the checked contents of the `design` struct, plus the script that produced it. */

3130 export interface Design {

3231 name: string;

3332 material: Material;

3433 Lx: number;

3534 Ly: number;

36- inside: OutlineFn | null;

37- thickness: PointFn;

38- profile: PointFn;

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;

3941 mount: Mount | null;

4042 loss: Loss;

43+ source: string;

4144 }

4245

4346 export const DEFAULT_LOSS: Loss = { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 };

4447

45-// ---------- Helpers available to the script ----------

48+/** Message of an error thrown by the script or by the checks. */

49+export class DesignError extends Error {}

4650

47-/** Deterministic uniform generator on [0, 1) (mulberry32). */

48-function mulberry32(seed: number): () => number {

49- let a = Math.floor(seed) >>> 0;

50- return () => {

51- a = (a + 0x6d2b79f5) >>> 0;

52- let t = a;

53- t = Math.imul(t ^ (t >>> 15), t | 1);

54- t ^= t + Math.imul(t ^ (t >>> 7), t | 61);

55- return ((t ^ (t >>> 14)) >>> 0) / 4294967296;

56- };

57-}

51+const FILE_NAME = "design.m";

5852

59-const NOISE_TERMS = 64;

53+/** Runs a design script (with the helpers in scope) and returns its workspace variables. */

54+export 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+}

6061

61-/**

62- * A smooth random field of unit RMS with the covariance exp(-d^2 / length^2):

63- * a sum of plane waves whose wavevectors are Gaussian with standard deviation

64- * sqrt(2) / length per component. The same seed gives the same field on any grid.

65- */

66-export function noise2d(seed: number, length: number): (x: number, y: number) => number {

67- const rnd = mulberry32(seed);

68- const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());

69- const L = Math.max(length, 1e-4);

70- const terms = Array.from({ length: NOISE_TERMS }, () => ({ kx: (gauss() * Math.SQRT2) / L, ky: (gauss() * Math.SQRT2) / L, ph: 2 * Math.PI * rnd() }));

71- const norm = 1 / Math.sqrt(NOISE_TERMS / 2);

72- return (x, y) => {

73- let v = 0;

74- for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);

75- return v * norm;

76- };

62+/** The message of a numbl error (a RuntimeError or a parser SyntaxError), with its line when it lies in the user's script. */

63+function 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}`;

7771 }

7872

79-export const HELPERS = {

80- inch: 0.0254,

81- mm: 1e-3,

82- /** Outline of a circle of the given diameter, with an optional center hole of diameter `hole`. */

83- circle: (diameter: number, hole = 0): OutlineFn => (_x, _y, r) => r <= diameter / 2 && r >= hole / 2,

84- /** Outline of the ellipse with the given width and height, with an optional center hole. */

85- ellipse: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => (x / (width / 2)) ** 2 + (y / (height / 2)) ** 2 <= 1 && r >= hole / 2,

86- /** Outline of the rectangle of the given width and height (centered), with an optional center hole. */

87- rectangle: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => Math.abs(x) <= width / 2 && Math.abs(y) <= height / 2 && r >= hole / 2,

88- clamp: (x: number, lo: number, hi: number): number => Math.min(hi, Math.max(lo, x)),

89- lerp: (a: number, b: number, t: number): number => a + (b - a) * t,

90- /** 0 below e0, 1 above e1, a smooth cubic between. */

91- smoothstep: (e0: number, e1: number, x: number): number => {

92- const t = Math.min(1, Math.max(0, (x - e0) / (e1 - e0)));

93- return t * t * (3 - 2 * t);

94- },

95- /** Gaussian bump of the given height, falling to 1/e at `radius` from the center; `d` is the distance from the center. */

96- bump: (d: number, radius: number, height: number): number => height * Math.exp(-(d * d) / (radius * radius)),

97- /** Cosine dome of the given height and base radius (zero slope at the rim and the apex); `d` is the distance from the center. */

98- dome: (d: number, radius: number, height: number): number => (d >= radius ? 0 : 0.5 * height * (1 + Math.cos((Math.PI * d) / radius))),

99- /** Parabolic bow of the given sag across a span L centered on 0; `s` is the coordinate along the span; the offset in -1..1 moves the crest toward one end. */

100- bow: (s: number, L: number, sag: number, offset = 0): number => {

101- const t = (2 * s) / L;

102- return sag * (1 - t * t) * (1 + offset * t);

103- },

104- noise2d,

105-};

73+// ---------- Reading runtime values ----------

10674

107-export type Helpers = typeof HELPERS;

75+function isStruct(v: RuntimeValue | undefined): v is RuntimeStruct {

76+ return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "struct";

77+}

78+function isTensor(v: RuntimeValue | undefined): v is RuntimeTensor {

79+ return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor";

80+}

81+function isFunction(v: RuntimeValue | undefined): boolean {

82+ return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "function";

83+}

84+function isEmpty(v: RuntimeValue | undefined): boolean {

85+ return v === undefined || (isTensor(v) && v.data.length === 0);

86+}

10887

109-// ---------- Compilation ----------

88+/** A scalar number from a runtime value, or null when it is not one. */

89+function 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+}

11095

111-/** Message of an error thrown by the script or by the checks. */

112-export class DesignError extends Error {}

96+function 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+}

113102

114-function num(v: unknown, what: string, min = -Infinity): number {

115- if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} must be a finite number`);

116- if (v < min) throw new DesignError(`${what} must be at least ${min}`);

117- return v;

103+function field(s: RuntimeStruct, name: string): RuntimeValue | undefined {

104+ return s.fields.get(name);

118105 }

119106

120-function fieldFn(v: unknown, what: string): PointFn {

121- if (typeof v === "number") {

122- if (!isFinite(v)) throw new DesignError(`${what} must be a finite number or a function`);

123- return () => v;

124- }

125- if (typeof v === "function") return v as PointFn;

126- throw new DesignError(`${what} must be a number or a function (x, y, r, theta) => number`);

107+function 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;

127111 }

128112

113+// ---------- Compilation ----------

114+

129115 /**

130- * Runs the script and checks the object it returns. The script body has the

131- * helpers in scope as free variables and must `return` the design object.

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.

132119 */

133120 export function compileDesign(source: string): Design {

134- let raw: unknown;

135- try {

136- const fn = new Function(...Object.keys(HELPERS), source) as (...args: unknown[]) => unknown;

137- raw = fn(...Object.values(HELPERS));

138- } catch (e) {

139- throw new DesignError(`script error: ${e instanceof Error ? e.message : String(e)}`);

140- }

141- if (!raw || typeof raw !== "object") throw new DesignError("the script must return an object (end it with `return { ... };`)");

142- const o = raw as Record<string, unknown>;

143- const name = typeof o.name === "string" && o.name.trim() ? o.name.trim() : "Untitled 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");

144125

145- const m = o.material as Record<string, unknown> | undefined;

146- if (!m || typeof m !== "object") throw new DesignError("material must be an object { E, nu, rho }");

147- const material: Material = { E: num(m.E, "material.E", 1e6), nu: num(m.nu, "material.nu", 0), rho: num(m.rho, "material.rho", 1) };

148- if (material.nu >= 0.5) throw new DesignError("material.nu must be below 0.5");

126+ const nameV = field(d, "name");

127+ const name = (nameV !== undefined ? text(nameV) : null)?.trim() || "Untitled design";

149128

150- const b = o.bounds as Record<string, unknown> | undefined;

151- if (!b || typeof b !== "object") throw new DesignError("bounds must be an object { width, height }");

152- const Lx = num(b.width, "bounds.width", 1e-3), Ly = num(b.height, "bounds.height", 1e-3);

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");

153133

154- let inside: OutlineFn | null = null;

155- if (o.inside !== undefined && o.inside !== null) {

156- if (typeof o.inside !== "function") throw new DesignError("inside must be a function (x, y, r, theta) => boolean");

157- inside = o.inside as OutlineFn;

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);

137+

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+ }

144+

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");

158152 }

159153

160- if (o.thickness === undefined) throw new DesignError("thickness is required (a number in meters, or a function)");

161- const thickness = fieldFn(o.thickness, "thickness");

162- const profile = o.profile === undefined || o.profile === null ? () => 0 : fieldFn(o.profile, "profile");

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+ }

163163

164164 let mount: Mount | null = null;

165- if (o.mount !== undefined && o.mount !== null && o.mount !== false) {

166- const mo = o.mount as Record<string, unknown>;

167- if (typeof mo !== "object") throw new DesignError("mount must be an object { x, y, radius, stiffness, damping } or 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));

168169 mount = {

169- x: mo.x === undefined ? 0 : num(mo.x, "mount.x"),

170- y: mo.y === undefined ? 0 : num(mo.y, "mount.y"),

171- radius: num(mo.radius, "mount.radius", 0),

172- stiffness: num(mo.stiffness, "mount.stiffness", 0),

173- damping: mo.damping === undefined ? 0 : num(mo.damping, "mount.damping", 0),

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),

174175 };

175176 }

176177

177- const lo = (o.loss ?? {}) as Record<string, unknown>;

178- if (typeof lo !== "object") throw new DesignError("loss must be an object { sigma0, sigma1, sigmaIn }");

179- const loss: Loss = {

180- sigma0: lo.sigma0 === undefined ? DEFAULT_LOSS.sigma0 : num(lo.sigma0, "loss.sigma0", 0),

181- sigma1: lo.sigma1 === undefined ? DEFAULT_LOSS.sigma1 : num(lo.sigma1, "loss.sigma1", 0),

182- sigmaIn: lo.sigmaIn === undefined ? DEFAULT_LOSS.sigmaIn : num(lo.sigmaIn, "loss.sigmaIn", 0),

183- };

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+ }

184184

185- return { name, material, Lx, Ly, inside, thickness, profile, mount, loss };

185+ return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, mount, loss, source };

186186 }

187187

188188 // ---------- Evaluation on a grid ----------

@@ -210,53 +210,78 @@ export function removePlane(f: Float64Array, g: Grid): void {

210210 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);

211211 }

212212

213-function callField(fn: PointFn, x: number, y: number, what: string): number {

214- let v: unknown;

215- try {

216- v = fn(x, y, Math.hypot(x, y), Math.atan2(y, x));

217- } catch (e) {

218- throw new DesignError(`${what} threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);

219- }

220- if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} returned ${String(v)} at (${x.toFixed(4)}, ${y.toFixed(4)})`);

221- return v;

213+/** The trailer appended to the script to sample the fields on the grid (all names start with steelsim_). */

214+function 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+}

236+

237+/** A real [rows, cols] matrix from the workspace as a Float64Array in column-major order. */

238+function 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;

222243 }

223244

224245 /**

225246 * Samples the design on a grid of about the requested spacing: the mask from

226- * the outline, the thickness at the nodes (which must be positive on the

227- * plate), and the profile at the nodes with its mean plane removed. Only the

228- * slopes of the profile enter the physics, so its mean is immaterial, but its

229- * tilt is not: the von Karman strains are not invariant under rotation, and a

230- * tilted stress-free shape would add a spurious membrane coupling.

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.

231253 */

232254 export function evaluatePlate(d: Design, spacing: number): Plate {

233- let inside: InsideFn | undefined;

234- if (d.inside) {

235- const fn = d.inside;

236- inside = (x, y) => {

237- try {

238- return !!fn(x, y, Math.hypot(x, y), Math.atan2(y, x));

239- } catch (e) {

240- throw new DesignError(`inside threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);

241- }

242- };

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);

243265 }

244- const g = makeGrid(d.Lx, d.Ly, spacing, inside);

245266 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);

246267 let any = false;

247- for (let j = 0; j < g.ny; j++) {

248- for (let i = 0; i < g.nx; i++) {

268+ for (let j = 0; j < ny; j++) {

269+ for (let i = 0; i < nx; i++) {

249270 if (!insideNode(g, i, j)) continue;

250271 any = true;

251- const x = nodeX(g, i), y = nodeY(g, j), k = j * g.nx + i;

252- h[k] = callField(d.thickness, x, y, "thickness");

253- if (h[k] <= 0) throw new DesignError(`thickness must be positive; it is ${h[k]} at (${x.toFixed(4)}, ${y.toFixed(4)})`);

254- w0[k] = callField(d.profile, x, y, "profile");

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];

255280 }

256281 }

257282 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");

258283 removePlane(w0, g);

259- for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (!insideNode(g, i, j)) w0[j * g.nx + i] = 0;

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;

260285 return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, mount: d.mount, loss: d.loss };

261286 }

262287

src/mhelpers.tsadded+143−0View file

@@ -0,0 +1,143 @@

1+// The MATLAB-syntax helper functions available to a design script, as numbl

2+// workspace files. Each is an ordinary .m function; the design script calls

3+// them like any other function. Everything is vectorized over grid matrices.

4+// steelsim_field is internal: evaluatePlate uses it to sample a field.

5+

6+import type { WorkspaceFile } from "numbl";

7+

8+const files: Record<string, string> = {

9+ "inch.m": `function v = inch()

10+% inch: one inch in meters.

11+v = 0.0254;

12+end

13+`,

14+ "mm.m": `function v = mm()

15+% mm: one millimeter in meters.

16+v = 1e-3;

17+end

18+`,

19+ "circle.m": `function f = circle(diameter, hole)

20+% circle(diameter, hole): outline of a disc with an optional center hole.

21+if nargin < 2, hole = 0; end

22+f = @(x, y, r, theta) (r <= diameter/2) & (r >= hole/2);

23+end

24+`,

25+ "ellipse.m": `function f = ellipse(width, height, hole)

26+% ellipse(width, height, hole): outline of an ellipse with an optional center hole.

27+if nargin < 3, hole = 0; end

28+f = @(x, y, r, theta) ((x/(width/2)).^2 + (y/(height/2)).^2 <= 1) & (r >= hole/2);

29+end

30+`,

31+ "rectangle.m": `function f = rectangle(width, height, hole)

32+% rectangle(width, height, hole): outline of a centered rectangle with an optional center hole.

33+if nargin < 3, hole = 0; end

34+f = @(x, y, r, theta) (abs(x) <= width/2) & (abs(y) <= height/2) & (r >= hole/2);

35+end

36+`,

37+ "clamp.m": `function v = clamp(x, lo, hi)

38+% clamp(x, lo, hi): x limited to [lo, hi], elementwise.

39+v = min(hi, max(lo, x));

40+end

41+`,

42+ "lerp.m": `function v = lerp(a, b, t)

43+% lerp(a, b, t): a at t = 0, b at t = 1, elementwise.

44+v = a + (b - a) .* t;

45+end

46+`,

47+ "smoothstep.m": `function v = smoothstep(e0, e1, x)

48+% smoothstep(e0, e1, x): 0 below e0, 1 above e1, a smooth cubic between.

49+t = clamp((x - e0) ./ (e1 - e0), 0, 1);

50+v = t.^2 .* (3 - 2*t);

51+end

52+`,

53+ "bump.m": `function v = bump(d, radius, height)

54+% bump(d, radius, height): Gaussian bump of the given height, falling to

55+% 1/e at the given radius; d is the distance from its center.

56+v = height * exp(-(d.^2) / radius^2);

57+end

58+`,

59+ "dome.m": `function v = dome(d, radius, height)

60+% dome(d, radius, height): cosine dome of the given height and base radius,

61+% with zero slope at the rim and the apex; d is the distance from its center.

62+v = 0.5 * height * (1 + cos(pi * d / radius)) .* (d < radius);

63+end

64+`,

65+ "bow.m": `function v = bow(s, L, sag, offset)

66+% bow(s, L, sag, offset): parabolic bow of the given sag across a span L

67+% centered on 0; s is the coordinate along the span, and the offset in

68+% -1..1 moves the crest toward one end.

69+if nargin < 4, offset = 0; end

70+t = 2*s/L;

71+v = sag * (1 - t.^2) .* (1 + offset*t);

72+end

73+`,

74+ "dents.m": `function v = dents(x, y, cx, cy, radius, height)

75+% dents(x, y, cx, cy, radius, height): the sum of Gaussian bumps centered

76+% at (cx(k), cy(k)) with the given radii and heights (vectors, or scalars

77+% shared by all the dents).

78+v = zeros(size(x));

79+for k = 1:numel(cx)

80+ rk = radius(min(k, numel(radius)));

81+ hk = height(min(k, numel(height)));

82+ v = v + bump(hypot(x - cx(k), y - cy(k)), rk, hk);

83+end

84+end

85+`,

86+ "seededrand.m": `function u = seededrand(seed, n)

87+% seededrand(seed, n): n uniform numbers in [0, 1) from a small seeded

88+% generator, the same on every run and every machine.

89+s = floor(seed);

90+u = zeros(1, n);

91+for k = 1:n

92+ s = mod(s * 1664525 + 1013904223, 4294967296);

93+ u(k) = s / 4294967296;

94+end

95+end

96+`,

97+ "noise2d.m": `function f = noise2d(seed, len)

98+% noise2d(seed, len): a smooth random field f(x, y) of unit RMS with the

99+% covariance exp(-d^2 / len^2), the same for a given seed on any grid.

100+n = 64;

101+u = reshape(seededrand(seed, 5*n), 5, n);

102+kx = sqrt(-2*log(1 - u(1,:))) .* cos(2*pi*u(2,:)) * sqrt(2) / max(len, 1e-4);

103+ky = sqrt(-2*log(1 - u(3,:))) .* cos(2*pi*u(4,:)) * sqrt(2) / max(len, 1e-4);

104+ph = 2*pi*u(5,:);

105+f = @(x, y) noise2d_eval(kx, ky, ph, x, y);

106+end

107+`,

108+ "noise2d_eval.m": `function v = noise2d_eval(kx, ky, ph, x, y)

109+v = zeros(size(x));

110+for k = 1:numel(kx)

111+ v = v + cos(kx(k) * x + ky(k) * y + ph(k));

112+end

113+v = v / sqrt(numel(kx) / 2);

114+end

115+`,

116+ "steelsim_field.m": `function v = steelsim_field(f, what, x, y, r, theta)

117+% Samples a design field (a number or a function handle) on the grid matrices.

118+if isa(f, 'function_handle')

119+ try

120+ v = f(x, y, r, theta);

121+ catch e

122+ error('%s: %s', what, e.message);

123+ end

124+elseif isnumeric(f) || islogical(f)

125+ v = f;

126+else

127+ error('%s must be a number or a function handle @(x, y, r, theta)', what);

128+end

129+if isscalar(v)

130+ v = double(v) * ones(size(x));

131+elseif ~isequal(size(v), size(x))

132+ error('%s returned a %d x %d array for a %d x %d grid; use elementwise operators (.*, ./, .^)', what, size(v, 1), size(v, 2), size(x, 1), size(x, 2));

133+else

134+ v = double(v);

135+end

136+end

137+`,

138+};

139+

140+export const HELPER_FILES: WorkspaceFile[] = Object.entries(files).map(([name, source]) => ({ name, source }));

141+

142+/** The helper names a script may call (for the interface's hint). */

143+export const HELPER_NAMES = Object.keys(files).filter((n) => !n.startsWith("steelsim_") && !n.startsWith("noise2d_eval")).map((n) => n.replace(/\.m$/, ""));

src/params.tsmodified+6−1View file

@@ -89,10 +89,15 @@ export function makeGrid(Lx: number, Ly: number, spacing: number, inside?: Insid

8989 * center hole) is handled by the same rule.

9090 */

9191 export function withMask(g: Grid, inside: InsideFn): Grid {

92+ return withCellMask(g, (ci, cj) => inside(nodeX(g, ci) + g.hx / 2, nodeY(g, cj) + g.hy / 2));

93+}

94+

95+/** The same restriction with the membership of each cell (ci, cj) given directly. */

96+export function withCellMask(g: Grid, cellInside: (ci: number, cj: number) => boolean): Grid {

9297 const mask = new Uint8Array(g.n);

9398 for (let cj = 0; cj < g.ny - 1; cj++) {

9499 for (let ci = 0; ci < g.nx - 1; ci++) {

95- if (inside(nodeX(g, ci) + g.hx / 2, nodeY(g, cj) + g.hy / 2)) {

100+ if (cellInside(ci, cj)) {

96101 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;

97102 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;

98103 }

src/presets.tsmodified+179−154View file

@@ -1,6 +1,6 @@

1-// Predefined design scripts, selectable from the interface. Each is a plain

2-// JavaScript script body that returns the design object; see design.ts for

3-// the fields and the helpers in scope.

1+// Predefined design scripts, selectable from the interface. Each is a

2+// MATLAB-syntax script that fills in the `design` struct; see design.ts for

3+// the fields and mhelpers.ts for the helper functions in scope.

44

55 export interface Preset {

66 id: string;

@@ -8,171 +8,196 @@ export interface Preset {

88 source: string;

99 }

1010

11-const flatRectangle = `// A flat rectangular sheet of stainless steel, hanging freely (no mount).

12-// Lengths are meters; inch and mm are conversion constants.

13-// Every field can be a number or a function (x, y, r, theta) of the point,

14-// measured from the center of the bounding box.

15-return {

16- name: "Flat steel rectangle",

17- material: { E: 170e9, nu: 0.29, rho: 7900 },

18- bounds: { width: 12 * inch, height: 18 * inch },

19- thickness: 0.79 * mm,

20- profile: 0,

21- mount: null,

22- loss: { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 },

23-};

11+const flatRectangle = `% A flat rectangular sheet of stainless steel, hanging freely (no mount).

12+% Lengths are meters; inch and mm are conversion constants. The script fills

13+% in a struct named design. A field can be a number or a function handle

14+% @(x, y, r, theta) of the point, measured from the center of the bounding

15+% box; the handles are applied to whole grid matrices, so use elementwise

16+% operators (.*, ./, .^) inside them.

17+design.name = 'Flat steel rectangle';

18+design.material.E = 170e9; % Young's modulus [Pa]

19+design.material.nu = 0.29; % Poisson ratio

20+design.material.rho = 7900; % density [kg/m^3]

21+design.bounds.width = 12*inch; % bounding box [m]

22+design.bounds.height = 18*inch;

23+design.thickness = 0.79*mm;

24+design.profile = 0; % stress-free height [m]: flat

25+design.mount = []; % free plate

26+design.loss.sigma0 = 1.15; % frequency-independent loss [1/s]

27+design.loss.sigma1 = 5e-4; % loss proportional to k^2 [m^2/s]

28+design.loss.sigmaIn = 20; % extra loss of the in-plane motion [1/s]

2429 `;

2530

26-const bowedRectangle = `// The same steel sheet with a 3 mm parabolic bow across its width, the

27-// kind of bow that rolling and handling leave in real sheet metal. The bow

28-// makes the sheet a shallow shell, which oil-cans and wobbles when struck.

29-const W = 12 * inch, H = 18 * inch;

30-return {

31- name: "Bowed steel rectangle",

32- material: { E: 170e9, nu: 0.29, rho: 7900 },

33- bounds: { width: W, height: H },

34- thickness: 0.79 * mm,

35- // bow(s, span, sag, offset): parabola of the given sag across the span; the

36- // offset in -1..1 moves the crest toward one edge.

37- profile: (x, y) => bow(x, W, 3 * mm, 0.3),

38- mount: null,

39- loss: { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 },

40-};

31+const bowedRectangle = `% The same steel sheet with a 3 mm parabolic bow across its width, the

32+% kind of bow that rolling and handling leave in real sheet metal. The bow

33+% makes the sheet a shallow shell, which oil-cans and wobbles when struck.

34+W = 12*inch; H = 18*inch;

35+design.name = 'Bowed steel rectangle';

36+design.material.E = 170e9;

37+design.material.nu = 0.29;

38+design.material.rho = 7900;

39+design.bounds.width = W;

40+design.bounds.height = H;

41+design.thickness = 0.79*mm;

42+% bow(s, span, sag, offset): parabola of the given sag across the span; the

43+% offset in -1..1 moves the crest toward one edge.

44+design.profile = @(x, y, r, theta) bow(x, W, 3*mm, 0.3);

45+design.mount = [];

46+design.loss.sigma0 = 1.15;

47+design.loss.sigma1 = 5e-4;

48+design.loss.sigmaIn = 20;

4149 `;

4250

43-const flatDisc = `// A flat 16 inch steel disc with a center hole, sitting on a stand: the

44-// nodes within the mount radius are tied to the ground by a spring (the felt)

45-// and a dashpot, spread over the mount area.

46-const D = 16 * inch;

47-return {

48- name: "Flat steel disc",

49- material: { E: 170e9, nu: 0.29, rho: 7900 },

50- bounds: { width: D, height: D },

51- inside: circle(D, 12 * mm), // diameter, hole diameter

52- thickness: 1.0 * mm,

53- profile: 0,

54- mount: { radius: 20 * mm, stiffness: 2e4, damping: 5 },

55- loss: { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 },

56-};

51+const flatDisc = `% A flat 16 inch steel disc with a center hole, sitting on a stand: the

52+% nodes within the mount radius are tied to the ground by a spring (the felt)

53+% and a dashpot, spread over the mount area.

54+D = 16*inch;

55+design.name = 'Flat steel disc';

56+design.material.E = 170e9;

57+design.material.nu = 0.29;

58+design.material.rho = 7900;

59+design.bounds.width = D;

60+design.bounds.height = D;

61+design.inside = circle(D, 12*mm); % diameter, hole diameter

62+design.thickness = 1.0*mm;

63+design.profile = 0;

64+design.mount.radius = 20*mm; % felt radius

65+design.mount.stiffness = 2e4; % [N/m]

66+design.mount.damping = 5; % [N s/m]

67+design.loss.sigma0 = 1.15;

68+design.loss.sigma1 = 5e-4;

69+design.loss.sigmaIn = 20;

5770 `;

5871

59-const crash16 = `// A 16 inch crash: B20 bronze, tapered from the bell to the edge, with a

60-// shallow bow and a cosine bell. The shallow-shell model behind the

61-// simulation is strained where the slope of the profile exceeds about 0.1

62-// to 0.2, so the bow and the bell are a little flatter than on a real crash.

63-const D = 16 * inch, R = D / 2;

64-const bellR = 60 * mm, bellH = 12 * mm; // base radius and height of the bell

65-const bowH = 16 * mm; // rise from the edge to the foot of the bell

66-return {

67- name: "16 inch crash",

68- material: { E: 110e9, nu: 0.34, rho: 8800 },

69- bounds: { width: D, height: D },

70- inside: circle(D, 12 * mm),

71- // Taper: 1.4 mm at the bell, 0.8 mm at the edge.

72- thickness: (x, y, r) => lerp(1.4 * mm, 0.8 * mm, clamp((r - bellR) / (R - bellR), 0, 1)),

73- // A paraboloidal bow plus the bell; dome(d, radius, height) has zero slope at its rim.

74- profile: (x, y, r) => bowH * (1 - (r / R) ** 2) + dome(r, bellR, bellH),

75- mount: { radius: 20 * mm, stiffness: 2e4, damping: 5 },

76- loss: { sigma0: 0.6, sigma1: 2.5e-4, sigmaIn: 20 },

77-};

72+const crash16 = `% A 16 inch crash: B20 bronze, tapered from the bell to the edge, with a

73+% shallow bow and a cosine bell. The shallow-shell model behind the

74+% simulation is strained where the slope of the profile exceeds about 0.1

75+% to 0.2, so the bow and the bell are a little flatter than on a real crash.

76+D = 16*inch; R = D/2;

77+bellR = 60*mm; bellH = 12*mm; % base radius and height of the bell

78+bowH = 16*mm; % rise from the edge to the foot of the bell

79+design.name = '16 inch crash';

80+design.material.E = 110e9;

81+design.material.nu = 0.34;

82+design.material.rho = 8800;

83+design.bounds.width = D;

84+design.bounds.height = D;

85+design.inside = circle(D, 12*mm);

86+% Taper: 1.4 mm at the bell, 0.8 mm at the edge.

87+design.thickness = @(x, y, r, theta) lerp(1.4*mm, 0.8*mm, clamp((r - bellR) ./ (R - bellR), 0, 1));

88+% A paraboloidal bow plus the bell; dome(d, radius, height) has zero slope at its rim.

89+design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH);

90+design.mount.radius = 20*mm;

91+design.mount.stiffness = 2e4;

92+design.mount.damping = 5;

93+design.loss.sigma0 = 0.6;

94+design.loss.sigma1 = 2.5e-4;

95+design.loss.sigmaIn = 20;

7896 `;

7997

80-const crashDented = `// A 16 inch crash that has been through a lot: the same bronze, taper,

81-// bow, and bell as the clean crash, plus a dozen dents scattered over the

82-// bow and a faint hammering texture. The dents are Gaussian bumps of a few

83-// centimeters (broad compared with the 8 mm grid) and a millimeter or two

84-// in height, up or down; change the seed for a different set. Dents break

85-// the symmetry of the disc, split the degenerate mode pairs, and give the

86-// cymbal a less pure, more trashy tone.

87-const D = 16 * inch, R = D / 2;

88-const bellR = 60 * mm, bellH = 12 * mm, bowH = 16 * mm;

89-// A small seeded generator, so the dents are the same every time.

90-let seed = 7;

91-const rnd = () => { seed = (seed * 1664525 + 1013904223) % 4294967296; return seed / 4294967296; };

92-const dents = [];

93-for (let k = 0; k < 12; k++) {

94- const rr = bellR + 25 * mm + (R - bellR - 45 * mm) * rnd(); // on the bow, clear of the bell and the edge

95- const th = 2 * Math.PI * rnd();

96- dents.push({

97- x: rr * Math.cos(th), y: rr * Math.sin(th),

98- radius: (18 + 14 * rnd()) * mm,

99- height: (rnd() < 0.5 ? -1 : 1) * (0.8 + 1.4 * rnd()) * mm,

100- });

101-}

102-const hammer = noise2d(3, 40 * mm); // unit RMS, correlation length 40 mm

103-return {

104- name: "16 inch crash, dented",

105- material: { E: 110e9, nu: 0.34, rho: 8800 },

106- bounds: { width: D, height: D },

107- inside: circle(D, 12 * mm),

108- thickness: (x, y, r) => lerp(1.4 * mm, 0.8 * mm, clamp((r - bellR) / (R - bellR), 0, 1)),

109- profile: (x, y, r) => {

110- let z = bowH * (1 - (r / R) ** 2) + dome(r, bellR, bellH);

111- for (const d of dents) z += bump(Math.hypot(x - d.x, y - d.y), d.radius, d.height);

112- z += 0.3 * mm * hammer(x, y) * smoothstep(bellR, bellR + 30 * mm, r); // no hammering on the bell

113- return z;

114- },

115- mount: { radius: 20 * mm, stiffness: 2e4, damping: 5 },

116- loss: { sigma0: 0.6, sigma1: 2.5e-4, sigmaIn: 20 },

117-};

98+const crashDented = `% A 16 inch crash that has been through a lot: the same bronze, taper,

99+% bow, and bell as the clean crash, plus a dozen dents scattered over the

100+% bow and a faint hammering texture. The dents are Gaussian bumps of a few

101+% centimeters (broad compared with the 8 mm grid) and a millimeter or two

102+% in height, up or down; change the seed for a different set. Dents break

103+% the symmetry of the disc, split the degenerate mode pairs, and give the

104+% cymbal a less pure, more trashy tone.

105+D = 16*inch; R = D/2;

106+bellR = 60*mm; bellH = 12*mm; bowH = 16*mm;

107+% Twelve dents from a seeded generator (the same on every run).

108+nd = 12;

109+u = reshape(seededrand(7, 4*nd), 4, nd); % rows: radius, angle, size, height

110+rr = bellR + 25*mm + (R - bellR - 45*mm) * u(1,:); % on the bow, clear of the bell and the edge

111+th = 2*pi * u(2,:);

112+cx = rr .* cos(th); cy = rr .* sin(th);

113+rad = (18 + 14*u(3,:)) * mm;

114+hgt = sign(u(4,:) - 0.5) .* (0.8 + 1.4*u(4,:)) * mm;

115+hammer = noise2d(3, 40*mm); % unit RMS, correlation length 40 mm

116+design.name = '16 inch crash, dented';

117+design.material.E = 110e9;

118+design.material.nu = 0.34;

119+design.material.rho = 8800;

120+design.bounds.width = D;

121+design.bounds.height = D;

122+design.inside = circle(D, 12*mm);

123+design.thickness = @(x, y, r, theta) lerp(1.4*mm, 0.8*mm, clamp((r - bellR) ./ (R - bellR), 0, 1));

124+design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH) ...

125+ + dents(x, y, cx, cy, rad, hgt) ...

126+ + 0.3*mm * hammer(x, y) .* smoothstep(bellR, bellR + 30*mm, r); % no hammering on the bell

127+design.mount.radius = 20*mm;

128+design.mount.stiffness = 2e4;

129+design.mount.damping = 5;

130+design.loss.sigma0 = 0.6;

131+design.loss.sigma1 = 2.5e-4;

132+design.loss.sigmaIn = 20;

118133 `;

119134

120-const ride20 = `// A 20 inch ride: heavier bronze, a larger bell, a flatter bow. The taper

121-// is stronger than on the crash, and the mount is a little stiffer.

122-const D = 20 * inch, R = D / 2;

123-const bellR = 70 * mm, bellH = 14 * mm;

124-const bowH = 15 * mm;

125-return {

126- name: "20 inch ride",

127- material: { E: 110e9, nu: 0.34, rho: 8800 },

128- bounds: { width: D, height: D },

129- inside: circle(D, 12 * mm),

130- thickness: (x, y, r) => lerp(2.0 * mm, 1.1 * mm, smoothstep(bellR, R, r)),

131- profile: (x, y, r) => bowH * (1 - (r / R) ** 2) + dome(r, bellR, bellH),

132- mount: { radius: 22 * mm, stiffness: 3e4, damping: 8 },

133- loss: { sigma0: 0.5, sigma1: 2e-4, sigmaIn: 20 },

134-};

135+const ride20 = `% A 20 inch ride: heavier bronze, a larger bell, a flatter bow. The taper

136+% is stronger than on the crash, and the mount is a little stiffer.

137+D = 20*inch; R = D/2;

138+bellR = 70*mm; bellH = 14*mm;

139+bowH = 15*mm;

140+design.name = '20 inch ride';

141+design.material.E = 110e9;

142+design.material.nu = 0.34;

143+design.material.rho = 8800;

144+design.bounds.width = D;

145+design.bounds.height = D;

146+design.inside = circle(D, 12*mm);

147+design.thickness = @(x, y, r, theta) lerp(2.0*mm, 1.1*mm, smoothstep(bellR, R, r));

148+design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH);

149+design.mount.radius = 22*mm;

150+design.mount.stiffness = 3e4;

151+design.mount.damping = 8;

152+design.loss.sigma0 = 0.5;

153+design.loss.sigma1 = 2e-4;

154+design.loss.sigmaIn = 20;

135155 `;

136156

137-const ellipseTwist = `// A steel ellipse with a saddle (twist) and a bow, hanging freely. The

138-// twist raises two opposite quadrants and lowers the other two, which

139-// breaks the symmetry of the disc modes.

140-const W = 14 * inch, H = 20 * inch;

141-return {

142- name: "Steel ellipse with a twist",

143- material: { E: 170e9, nu: 0.29, rho: 7900 },

144- bounds: { width: W, height: H },

145- inside: ellipse(W, H),

146- thickness: 0.79 * mm,

147- profile: (x, y) => {

148- const X = (2 * x) / W, Y = (2 * y) / H; // -1..1 over the bounding box

149- return 2 * mm * X * Y + bow(y, H, 2 * mm);

150- },

151- mount: null,

152- loss: { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 },

153-};

157+const ellipseTwist = `% A steel ellipse with a saddle (twist) and a bow, hanging freely. The

158+% twist raises two opposite quadrants and lowers the other two, which

159+% breaks the symmetry of the disc modes.

160+W = 14*inch; H = 20*inch;

161+design.name = 'Steel ellipse with a twist';

162+design.material.E = 170e9;

163+design.material.nu = 0.29;

164+design.material.rho = 7900;

165+design.bounds.width = W;

166+design.bounds.height = H;

167+design.inside = ellipse(W, H);

168+design.thickness = 0.79*mm;

169+% 2x/W and 2y/H run from -1 to 1 over the bounding box.

170+design.profile = @(x, y, r, theta) 2*mm * (2*x/W) .* (2*y/H) + bow(y, H, 2*mm);

171+design.mount = [];

172+design.loss.sigma0 = 1.15;

173+design.loss.sigma1 = 5e-4;

174+design.loss.sigmaIn = 20;

154175 `;

155176

156-const china18 = `// An 18 inch china: a conical bow that falls from the bell toward a flange,

157-// and an upturned edge that rises again over the outer 45 mm. The bell is

158-// square-topped (a plateau with rounded shoulders). The edge slope is well

159-// beyond the shallow-shell range, so this preset is the least trustworthy.

160-const D = 18 * inch, R = D / 2;

161-const bellR = 65 * mm, bellH = 14 * mm;

162-const flangeR = R - 45 * mm, bowH = 14 * mm, edgeH = 12 * mm;

163-return {

164- name: "18 inch china",

165- material: { E: 110e9, nu: 0.34, rho: 8800 },

166- bounds: { width: D, height: D },

167- inside: circle(D, 12 * mm),

168- thickness: (x, y, r) => lerp(1.6 * mm, 0.9 * mm, smoothstep(bellR, flangeR, r)),

169- profile: (x, y, r) =>

170- bowH * Math.max(0, 1 - r / flangeR) // cone down to the flange

171- + edgeH * smoothstep(flangeR, R, r) // upturned edge

172- + bellH * (1 - smoothstep(bellR - 20 * mm, bellR, r)), // flat-topped bell

173- mount: { radius: 20 * mm, stiffness: 2e4, damping: 5 },

174- loss: { sigma0: 0.7, sigma1: 3e-4, sigmaIn: 20 },

175-};

177+const china18 = `% An 18 inch china: a conical bow that falls from the bell toward a flange,

178+% and an upturned edge that rises again over the outer 45 mm. The bell is

179+% square-topped (a plateau with rounded shoulders). The edge slope is well

180+% beyond the shallow-shell range, so this preset is the least trustworthy.

181+D = 18*inch; R = D/2;

182+bellR = 65*mm; bellH = 14*mm;

183+flangeR = R - 45*mm; bowH = 14*mm; edgeH = 12*mm;

184+design.name = '18 inch china';

185+design.material.E = 110e9;

186+design.material.nu = 0.34;

187+design.material.rho = 8800;

188+design.bounds.width = D;

189+design.bounds.height = D;

190+design.inside = circle(D, 12*mm);

191+design.thickness = @(x, y, r, theta) lerp(1.6*mm, 0.9*mm, smoothstep(bellR, flangeR, r));

192+design.profile = @(x, y, r, theta) bowH * max(0, 1 - r/flangeR) ... % cone down to the flange

193+ + edgeH * smoothstep(flangeR, R, r) ... % upturned edge

194+ + bellH * (1 - smoothstep(bellR - 20*mm, bellR, r)); % flat-topped bell

195+design.mount.radius = 20*mm;

196+design.mount.stiffness = 2e4;

197+design.mount.damping = 5;

198+design.loss.sigma0 = 0.7;

199+design.loss.sigma1 = 3e-4;

200+design.loss.sigmaIn = 20;

176201 `;

177202

178203 export const PRESETS: Preset[] = [

src/ui/editor.tsmodified+2−1View file

@@ -4,6 +4,7 @@

44

55 import { type Design, DesignError, compileDesign } from "../design.ts";

66 import { DEFAULT_PRESET, PRESETS } from "../presets.ts";

7+import { HELPER_NAMES } from "../mhelpers.ts";

78

89 const CUSTOM = "custom";

910 const DEBOUNCE_MS = 400;

@@ -79,7 +80,7 @@ export class DesignEditor {

7980 this.summary.className = "editor-summary";

8081 const hint = document.createElement("p");

8182 hint.className = "sub";

82- hint.textContent = "The script returns the design: material, bounds, outline, thickness, profile, mount, and losses, in SI units. It is applied as you type (or with Ctrl+Enter). Helpers: inch, mm, circle, ellipse, rectangle, clamp, lerp, smoothstep, bump, dome, bow, noise2d; every field function receives (x, y, r, theta).";

83+ hint.textContent = `MATLAB syntax (run by numbl). The script fills in the struct design: material, bounds, inside, thickness, profile, mount, and loss, in SI units; it is applied as you type (or with Ctrl+Enter). The thickness, profile, and inside fields are numbers or handles @(x, y, r, theta) applied to whole grid matrices, so use .*, ./, .^ in them. Helpers: ${HELPER_NAMES.join(", ")}.`;

8384 root.append(head, this.textarea, this.status, this.summary, hint);

8485 }

8586

tests/design.test.tsmodified+99−43View file

@@ -1,65 +1,104 @@

1-// Tests of the design script: compilation, the checks, the helpers, the

2-// evaluation on the grid, and the presets.

1+// Tests of the design script: the MATLAB-syntax helpers, compilation and

2+// the checks, the evaluation on the grid, and the presets.

33

44 import { test } from "node:test";

55 import assert from "node:assert/strict";

6-import { DesignError, HELPERS, compileDesign, evaluatePlate, noise2d, profileStats, removePlane } from "../src/design.ts";

6+import { DesignError, compileDesign, evaluatePlate, profileStats, removePlane, runDesignScript } from "../src/design.ts";

77 import { PRESETS } from "../src/presets.ts";

88 import { countInside, insideNode, mountPatch, nodeWeight, nodeX, nodeY, thicknessRange } from "../src/params.ts";

99

10-const minimal = `return { material: { E: 170e9, nu: 0.29, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3 };`;

10+const minimal = `design.material.E = 170e9; design.material.nu = 0.29; design.material.rho = 7900;

11+design.bounds.width = 0.3; design.bounds.height = 0.2;

12+design.thickness = 1e-3;`;

13+

14+/** A scalar from the workspace. */

15+function num(vars: Record<string, unknown>, name: string): number {

16+ const v = vars[name] as number | { data: Float64Array };

17+ return typeof v === "number" ? v : v.data[0];

18+}

19+function data(vars: Record<string, unknown>, name: string): Float64Array {

20+ return (vars[name] as { data: Float64Array }).data;

21+}

1122

1223 test("a minimal script compiles with defaults filled in", () => {

1324 const d = compileDesign(minimal);

1425 assert.equal(d.name, "Untitled design");

15- assert.equal(d.inside, null);

26+ assert.equal(d.hasInside, false);

1627 assert.equal(d.mount, null);

17- assert.equal(d.thickness(0, 0, 0, 0), 1e-3);

18- assert.equal(d.profile(0.1, 0.05, 0, 0), 0);

28+ assert.equal(d.thicknessConst, 1e-3);

29+ assert.equal(d.profileConst, 0);

1930 assert.deepEqual(d.loss, { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 });

31+ assert.ok(Math.abs(d.Lx - 0.3) < 1e-12 && Math.abs(d.Ly - 0.2) < 1e-12);

2032 });

2133

22-test("the helpers are in scope and lengths are SI", () => {

23- const d = compileDesign(`return { name: "x", material: { E: 1e11, nu: 0.3, rho: 8000 }, bounds: { width: 10 * inch, height: 100 * mm }, thickness: (x, y, r) => lerp(2 * mm, 1 * mm, clamp(r / 0.1, 0, 1)), profile: (x, y, r) => dome(r, 0.05, 0.01) + bump(r, 0.02, 0.001) + bow(x, 0.254, 0.002) };`);

24- assert.ok(Math.abs(d.Lx - 0.254) < 1e-12 && Math.abs(d.Ly - 0.1) < 1e-12);

25- assert.ok(Math.abs(d.thickness(0, 0, 0, 0) - 2e-3) < 1e-12);

26- assert.ok(Math.abs(d.thickness(0, 0, 0.1, 0) - 1e-3) < 1e-12);

27- assert.ok(Math.abs(d.profile(0, 0, 0, 0) - (0.01 + 0.001 + 0.002)) < 1e-12);

28- assert.ok(Math.abs(d.profile(0, 0, 0.06, 0) - (0.002 + 0.001 * Math.exp(-9))) < 1e-15);

29- assert.equal(HELPERS.smoothstep(0, 1, -1), 0);

30- assert.equal(HELPERS.smoothstep(0, 1, 2), 1);

31- assert.equal(HELPERS.smoothstep(0, 1, 0.5), 0.5);

34+test("the helpers are in scope, vectorized, and in SI", () => {

35+ const v = runDesignScript(`

36+a = 10*inch; b = 100*mm;

37+t1 = lerp(2*mm, 1*mm, clamp(0.05/0.1, 0, 1));

38+t2 = lerp(2*mm, 1*mm, clamp(0.2/0.1, 0, 1));

39+p0 = dome(0, 0.05, 0.01) + bump(0, 0.02, 0.001) + bow(0, 0.254, 0.002);

40+p1 = dome(0.06, 0.05, 0.01) + bump(0.06, 0.02, 0.001);

41+s = [smoothstep(0, 1, -1), smoothstep(0, 1, 2), smoothstep(0, 1, 0.5)];

42+f = circle(0.2, 0.03);

43+[X, Y] = meshgrid(-0.1:0.05:0.1, -0.1:0.05:0.1);

44+inside = f(X, Y, hypot(X, Y), atan2(Y, X));

45+dd = dents(X, Y, [0 0.05], [0 0], 0.02, [1e-3 -1e-3]);

46+u = seededrand(5, 3); u2 = seededrand(5, 3);

47+`);

48+ assert.ok(Math.abs(num(v, "a") - 0.254) < 1e-12 && Math.abs(num(v, "b") - 0.1) < 1e-12);

49+ assert.ok(Math.abs(num(v, "t1") - 1.5e-3) < 1e-12);

50+ assert.ok(Math.abs(num(v, "t2") - 1e-3) < 1e-12);

51+ assert.ok(Math.abs(num(v, "p0") - (0.01 + 0.001 + 0.002)) < 1e-12);

52+ assert.ok(Math.abs(num(v, "p1") - 0.001 * Math.exp(-9)) < 1e-15);

53+ assert.deepEqual(Array.from(data(v, "s")), [0, 1, 0.5]);

54+ const inside = data(v, "inside"); // 5 x 5, column-major; center at index 12, corner at 0, hole at the center

55+ assert.equal(inside[12], 0, "the center lies in the hole");

56+ assert.equal(inside[0], 0, "the corner lies outside");

57+ assert.equal(inside[2], 1, "the point (-0.1, 0) lies on the rim");

58+ const dd = data(v, "dd");

59+ assert.ok(Math.abs(dd[12] - (1e-3 - 1e-3 * Math.exp(-(0.05 ** 2) / 0.02 ** 2))) < 1e-15, "two dents of opposite sign");

60+ assert.deepEqual(Array.from(data(v, "u")), Array.from(data(v, "u2")), "seededrand is deterministic");

61+ assert.ok(data(v, "u").every((x) => x >= 0 && x < 1));

3262 });

3363

3464 test("scripts that do not describe a plate are rejected with a message", () => {

3565 const bad: [string, RegExp][] = [

36- [`return 5;`, /must return an object/],

37- [`this is not javascript`, /script error/],

38- [`throw new Error("boom");`, /boom/],

39- [`return { bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3 };`, /material/],

40- [`return { material: { E: 170e9, nu: 0.6, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3 };`, /nu/],

41- [`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3 }, thickness: 1e-3 };`, /bounds.height/],

42- [`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 } };`, /thickness/],

43- [`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: "thin" };`, /thickness/],

44- [`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, inside: 5 };`, /inside/],

45- [`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, mount: { stiffness: 1e4 } };`, /mount.radius/],

66+ [`x = 5;`, /struct named design/],

67+ [`design = 5;`, /design must be a struct/],

68+ [`x = (1 + ;`, /line 1/],

69+ [`error('boom');`, /boom/],

70+ [`design.bounds.width = 0.3; design.bounds.height = 0.2; design.thickness = 1e-3;`, /material/],

71+ [`design.material.E = 170e9; design.material.nu = 0.6; design.material.rho = 7900; design.bounds.width = 0.3; design.bounds.height = 0.2; design.thickness = 1e-3;`, /nu/],

72+ [`design.material.E = 170e9; design.material.nu = 0.3; design.material.rho = 7900; design.bounds.width = 0.3; design.thickness = 1e-3;`, /bounds.height/],

73+ [`design.material.E = 170e9; design.material.nu = 0.3; design.material.rho = 7900; design.bounds.width = 0.3; design.bounds.height = 0.2;`, /thickness/],

74+ [`${minimal} design.thickness = 'thin';`, /thickness must be a number or a function handle/],

75+ [`${minimal} design.thickness = -1;`, /thickness must be positive/],

76+ [`${minimal} design.inside = 5;`, /inside must be a function handle/],

77+ [`${minimal} design.mount.stiffness = 1e4;`, /mount.radius/],

78+ [`${minimal} design.loss = 3;`, /loss must be a struct/],

4679 ];

4780 for (const [src, re] of bad) assert.throws(() => compileDesign(src), (e: unknown) => e instanceof DesignError && re.test(e.message), src);

4881 });

4982

50-test("evaluation rejects a thickness that is not positive and a profile that is not finite", () => {

51- const zero = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: (x) => x };`);

83+test("evaluation rejects fields that are not positive, not finite, not vectorized, or that throw", () => {

84+ const zero = compileDesign(`${minimal} design.thickness = @(x, y, r, theta) x;`);

5285 assert.throws(() => evaluatePlate(zero, 0.02), /thickness must be positive/);

53- const nan = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, profile: (x, y) => Math.sqrt(x) };`);

54- assert.throws(() => evaluatePlate(nan, 0.02), /profile returned NaN/);

55- const thrower = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, inside: () => { throw new Error("no"); } };`);

56- assert.throws(() => evaluatePlate(thrower, 0.02), /inside threw/);

57- const empty = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, inside: () => false };`);

86+ const nan = compileDesign(`${minimal} design.profile = @(x, y, r, theta) 0*x./0;`);

87+ assert.throws(() => evaluatePlate(nan, 0.02), /profile is NaN/);

88+ const shape = compileDesign(`${minimal} design.profile = @(x, y, r, theta) [1 2 3];`);

89+ assert.throws(() => evaluatePlate(shape, 0.02), /profile returned a 1 x 3 array/);

90+ const thrower = compileDesign(`${minimal} design.inside = @(x, y, r, theta) no_such_function(x);`);

91+ assert.throws(() => evaluatePlate(thrower, 0.02), /inside: .*no_such_function/);

92+ const empty = compileDesign(`${minimal} design.inside = @(x, y, r, theta) false & x;`);

5893 assert.throws(() => evaluatePlate(empty, 0.02), /no grid cells/);

94+ const scalarFn = evaluatePlate(compileDesign(`${minimal} design.thickness = @(x, y, r, theta) 2e-3;`), 0.02);

95+ assert.equal(thicknessRange(scalarFn).min, 2e-3, "a handle returning a scalar is spread over the grid");

5996 });

6097

6198 test("a circle with a hole excludes the center and the corners", () => {

62- const d = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.2, height: 0.2 }, thickness: 1e-3, inside: circle(0.2, 0.03) };`);

99+ const d = compileDesign(`design.material.E = 170e9; design.material.nu = 0.3; design.material.rho = 7900;

100+design.bounds.width = 0.2; design.bounds.height = 0.2; design.thickness = 1e-3;

101+design.inside = circle(0.2, 0.03);`);

63102 const pl = evaluatePlate(d, 0.01);

64103 const g = pl.grid;

65104 assert.equal(g.nx, 21);

@@ -75,8 +114,21 @@ test("a circle with a hole excludes the center and the corners", () => {

75114 assert.equal(thicknessRange(pl).max, 1e-3);

76115 });

77116

117+test("fields are sampled at the right nodes (meshgrid orientation)", () => {

118+ const d = compileDesign(`${minimal} design.thickness = @(x, y, r, theta) 10 + x + 10*y;`);

119+ const pl = evaluatePlate(d, 0.05);

120+ const g = pl.grid;

121+ assert.ok(g.nx !== g.ny, "the grid is not square, so a transposed sampling would show");

122+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {

123+ const want = 10 + nodeX(g, i) + 10 * nodeY(g, j);

124+ assert.ok(Math.abs(pl.h[j * g.nx + i] - want) < 1e-12, `node ${i}, ${j}: ${pl.h[j * g.nx + i]} vs ${want}`);

125+ }

126+});

127+

78128 test("the profile is made plane-free and the mount patch is found", () => {

79- const d = compileDesign(`return { material: { E: 170e9, nu: 0.3, rho: 7900 }, bounds: { width: 0.3, height: 0.2 }, thickness: 1e-3, profile: (x, y) => 0.01 + 0.05 * x - 0.02 * y + 0.001 * x * x, mount: { radius: 0.02, stiffness: 1e4, damping: 1 } };`);

129+ const d = compileDesign(`${minimal}

130+design.profile = @(x, y, r, theta) 0.01 + 0.05*x - 0.02*y + 0.001*x.^2;

131+design.mount.radius = 0.02; design.mount.stiffness = 1e4; design.mount.damping = 1;`);

80132 const pl = evaluatePlate(d, 0.01);

81133 const g = pl.grid;

82134 let s = 0, sx = 0, sy = 0;

@@ -101,17 +153,21 @@ test("removePlane leaves a plane-free field unchanged", () => {

101153 const before = Float64Array.from(f);

102154 removePlane(f, g);

103155 // The parabola x^2 has a nonzero mean but no tilt; only the mean should move.

104- let shift = f[0] - before[0];

156+ const shift = f[0] - before[0];

105157 for (let k = 0; k < g.n; k++) assert.ok(Math.abs(f[k] - before[k] - shift) < 1e-12);

106158 });

107159

108160 test("noise2d is deterministic with unit RMS", () => {

109- const a = noise2d(7, 0.05), b = noise2d(7, 0.05), c = noise2d(8, 0.05);

110- assert.equal(a(0.01, 0.02), b(0.01, 0.02));

111- assert.notEqual(a(0.01, 0.02), c(0.01, 0.02));

112- let s2 = 0, n = 0;

113- for (let x = -0.5; x < 0.5; x += 0.01) for (let y = -0.5; y < 0.5; y += 0.01) { s2 += a(x, y) ** 2; n++; }

114- const rms = Math.sqrt(s2 / n);

161+ const v = runDesignScript(`

162+a = noise2d(7, 0.05); b = noise2d(7, 0.05); c = noise2d(8, 0.05);

163+[X, Y] = meshgrid(-0.5:0.01:0.49, -0.5:0.01:0.49);

164+A = a(X, Y); B = b(X, Y); C = c(X, Y);

165+same = isequal(A, B); different = ~isequal(A, C);

166+rms = sqrt(mean(A(:).^2));

167+`);

168+ assert.equal(v.same, true);

169+ assert.equal(v.different, true);

170+ const rms = num(v, "rms");

115171 assert.ok(Math.abs(rms - 1) < 0.2, `rms ${rms}`);

116172 });

117173

tests/gpu.test.tsmodified+1−1View file

@@ -124,7 +124,7 @@ test("a script with an error is reported and does not start", { timeout: 120_000

124124 const ready = await page.evaluate<boolean>("!!window.steelsim");

125125 if (ready) break;

126126 }

127- await page.evaluate("window.steelsim.setScript('return { material: { E: 1e11, nu: 0.3, rho: 8000 }, bounds: { width: 0.3, height: 0.2 }, thickness: -1 };')");

127+ await page.evaluate("window.steelsim.setScript('design.material.E = 1e11; design.material.nu = 0.3; design.material.rho = 8000; design.bounds.width = 0.3; design.bounds.height = 0.2; design.thickness = -1;')");

128128 await page.waitForTimeout(200);

129129 const d = await page.evaluate<Status["design"]>("window.steelsim.design()");

130130 assert.equal(d.name, null);