Add elliptical sheet shape on a masked grid
A shape selector makes the sheet the ellipse inscribed in the width and
height. Cells belong to the sheet when their centre lies inside the
ellipse; the variational finite-difference bending, the Q4 membrane,
the lumped mass, the rigid-body projection and the radiation sum all
follow from a per-node weight (fraction of surrounding cells in the
sheet), which reduces to the trapezoidal weights on a rectangle. The
hole and pushed patch move to the rim. Tests add masked-grid symmetry,
rigid-mode and energy-gradient checks and a free circular plate against
Leissa's values (about 2 percent at 24 cells across the diameter).
13 changed files+424−100
README.mdmodified+4−3View file
@@ -1,6 +1,6 @@
11 # steelsim21
22
3-A browser simulation of the sound of a thin rectangular stainless steel sheet that hangs from a thread, is pushed by hand at the far corner, held until it settles, and released. The sheet rings like a gong, so the simulator solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia retained, stepped explicitly on the GPU through WebGPU. The model and the numerical method are described in [THEORY.md](THEORY.md).
3+A browser simulation of the sound of a thin stainless steel sheet, rectangular or elliptical, that hangs from a thread, is pushed by hand at the far corner, held until it settles, and released. The sheet rings like a gong, so the simulator solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia retained, stepped explicitly on the GPU through WebGPU. The model and the numerical method are described in [THEORY.md](THEORY.md).
44
55 The app and the gallery are deployed by the `deploy` workflow on every push to `main`: https://steelsim21--proof-of-concept.vault1-sites.magland.org/ (gallery at `/docs/gallery/index.html`).
66
@@ -17,13 +17,14 @@ Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefo
1717
1818 Press Start. With the default excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles, and then released by ramping the finger force to zero over the release time; the recording of audio and movie frames starts at the release. The other excitation is a drumstick strike (Hertz contact with a free mass) at a configurable point and speed; the recording then starts at the impact. The view shows the transverse displacement as a color map (2D) or a lit height field (3D, drag to orbit and wheel to zoom), the spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.
1919
20-- The sheet has a stress-free initial bow (default 3 mm across the width; set the shape to Flat for the ideal plate). A negative push amount pushes against the bow.
20+- The sheet is a rectangle by default; the Shape selector makes it the ellipse inscribed in the given width and height. The ellipse is represented on the same grid with a stepped edge (see THEORY.md, Section 4), the hole sits near the upper-left part of the rim and the push acts on the rim point opposite it.
21+- The sheet has a stress-free initial bow (default 3 mm across the width; set the bow shape to Flat for the ideal plate). A negative push amount pushes against the bow.
2122 - Plate, initial shape, excitation, and numerics parameters take effect at the next Start. Damping, listener position, and the radiation model can be changed while a simulation runs.
2223 - Pause halts the run temporarily; Stop ends it, keeps the audio and movie, and unlocks the parameter fields for the next Start.
2324 - Playback level: "Absolute" maps 1 Pa to full scale, shifted by the gain slider, so different runs can be compared by ear (the peak level in dBFS is shown, with a clipping warning); "Loudness" and "Peak" normalize each run instead. The WAV download uses the same setting.
2425 - After the simulation completes, drag the time slider to explore the movie; Play starts audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.
2526 - To extend a finished simulation, raise the duration and press Continue; the state is kept on the GPU.
26-- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1`.
27+- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?shape=ellipse`.
2728
2829 The default grid spacing of 4 mm (about 8.9k nodes) resolves bending waves to about 20 kHz. With in-plane inertia the stable time step is set by the in-plane wave speed, about 0.65 microseconds, so four seconds of sound take about six million steps. This runs at roughly 35k steps per second on an Intel Iris Xe integrated GPU (about three minutes for four seconds) and considerably faster on a discrete GPU. A spacing of 6 mm is about three times faster and adequate below 10 kHz.
2930
THEORY.mdmodified+9−7View file
@@ -1,10 +1,10 @@
11 # Theory and numerical method
22
3-This note describes the model behind the simulator: a thin stainless steel sheet, suspended by a thread, pushed at a corner, held, and released. The sheet rings like a gong, which means that its vibration amplitude (about ten thicknesses here) is far beyond the range where a linear plate model is accurate. We therefore solve the geometrically nonlinear plate equations of Foppl and von Karman with the in-plane inertia retained, discretize them on a regular grid, and step them explicitly on the GPU.
3+This note describes the model behind the simulator: a thin stainless steel sheet, rectangular or elliptical, suspended by a thread, pushed at a point of its rim, held, and released. The sheet rings like a gong, which means that its vibration amplitude (about ten thicknesses here) is far beyond the range where a linear plate model is accurate. We therefore solve the geometrically nonlinear plate equations of Foppl and von Karman with the in-plane inertia retained, discretize them on a regular grid, and step them explicitly on the GPU.
44
55 ## 1. Plate model
66
7-The sheet occupies the rectangle $[-L_x/2, L_x/2] \times [-L_y/2, L_y/2]$ and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. Its motion is described by three displacement fields, the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$. With the membrane strains
7+The sheet occupies the rectangle $[-L_x/2, L_x/2] \times [-L_y/2, L_y/2]$, or the ellipse inscribed in it, and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. Its motion is described by three displacement fields, the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$. With the membrane strains
88
99 $$
1010 \varepsilon_{xx} = u_x + \tfrac12 w_x^2, \qquad
@@ -41,11 +41,11 @@ the membrane forces act through the total slope $(w_0 + w)$, and bending acts on
4141
4242 ## 2. Boundary conditions and suspension
4343
44-All four edges are free: the bending moment and the Kirchhoff shear vanish, and so does the in-plane traction. Rather than imposing these conditions explicitly, the discretization of Section 4 is derived from a discrete energy, and free edges are then the natural boundary conditions of that energy. This avoids the intricate ghost-point formulas that free edges otherwise require, and it guarantees a symmetric stiffness operator.
44+The whole edge is free: the bending moment and the Kirchhoff shear vanish, and so does the in-plane traction. Rather than imposing these conditions explicitly, the discretization of Section 4 is derived from a discrete energy, and free edges are then the natural boundary conditions of that energy. This avoids the intricate ghost-point formulas that free edges otherwise require, and it guarantees a symmetric stiffness operator.
4545
46-The thread through the hole near one corner is not modelled as a constraint. At audio frequencies the thread is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the plate's average position. We therefore treat the plate as free and remove the six rigid-body modes (three translations, three rotations) from the displacement and velocity fields by a mass-weighted projection, applied every few time steps. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$, which is negligible.
46+The thread through the hole near one corner (for an ellipse, near the rim in the upper-left direction, at a parametric angle of 135 degrees) is not modelled as a constraint. At audio frequencies the thread is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the plate's average position. We therefore treat the plate as free and remove the six rigid-body modes (three translations, three rotations) from the displacement and velocity fields by a mass-weighted projection, applied every few time steps. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$, which is negligible.
4747
48-The same projection defines what "holding" the sheet means during the push. The finger prescribes the displacement of a small patch at the far corner. The reaction to the finger force is distributed over the plate as an inertial load (this is the inertia-relief equilibrium of a free body), so the push produces pure bending of the prescribed amount. This is a modelling choice: a sheet hanging freely from a thread and pushed 10 mm at a corner would mostly swing as a pendulum and bend by only about 1.5 mm.
48+The same projection defines what "holding" the sheet means during the push. The finger prescribes the displacement of a small patch at the far corner (for an ellipse, at the rim point opposite the hole). The reaction to the finger force is distributed over the plate as an inertial load (this is the inertia-relief equilibrium of a free body), so the push produces pure bending of the prescribed amount. This is a modelling choice: a sheet hanging freely from a thread and pushed 10 mm at a corner would mostly swing as a pendulum and bend by only about 1.5 mm.
4949
5050 ## 3. Excitation
5151
@@ -69,7 +69,9 @@ The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $
6969
7070 **Mass.** The mass matrix is lumped with trapezoidal weights: a full cell mass in the interior, one half on an edge, one quarter at a corner, so that the total mass is exact.
7171
72-**Validation.** The tests in `tests/` check that both stiffness blocks are symmetric with exactly three rigid modes, that the forces are the exact gradient of the discrete energy at a strongly nonlinear state, and that the first five frequencies of a completely free square plate agree with the reference values of Leissa and Gorman ($\lambda^2 = 13.47, 19.60, 24.27, 34.80, 34.80$) to within 1 percent at 20 cells per side.
72+**Elliptical sheets.** The same grid serves an elliptical sheet through a mask. A cell belongs to the sheet when its centre lies inside the ellipse, and a node when it is a corner of such a cell; the edge is therefore a staircase of whole cells. Every rule above carries over unchanged: a second difference exists where the node and both of its neighbours belong to the sheet, a cell contributes membrane energy and the mixed bending term when its four corners do, and the weight of a node is the fraction of its four cells that belong to the sheet, which reduces to the trapezoidal weights on a rectangle. Because the forces are still the gradient of the discrete energy, the staircase edge is free, the stiffness stays symmetric, and exactly three rigid modes remain in each block. Nodes outside the sheet are frozen and carry no mass, so the rigid-body projection and the radiation sum ignore them. Deciding membership by cell centres rather than by nodes matters: a node-based mask loses about half a cell all around the edge and raises every frequency by several percent (13 percent at 24 cells across the diameter of a circle), whereas the cell-centre mask keeps the area unbiased. The hole and the pushed patch are placed on the rim in the upper-left and lower-right directions, so the push remains asymmetric as it is at the corner of a rectangle. The one-parameter description of the shape by its bounding box is convenient in the interface, since width and height keep their meaning.
73+
74+**Validation.** The tests in `tests/` check that both stiffness blocks are symmetric with exactly three rigid modes, that the forces are the exact gradient of the discrete energy at a strongly nonlinear state, and that the first five frequencies of a completely free square plate agree with the reference values of Leissa and Gorman ($\lambda^2 = 13.47, 19.60, 24.27, 34.80, 34.80$) to within 1 percent at 20 cells per side. The same checks are repeated on an elliptical mask, and the first nine frequencies of a completely free circular plate are compared with Leissa's values for $\nu = 0.3$ ($\lambda^2 = 5.36, 9.00, 12.44, 20.48, 21.84$, the modes with nodal diameters being double). The staircase edge is first-order accurate: the frequencies are within about 2 percent at 24 cells across the diameter and within 1 percent at 60, and the degenerate pairs, which the staircase splits, differ by 1 to 2 percent.
7375
7476 ## 5. Time integration
7577
@@ -93,4 +95,4 @@ The pressure is reported in pascals, and the peak sound pressure level is displa
9395
9496 ## 7. Limitations
9597
96-The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The Marguerre treatment of the initial shape is likewise a shallow-shell approximation, adequate for a bow of a few millimetres, and the real shape of a given sheet is of course unknown; the three offered shapes are idealizations. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.
98+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The Marguerre treatment of the initial shape is likewise a shallow-shell approximation, adequate for a bow of a few millimetres, and the real shape of a given sheet is of course unknown; the three offered shapes are idealizations. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate, and for an elliptical sheet the edge is a staircase of whole cells, which places the low frequencies within a few percent at the default 8 mm spacing and splits modes that would be degenerate on a true ellipse. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.
src/cli/main.tsmodified+11−6View file
@@ -15,7 +15,7 @@ import { encodeWavBytes } from "../audio/wav.ts";
1515 import { Spectrogram, spectrogramPixels } from "../audio/spectrogram.ts";
1616 import { encodePng } from "./png.ts";
1717 import { clipMetrics, type ClipMetrics } from "./metrics.ts";
18-import { defaultParams, type PhysParams } from "../params.ts";
18+import { defaultParams, holePoint, pushPoint, type PhysParams } from "../params.ts";
1919 import { peakDbfs } from "../audio/level.ts";
2020
2121 const args = process.argv.slice(2);
@@ -207,9 +207,14 @@ function setupDiagram(): string {
207207 const p = defaultParams;
208208 const k = 0.5; // px per mm
209209 const x0 = 60, y0 = 80, W = p.Lx * 1e3 * k, H = p.Ly * 1e3 * k;
210- const hole = [x0 + p.holeInset * 1e3 * k + 4, y0 + p.holeInset * 1e3 * k + 4];
211- const push = [x0 + W, y0 + H];
210+ // Plate coordinates (m, origin at the center) to SVG pixels.
211+ const toSvg = ([X, Y]: [number, number]) => [x0 + (X + p.Lx / 2) * 1e3 * k, y0 + (p.Ly / 2 - Y) * 1e3 * k];
212+ const hole = toSvg(holePoint(p));
213+ const push = toSvg(pushPoint(p));
212214 const strike = [x0 + p.strikeFracX * W, y0 + (1 - p.strikeFracY) * H];
215+ const outline = p.shape === "ellipse"
216+ ? `<ellipse cx="${x0 + W / 2}" cy="${y0 + H / 2}" rx="${W / 2}" ry="${H / 2}" fill="#2a3340" stroke="${"#c8cbd0"}" stroke-width="1.5"/>`
217+ : `<rect x="${x0}" y="${y0}" width="${W}" height="${H}" fill="#2a3340" stroke="${"#c8cbd0"}" stroke-width="1.5"/>`;
213218 const gray = "#c8cbd0", dim = "#9aa0a8", accent = "#f0b429", blue = "#6fa8ea";
214219 const tx = 430, ty = y0 + 10, tw = W, sag = Math.max(10, p.bowAmplitude * 1e3 * k * 6);
215220 const listenerY = ty + 200;
@@ -221,7 +226,7 @@ function setupDiagram(): string {
221226 <text x="${x0}" y="${y0 - 24}" fill="${dim}">${(p.Lx / 0.0254).toFixed(0)} x ${(p.Ly / 0.0254).toFixed(0)} in, ${(p.h * 1e3).toFixed(2)} mm stainless steel</text>
222227 <!-- thread and hanging sheet -->
223228 <line x1="${hole[0]}" y1="${y0 - 8}" x2="${hole[0]}" y2="${hole[1]}" stroke="${gray}" stroke-width="1" stroke-dasharray="3 2"/>
224- <rect x="${x0}" y="${y0}" width="${W}" height="${H}" fill="#2a3340" stroke="${gray}" stroke-width="1.5"/>
229+ ${outline}
225230 <circle cx="${hole[0]}" cy="${hole[1]}" r="3" fill="#15171a" stroke="${gray}"/>
226231 <text x="${hole[0] + 10}" y="${hole[1] + 16}" fill="${dim}">hangs from a thread through a hole</text>
227232 <!-- strike point -->
@@ -231,8 +236,8 @@ function setupDiagram(): string {
231236 <text x="${strike[0] + 68}" y="${strike[1] - 46}" fill="${dim}">${(p.strikeFracX * 100).toFixed(0)}% of width, ${(p.strikeFracY * 100).toFixed(0)}% of height</text>
232237 <!-- pushed corner -->
233238 <circle cx="${push[0]}" cy="${push[1]}" r="7" fill="none" stroke="${accent}" stroke-width="2"/>
234- <text x="${x0}" y="${push[1] + 24}" fill="${accent}">finger pushes this corner ${(p.pushAmount * 1e3).toFixed(0)} mm toward the listener,</text>
235- <text x="${x0}" y="${push[1] + 38}" fill="${dim}">holds it until the sheet settles, then lets go</text>
239+ <text x="${x0}" y="${y0 + H + 24}" fill="${accent}">finger pushes this ${p.shape === "ellipse" ? "edge point" : "corner"} ${(p.pushAmount * 1e3).toFixed(0)} mm toward the listener,</text>
240+ <text x="${x0}" y="${y0 + H + 38}" fill="${dim}">holds it until the sheet settles, then lets go</text>
236241 <!-- top view -->
237242 <text x="${tx}" y="${y0 - 42}" font-weight="600">Top view</text>
238243 <text x="${tx}" y="${y0 - 24}" fill="${dim}">the sheet seen from above</text>
src/main.tsmodified+43−5View file
@@ -1,6 +1,6 @@
11 // UI wiring: parameters, run loop, live and movie views, playback, downloads.
22
3-import { defaultParams, makeGrid, makeMaterial, t60, type BowShape, type Excitation, type PhysParams } from "./params.ts";
3+import { countInside, defaultParams, makeGrid, makeMaterial, t60, type BowShape, type Excitation, type PhysParams, type PlateShape } from "./params.ts";
44 import { Simulation } from "./sim/simulation.ts";
55 import { Renderer, type ViewMode } from "./render/renderer.ts";
66 import { Spectrogram } from "./audio/spectrogram.ts";
@@ -38,7 +38,7 @@ const groups: { title: string; open: boolean; fields: FieldSpec[] }[] = [
3838 title: "Excitation: push and release",
3939 open: true,
4040 fields: [
41- { key: "pushAmount", label: "Push at corner", unit: "mm", scale: 1e-3, step: 0.5 },
41+ { key: "pushAmount", label: "Push amount", unit: "mm", scale: 1e-3, step: 0.5 },
4242 { key: "fingerRadius", label: "Finger radius", unit: "mm", scale: 1e-3, step: 1 },
4343 { key: "rampTime", label: "Push ramp", unit: "s", scale: 1, step: 0.01 },
4444 { key: "holdMax", label: "Max settle time", unit: "s", scale: 1, step: 0.1 },
@@ -122,6 +122,23 @@ function buildForm(): void {
122122 sum.textContent = g.title;
123123 det.appendChild(sum);
124124 for (const f of g.fields) {
125+ if (f.key === "Lx") {
126+ const r0 = document.createElement("div");
127+ r0.className = "field";
128+ const l0 = document.createElement("label");
129+ l0.textContent = "Shape";
130+ const sel = document.createElement("select");
131+ for (const [v, t] of [["rectangle", "Rectangle"], ["ellipse", "Ellipse (inscribed)"]]) {
132+ const o = document.createElement("option");
133+ o.value = v; o.textContent = t;
134+ sel.appendChild(o);
135+ }
136+ sel.value = params.shape;
137+ sel.addEventListener("change", () => { params.shape = sel.value as PlateShape; updateGridNote(); });
138+ r0.append(l0, sel, document.createElement("span"));
139+ det.appendChild(r0);
140+ inputs.set("shape", sel);
141+ }
125142 const row = document.createElement("div");
126143 row.className = "field";
127144 const label = document.createElement("label");
@@ -184,10 +201,24 @@ function buildForm(): void {
184201 note.id = "dampNote";
185202 det.appendChild(note);
186203 }
204+ if (f.key === "rho") {
205+ const note = document.createElement("p");
206+ note.className = "sub";
207+ note.id = "gridNote";
208+ det.appendChild(note);
209+ }
187210 }
188211 root.appendChild(det);
189212 }
190213 updateDampingNote();
214+ updateGridNote();
215+}
216+
217+function updateGridNote(): void {
218+ const g = makeGrid(params);
219+ const inside = countInside(g);
220+ const shape = params.shape === "ellipse" ? `ellipse of ${(params.Lx / 0.0254).toFixed(1)} by ${(params.Ly / 0.0254).toFixed(1)} in, edge stepped in whole cells` : "rectangle";
221+ $("gridNote").textContent = `${shape}; ${g.nx} × ${g.ny} grid, ${inside} nodes in the sheet.`;
191222 }
192223
193224 function fmt(v: number): string {
@@ -201,6 +232,7 @@ function onFieldChange(f: FieldSpec, inp: HTMLInputElement): void {
201232 const v = parseFloat(inp.value);
202233 if (!isFinite(v)) return;
203234 (params as unknown as Record<string, number>)[f.key] = v * f.scale;
235+ if (f.key === "Lx" || f.key === "Ly" || f.key === "spacing") updateGridNote();
204236 if (f.key === "sigma0" || f.key === "sigma1" || f.key === "sigmaIn") {
205237 updateDampingNote();
206238 sim?.setDamping(params.sigma0, params.sigma1, params.sigmaIn);
@@ -222,6 +254,7 @@ function setFormEnabled(enabled: boolean): void {
222254 for (const g of groups) for (const f of g.fields) if (!f.live) inputs.get(f.key)!.disabled = !enabled;
223255 inputs.get("bowShape")!.disabled = !enabled;
224256 inputs.get("excitation")!.disabled = !enabled;
257+ inputs.get("shape")!.disabled = !enabled;
225258 }
226259
227260 // ---------- State ----------
@@ -485,7 +518,7 @@ async function start(): Promise<void> {
485518 btnPause.textContent = "Pause";
486519 params.duration = duration;
487520 sim = new Simulation(device, params);
488- renderer = new Renderer(device, $("view"), sim.grid, params, sim.w0);
521+ renderer = new Renderer(device, $("view"), sim.grid, params, sim.w0, sim.nw);
489522 movieBuf?.destroy();
490523 movieBuf = device.createBuffer({ size: Math.max(16, sim.grid.n * 4), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
491524 movieBufFrame = -1;
@@ -600,7 +633,7 @@ function updateInfo(st: ReturnType<Simulation["status"]>): void {
600633 const rows: [string, string][] = [
601634 ["Phase", phase],
602635 ["Time", `${st.t.toFixed(3)} s`],
603- ["Grid", `${g.nx} × ${g.ny} nodes, ${(g.hx * 1e3).toFixed(2)} × ${(g.hy * 1e3).toFixed(2)} mm`],
636+ ["Grid", `${g.nx} × ${g.ny} nodes${g.mask ? ` (${countInside(g)} in the sheet)` : ""}, ${(g.hx * 1e3).toFixed(2)} × ${(g.hy * 1e3).toFixed(2)} mm`],
604637 ["Time step", `${(sim.dt * 1e6).toFixed(3)} µs (${sim.stepsPerSample} per sample)`],
605638 ["Speed", `${(stepsPerSec / 1e3).toFixed(0)} k steps/s (${(stepsPerSec * sim.dt).toFixed(3)}× real time)`],
606639 ["Kinetic energy", `${st.ke.toExponential(2)} J`],
@@ -652,7 +685,7 @@ async function init(): Promise<void> {
652685 function applyUrlParams(): void {
653686 const q = new URLSearchParams(location.search);
654687 for (const [k, v] of q) {
655- if (k in params && k !== "radiation" && k !== "bowShape" && k !== "excitation") {
688+ if (k in params && k !== "radiation" && k !== "bowShape" && k !== "excitation" && k !== "shape") {
656689 const num = parseFloat(v);
657690 if (isFinite(num)) {
658691 (params as unknown as Record<string, number>)[k] = num;
@@ -670,12 +703,17 @@ function applyUrlParams(): void {
670703 params.bowShape = v as BowShape;
671704 (inputs.get("bowShape") as HTMLSelectElement).value = v;
672705 }
706+ if (k === "shape" && (v === "rectangle" || v === "ellipse")) {
707+ params.shape = v;
708+ (inputs.get("shape") as HTMLSelectElement).value = v;
709+ }
673710 if (k === "duration") {
674711 duration = params.duration;
675712 durInput.value = String(duration);
676713 }
677714 }
678715 updateDampingNote();
716+ updateGridNote();
679717 (window as unknown as { steelsim: unknown }).steelsim = {
680718 status: () => (sim ? { ...sim.status(), dt: sim.dt, stepsPerSample: sim.stepsPerSample, nx: sim.grid.nx, ny: sim.grid.ny, stepsPerSec } : null),
681719 audio: () => (sim ? Array.from(sim.processedAudio()) : null),
src/params.tsmodified+108−7View file
@@ -1,7 +1,8 @@
11 // Physical and numerical parameters of the simulation, plus derived quantities.
22
33 export interface PhysParams {
4- Lx: number; // plate width along x [m]
4+ shape: PlateShape; // outline of the sheet within its bounding box
5+ Lx: number; // plate width along x [m] (bounding box for an ellipse)
56 Ly: number; // plate height along y [m]
67 h: number; // thickness [m]
78 E: number; // Young's modulus [Pa]
@@ -40,9 +41,13 @@ export interface PhysParams {
4041
4142 export type Excitation = "push" | "strike";
4243
44+/** A rectangle Lx by Ly, or the ellipse inscribed in that rectangle. */
45+export type PlateShape = "rectangle" | "ellipse";
46+
4347 export type BowShape = "none" | "width" | "height" | "dish";
4448
4549 export const defaultParams: PhysParams = {
50+ shape: "rectangle",
4651 Lx: 12 * 0.0254,
4752 Ly: 18 * 0.0254,
4853 h: 0.79e-3,
@@ -85,6 +90,28 @@ export function strikePoint(p: PhysParams): [number, number] {
8590 return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
8691 }
8792
93+/**
94+ * Center of the thread hole. Rectangle: inset from the two edges of the
95+ * top-left corner. Ellipse: inset along the radius from the boundary point
96+ * in the top-left direction (parametric angle 135 degrees).
97+ */
98+export function holePoint(p: PhysParams): [number, number] {
99+ if (p.shape === "ellipse") {
100+ const c = Math.SQRT1_2;
101+ return [-(p.Lx / 2 - p.holeInset) * c, (p.Ly / 2 - p.holeInset) * c];
102+ }
103+ return [-p.Lx / 2 + p.holeInset, p.Ly / 2 - p.holeInset];
104+}
105+
106+/** Center of the pushed patch: the bottom-right corner, or the boundary point opposite the hole. */
107+export function pushPoint(p: PhysParams): [number, number] {
108+ if (p.shape === "ellipse") {
109+ const c = Math.SQRT1_2;
110+ return [(p.Lx / 2) * c, -(p.Ly / 2) * c];
111+ }
112+ return [p.Lx / 2, -p.Ly / 2];
113+}
114+
88115 /** Stress-free initial shape w0 at every node. */
89116 export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
90117 const w0 = new Float32Array(g.n);
@@ -115,6 +142,8 @@ export interface Grid {
115142 hx: number;
116143 hy: number;
117144 n: number;
145+ /** 1 for nodes that belong to the sheet, 0 otherwise; absent for a full rectangle. */
146+ mask?: Uint8Array;
118147 }
119148
120149 export interface Material {
@@ -128,7 +157,63 @@ export interface Material {
128157 export function makeGrid(p: PhysParams): Grid {
129158 const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
130159 const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
131- return { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
160+ const g: Grid = { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
161+ return p.shape === "ellipse" ? withEllipseMask(g) : g;
162+}
163+
164+/**
165+ * Restricts a grid to the ellipse inscribed in its bounding box. The edge is a
166+ * staircase of whole cells: a cell belongs to the sheet when its center lies
167+ * inside the ellipse, and a node when it is a corner of such a cell. Deciding
168+ * by cell centers rather than by nodes keeps the covered area unbiased (a
169+ * node-based mask loses about half a cell all around the edge, which raises
170+ * every frequency by a few percent). Nodes without a cell of their own would
171+ * carry no mass; the trailing loop removes any such node.
172+ */
173+export function withEllipseMask(g: Grid): Grid {
174+ const mask = new Uint8Array(g.n);
175+ const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
176+ for (let cj = 0; cj < g.ny - 1; cj++) {
177+ for (let ci = 0; ci < g.nx - 1; ci++) {
178+ const X = (nodeX(g, ci) + g.hx / 2) / a, Y = (nodeY(g, cj) + g.hy / 2) / b;
179+ if (X * X + Y * Y <= 1) {
180+ mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
181+ mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
182+ }
183+ }
184+ }
185+ const gm: Grid = { ...g, mask };
186+ for (let changed = true; changed; ) {
187+ changed = false;
188+ for (let j = 0; j < g.ny; j++) {
189+ for (let i = 0; i < g.nx; i++) {
190+ if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
191+ mask[j * g.nx + i] = 0;
192+ changed = true;
193+ }
194+ }
195+ }
196+ }
197+ return gm;
198+}
199+
200+/** True when node (i, j) exists and belongs to the sheet. */
201+export function insideNode(g: Grid, i: number, j: number): boolean {
202+ if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
203+ return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
204+}
205+
206+/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the sheet. */
207+export function insideCell(g: Grid, ci: number, cj: number): boolean {
208+ return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
209+}
210+
211+/** Number of inside nodes (the degrees of freedom of the sheet). */
212+export function countInside(g: Grid): number {
213+ if (!g.mask) return g.n;
214+ let c = 0;
215+ for (let k = 0; k < g.n; k++) c += g.mask[k];
216+ return c;
132217 }
133218
134219 export function makeMaterial(p: PhysParams): Material {
@@ -149,11 +234,27 @@ export function nodeY(g: Grid, j: number): number {
149234 return (j - (g.ny - 1) / 2) * g.hy;
150235 }
151236
152-/** Trapezoidal weight of a node: 1 in the interior, 1/2 on an edge, 1/4 at a corner. */
153-export function trapWeight(g: Grid, i: number, j: number): number {
154- const tx = i === 0 || i === g.nx - 1 ? 0.5 : 1;
155- const ty = j === 0 || j === g.ny - 1 ? 0.5 : 1;
156- return tx * ty;
237+/**
238+ * Quadrature (and lumped-mass) weight of a node: the fraction of its four
239+ * surrounding cells that belong to the sheet. On a full rectangle this is the
240+ * trapezoidal rule: 1 in the interior, 1/2 on an edge, 1/4 at a corner. Zero
241+ * for nodes outside the sheet.
242+ */
243+export function nodeWeight(g: Grid, i: number, j: number): number {
244+ if (!insideNode(g, i, j)) return 0;
245+ let c = 0;
246+ if (insideCell(g, i, j)) c++;
247+ if (insideCell(g, i - 1, j)) c++;
248+ if (insideCell(g, i, j - 1)) c++;
249+ if (insideCell(g, i - 1, j - 1)) c++;
250+ return c / 4;
251+}
252+
253+/** Per-node weights as a flat array (0 outside the sheet). */
254+export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
255+ const w = new Float32Array(g.n);
256+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w[j * g.nx + i] = nodeWeight(g, i, j);
257+ return w;
157258 }
158259
159260 /** Decay time T60 [s] of a mode with angular frequency omega under the two-parameter loss model. */
src/render/renderer.tsmodified+14−7View file
@@ -1,6 +1,6 @@
11 // WebGPU rendering of the displacement field: a 2D color map or a lit 3D height field.
22
3-import { type Grid, type PhysParams, strikePoint } from "../params.ts";
3+import { type Grid, type PhysParams, holePoint, insideCell, pushPoint, strikePoint } from "../params.ts";
44 import { lookAt, multiply, perspective } from "./mat4.ts";
55 import shaderCode from "./shaders.wgsl?raw";
66
@@ -32,10 +32,12 @@ export class Renderer {
3232 private bindGroups = new Map<GPUBuffer, GPUBindGroup>();
3333
3434 private w0: GPUBuffer;
35+ private nw: GPUBuffer;
3536
36- constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams, w0: GPUBuffer) {
37+ constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams, w0: GPUBuffer, nw: GPUBuffer) {
3738 this.device = device;
3839 this.w0 = w0;
40+ this.nw = nw;
3941 this.canvas = canvas;
4042 this.grid = grid;
4143 this.params = params;
@@ -49,6 +51,7 @@ export class Renderer {
4951 { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "uniform" } },
5052 { binding: 1, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
5153 { binding: 2, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
54+ { binding: 3, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
5255 ],
5356 });
5457 const layout = device.createPipelineLayout({ bindGroupLayouts: [this.bgLayout] });
@@ -65,19 +68,20 @@ export class Renderer {
6568 primitive: { topology: "triangle-list", cullMode: "none" },
6669 depthStencil: { format: "depth24plus", depthWriteEnabled: true, depthCompare: "less" },
6770 });
68- // Triangle indices of the height field.
71+ // Triangle indices of the height field (cells of the sheet only).
6972 const { nx, ny } = grid;
7073 const idx = new Uint32Array((nx - 1) * (ny - 1) * 6);
7174 let k = 0;
7275 for (let j = 0; j < ny - 1; j++) {
7376 for (let i = 0; i < nx - 1; i++) {
77+ if (!insideCell(grid, i, j)) continue;
7478 const a = j * nx + i, b = a + 1, c = a + nx, d = c + 1;
7579 idx[k++] = a; idx[k++] = b; idx[k++] = d;
7680 idx[k++] = a; idx[k++] = d; idx[k++] = c;
7781 }
7882 }
79- this.indexCount = idx.length;
80- this.index = device.createBuffer({ size: idx.byteLength, usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
83+ this.indexCount = k;
84+ this.index = device.createBuffer({ size: Math.max(16, idx.byteLength), usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
8185 device.queue.writeBuffer(this.index, 0, idx);
8286 this.installControls();
8387 }
@@ -134,15 +138,17 @@ export class Renderer {
134138 u[20] = this.exaggeration; u[21] = range; u[22] = aspect;
135139 ui[23] = this.mode === "3d" ? 1 : 0;
136140 const p = this.params;
137- u[24] = -p.Lx / 2 + p.holeInset; u[25] = p.Ly / 2 - p.holeInset;
141+ [u[24], u[25]] = holePoint(p);
138142 if (p.excitation === "strike") {
139143 const [sx, sy] = strikePoint(p);
140144 u[26] = sx; u[27] = sy; u[28] = p.contactRadius;
141145 } else {
142- u[26] = p.Lx / 2; u[27] = -p.Ly / 2; u[28] = p.fingerRadius;
146+ [u[26], u[27]] = pushPoint(p);
147+ u[28] = p.fingerRadius;
143148 }
144149 u[29] = 0.4; u[30] = 0.5; u[31] = 1.0;
145150 ui[32] = this.showShape ? 1 : 0;
151+ ui[33] = p.shape === "ellipse" ? 1 : 0;
146152 this.device.queue.writeBuffer(this.uniform, 0, u);
147153
148154 let bg = this.bindGroups.get(wBuffer);
@@ -153,6 +159,7 @@ export class Renderer {
153159 { binding: 0, resource: { buffer: this.uniform } },
154160 { binding: 1, resource: { buffer: wBuffer } },
155161 { binding: 2, resource: { buffer: this.w0 } },
162+ { binding: 3, resource: { buffer: this.nw } },
156163 ],
157164 });
158165 this.bindGroups.set(wBuffer, bg);
src/render/shaders.wgslmodified+28−4View file
@@ -17,7 +17,7 @@ struct RU {
1717 lightY: f32,
1818 lightZ: f32,
1919 showShape: u32, // 3D: add the stress-free initial shape to the height field
20- pad0: u32,
20+ shape: u32, // 0: rectangle, 1: ellipse inscribed in the bounding box
2121 pad1: u32,
2222 pad2: u32,
2323 };
@@ -25,6 +25,7 @@ struct RU {
2525 @group(0) @binding(0) var<uniform> R: RU;
2626 @group(0) @binding(1) var<storage, read> w: array<f32>;
2727 @group(0) @binding(2) var<storage, read> w0: array<f32>;
28+@group(0) @binding(3) var<storage, read> nw: array<f32>; // node weights, 0 outside the sheet
2829
2930 fn colormap(t: f32) -> vec3f {
3031 let s = clamp(t, -1.0, 1.0);
@@ -40,6 +41,10 @@ fn wAt(i: i32, j: i32) -> f32 {
4041 let jj = clamp(j, 0, i32(R.ny) - 1);
4142 return w[u32(jj) * R.nx + u32(ii)];
4243 }
44+fn nwAt(i: i32, j: i32) -> f32 {
45+ if (i < 0 || i >= i32(R.nx) || j < 0 || j >= i32(R.ny)) { return 0.0; }
46+ return nw[u32(j) * R.nx + u32(i)];
47+}
4348 // Height used by the 3D view: w, plus the initial shape when requested.
4449 fn zAt(i: i32, j: i32) -> f32 {
4550 let ii = clamp(i, 0, i32(R.nx) - 1);
@@ -48,6 +53,15 @@ fn zAt(i: i32, j: i32) -> f32 {
4853 if (R.showShape != 0u) { z += w0[u32(jj) * R.nx + u32(ii)]; }
4954 return z;
5055 }
56+// Slope along one axis, one-sided where a neighbour lies outside the sheet.
57+fn slope(i: i32, j: i32, di: i32, dj: i32, h: f32) -> f32 {
58+ let hasP = nwAt(i + di, j + dj) > 0.0;
59+ let hasM = nwAt(i - di, j - dj) > 0.0;
60+ if (hasP && hasM) { return (zAt(i + di, j + dj) - zAt(i - di, j - dj)) / (2.0 * h); }
61+ if (hasP) { return (zAt(i + di, j + dj) - zAt(i, j)) / h; }
62+ if (hasM) { return (zAt(i, j) - zAt(i - di, j - dj)) / h; }
63+ return 0.0;
64+}
5165
5266 // ---------- 2D view: full-screen triangle, bilinear lookup in the fragment shader ----------
5367
@@ -75,11 +89,21 @@ fn fs2d(in: V2) -> @location(0) vec4f {
7589 let py = in.uv.y / sy * 0.5 * Ly;
7690 let bg = vec4f(0.11, 0.12, 0.14, 1.0);
7791 if (abs(px) > 0.5 * Lx || abs(py) > 0.5 * Ly) { return bg; }
92+ if (R.shape == 1u) {
93+ let ex = px / (0.5 * Lx); let ey = py / (0.5 * Ly);
94+ if (ex * ex + ey * ey > 1.0) { return bg; }
95+ }
7896 let fi = (px + 0.5 * Lx) / R.hx;
7997 let fj = (py + 0.5 * Ly) / R.hy;
8098 let i0 = i32(floor(fi)); let j0 = i32(floor(fj));
8199 let a = fi - f32(i0); let b = fj - f32(j0);
82- let val = mix(mix(wAt(i0, j0), wAt(i0 + 1, j0), a), mix(wAt(i0, j0 + 1), wAt(i0 + 1, j0 + 1), a), b);
100+ // Bilinear interpolation over the nodes of the sheet only (outside nodes carry no value).
101+ let c = vec4f((1.0 - a) * (1.0 - b), a * (1.0 - b), (1.0 - a) * b, a * b);
102+ let m = vec4f(nwAt(i0, j0), nwAt(i0 + 1, j0), nwAt(i0, j0 + 1), nwAt(i0 + 1, j0 + 1));
103+ let wgt = c * step(vec4f(1e-6), m);
104+ let den = wgt.x + wgt.y + wgt.z + wgt.w;
105+ if (den <= 0.0) { return bg; }
106+ let val = dot(wgt, vec4f(wAt(i0, j0), wAt(i0 + 1, j0), wAt(i0, j0 + 1), wAt(i0 + 1, j0 + 1))) / den;
83107 var col = colormap(val / R.colorRange);
84108 // Markers: hole (ring) and pushed patch (dashed-looking ring).
85109 let dh = distance(vec2f(px, py), vec2f(R.holeX, R.holeY));
@@ -101,8 +125,8 @@ fn vs3d(@builtin(vertex_index) vi: u32) -> V3 {
101125 let y = (f32(j) - 0.5 * f32(R.ny - 1u)) * R.hy;
102126 let val = wAt(i, j);
103127 let z = zAt(i, j) * R.wScale;
104- let dzdx = (zAt(i + 1, j) - zAt(i - 1, j)) * R.wScale / (2.0 * R.hx);
105- let dzdy = (zAt(i, j + 1) - zAt(i, j - 1)) * R.wScale / (2.0 * R.hy);
128+ let dzdx = slope(i, j, 1, 0, R.hx) * R.wScale;
129+ let dzdy = slope(i, j, 0, 1, R.hy) * R.wScale;
106130 var o: V3;
107131 o.pos = R.mvp * vec4f(x, y, z, 1.0);
108132 o.val = val;
src/sim/cpuRef.tsmodified+25−13View file
@@ -2,7 +2,7 @@
22 // It mirrors step.wgsl exactly and is used for the stable time-step estimate
33 // and for the tests. It is not used in the time loop of the browser simulation.
44
5-import { type Grid, type Material, trapWeight } from "../params.ts";
5+import { type Grid, type Material, insideCell, insideNode, nodeWeight } from "../params.ts";
66
77 export interface Fields {
88 u: Float64Array;
@@ -26,24 +26,28 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?:
2626 const nu = m.nu;
2727 const id = (i: number, j: number) => j * nx + i;
2828
29+ // A second difference exists where the node and both of its neighbours belong to the sheet.
30+ const ins = (i: number, j: number) => insideNode(g, i, j);
31+ const hasA = (i: number, j: number) => ins(i - 1, j) && ins(i, j) && ins(i + 1, j);
32+ const hasB = (i: number, j: number) => ins(i, j - 1) && ins(i, j) && ins(i, j + 1);
2933 // Bending moment-like quantities at nodes (zero where the second difference is undefined).
3034 const Ma = (i: number, j: number): number => {
31- if (i < 1 || i > nx - 2 || j < 0 || j > ny - 1) return 0;
35+ if (!hasA(i, j)) return 0;
3236 const a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
3337 let b = 0;
34- if (j >= 1 && j <= ny - 2) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
35- return trapWeight(g, i, j) * (a + nu * b);
38+ if (hasB(i, j)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
39+ return nodeWeight(g, i, j) * (a + nu * b);
3640 };
3741 const Mb = (i: number, j: number): number => {
38- if (j < 1 || j > ny - 2 || i < 0 || i > nx - 1) return 0;
42+ if (!hasB(i, j)) return 0;
3943 const b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
4044 let a = 0;
41- if (i >= 1 && i <= nx - 2) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
42- return trapWeight(g, i, j) * (b + nu * a);
45+ if (hasA(i, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
46+ return nodeWeight(g, i, j) * (b + nu * a);
4347 };
4448 // Mixed second difference at a cell center (zero for cells outside the plate).
4549 const Cc = (ci: number, cj: number): number => {
46- if (ci < 0 || ci > nx - 2 || cj < 0 || cj > ny - 2) return 0;
50+ if (!insideCell(g, ci, cj)) return 0;
4751 return (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
4852 };
4953
@@ -57,7 +61,7 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?:
5761
5862 // Membrane force on local node `a` of cell (ci, cj); returns [Fu, Fv, Fw].
5963 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
60- if (ci < 0 || ci > nx - 2 || cj < 0 || cj > ny - 2) return;
64+ if (!insideCell(g, ci, cj)) return;
6165 const n0 = id(ci, cj), n1 = id(ci + 1, cj), n2 = id(ci + 1, cj + 1), n3 = id(ci, cj + 1);
6266 U[0] = u[n0]; U[1] = u[n1]; U[2] = u[n2]; U[3] = u[n3];
6367 V[0] = v[n0]; V[1] = v[n1]; V[2] = v[n2]; V[3] = v[n3];
@@ -98,6 +102,11 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?:
98102 for (let j = 0; j < ny; j++) {
99103 for (let i = 0; i < nx; i++) {
100104 const k = id(i, j);
105+ const weight = nodeWeight(g, i, j);
106+ if (weight === 0) {
107+ out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
108+ continue;
109+ }
101110 const bend =
102111 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
103112 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
@@ -107,7 +116,7 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?:
107116 cellForce(i - 1, j, 1, F);
108117 cellForce(i - 1, j - 1, 2, F);
109118 cellForce(i, j - 1, 3, F);
110- const mass = m.rhoh * hx * hy * trapWeight(g, i, j);
119+ const mass = m.rhoh * hx * hy * weight;
111120 out.u[k] = F[0] / mass;
112121 out.v[k] = F[1] / mass;
113122 out.w[k] = F[2] / mass;
@@ -159,17 +168,19 @@ export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Arr
159168 const { u, v, w } = f;
160169 const nu = m.nu;
161170 const id = (i: number, j: number) => j * nx + i;
171+ const ins = (i: number, j: number) => insideNode(g, i, j);
162172 let Ub = 0;
163173 for (let j = 0; j < ny; j++) {
164174 for (let i = 0; i < nx; i++) {
165175 let a = 0, b = 0;
166- if (i >= 1 && i <= nx - 2) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
167- if (j >= 1 && j <= ny - 2) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
168- Ub += trapWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
176+ if (ins(i - 1, j) && ins(i, j) && ins(i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
177+ if (ins(i, j - 1) && ins(i, j) && ins(i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
178+ Ub += nodeWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
169179 }
170180 }
171181 for (let cj = 0; cj < ny - 1; cj++) {
172182 for (let ci = 0; ci < nx - 1; ci++) {
183+ if (!insideCell(g, ci, cj)) continue;
173184 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
174185 Ub += 2 * (1 - nu) * c * c;
175186 }
@@ -180,6 +191,7 @@ export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Arr
180191 const q = (hx * hy) / 4;
181192 for (let cj = 0; cj < ny - 1; cj++) {
182193 for (let ci = 0; ci < nx - 1; ci++) {
194+ if (!insideCell(g, ci, cj)) continue;
183195 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
184196 for (let gp = 0; gp < 4; gp++) {
185197 const xi = gp & 1 ? G3 : -G3;
src/sim/shaders/common.wgslmodified+7−5View file
@@ -49,9 +49,11 @@ struct StepParams {
4949
5050 fn nodeX(i: i32) -> f32 { return (f32(i) - 0.5 * f32(P.nx - 1u)) * P.hx; }
5151 fn nodeY(j: i32) -> f32 { return (f32(j) - 0.5 * f32(P.ny - 1u)) * P.hy; }
52-fn trapWeight(i: i32, j: i32) -> f32 {
53- var t = 1.0;
54- if (i == 0 || i == i32(P.nx) - 1) { t = t * 0.5; }
55- if (j == 0 || j == i32(P.ny) - 1) { t = t * 0.5; }
56- return t;
52+// Quadrature / lumped-mass weight of every node: the fraction of its four
53+// cells that belong to the sheet (trapezoidal weights on a full rectangle),
54+// zero for nodes outside the sheet. Written once by the host.
55+@group(0) @binding(9) var<storage, read> nw: array<f32>;
56+fn nodeWeight(i: i32, j: i32) -> f32 {
57+ if (i < 0 || i >= i32(P.nx) || j < 0 || j >= i32(P.ny)) { return 0.0; }
58+ return nw[u32(j) * P.nx + u32(i)];
5759 }
src/sim/shaders/rigid.wgslmodified+4−3View file
@@ -18,7 +18,7 @@ fn sums(@builtin(local_invocation_id) lid: vec3u) {
1818 for (var k = t; k < P.n; k += 256u) {
1919 let i = i32(k % P.nx); let j = i32(k / P.nx);
2020 let x = nodeX(i); let y = nodeY(j);
21- let m = P.rhoh * P.hx * P.hy * trapWeight(i, j);
21+ let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
2222 let n = P.n;
2323 let p = vec3f(pos[k], pos[n + k], pos[2u * n + k]);
2424 let v = vec3f(vel[k], vel[n + k], vel[2u * n + k]);
@@ -57,6 +57,7 @@ fn apply(@builtin(global_invocation_id) gid: vec3u) {
5757 let k = gid.x;
5858 if (k >= P.n) { return; }
5959 let i = i32(k % P.nx); let j = i32(k / P.nx);
60+ if (nodeWeight(i, j) == 0.0) { return; } // outside the sheet
6061 let x = nodeX(i); let y = nodeY(j);
6162 let n = P.n;
6263 pos[k] -= stats[0] - stats[3] * y;
@@ -78,7 +79,7 @@ fn projectField(@builtin(local_invocation_id) lid: vec3u) {
7879 var s = vec3f(0.0);
7980 for (var k = t; k < P.n; k += 256u) {
8081 let i = i32(k % P.nx); let j = i32(k / P.nx);
81- let m = P.rhoh * P.hx * P.hy * trapWeight(i, j);
82+ let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
8283 let f = field[k];
8384 s += vec3f(m * f, m * nodeY(j) * f, -m * nodeX(i) * f);
8485 }
@@ -91,6 +92,6 @@ fn projectField(@builtin(local_invocation_id) lid: vec3u) {
9192 let c = vec3f(sh3[0].x / P.mass, sh3[0].y / P.Ix, sh3[0].z / P.Iy);
9293 for (var k = t; k < P.n; k += 256u) {
9394 let i = i32(k % P.nx); let j = i32(k / P.nx);
94- field[k] -= c.x + c.y * nodeY(j) - c.z * nodeX(i);
95+ if (nodeWeight(i, j) > 0.0) { field[k] -= c.x + c.y * nodeY(j) - c.z * nodeX(i); }
9596 }
9697 }
src/sim/shaders/step.wgslmodified+47−20View file
@@ -3,8 +3,12 @@
33 // State layout (structure of arrays, n = nx * ny nodes):
44 // pos: [u | v | w | wLo] wLo carries the rounding error of the w update
55 // vel: [du | dv | dw]
6-// Each 16x16 workgroup loads its tile plus a 2-node halo of u, v, w, dw into
7-// shared memory once, so the stencils below never touch global memory.
6+// Each 16x16 workgroup loads its tile plus a 2-node halo of w and of the node
7+// weights into shared memory once, so the stencils below never touch global
8+// memory. Nodes outside the sheet (weight 0) are skipped: a second difference
9+// exists only where the node and both neighbours belong to the sheet, and a
10+// cell contributes only when its four corners do. That reproduces the free
11+// edge of the rectangle and gives a staircase free edge for an ellipse.
812
913 @group(0) @binding(0) var<uniform> P: Params;
1014 @group(0) @binding(1) var<uniform> SP: StepParams;
@@ -22,6 +26,9 @@ const TILE = 16;
2226 const HALO = 2;
2327 const TW = 20; // TILE + 2 * HALO
2428 var<workgroup> tw: array<f32, 400>;
29+// Node weights of the same tile in quarters (0..4), packed four to a word to stay
30+// within the default 16 KB of workgroup storage.
31+var<workgroup> tn: array<u32, 100>;
2532 // Membrane force contributions of the 17x17 cells around the tile, 4 nodes x 3 components each.
2633 const CT = 17;
2734 var<workgroup> cellF: array<f32, 3468>;
@@ -37,6 +44,10 @@ fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
3744 fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
3845 // Global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
3946 fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
47+fn NWq(i: i32, j: i32) -> u32 { let l = lidx(i, j); return (tn[l >> 2u] >> ((l & 3u) * 8u)) & 0xFFu; }
48+fn NW(i: i32, j: i32) -> f32 { return 0.25 * f32(NWq(i, j)); }
49+fn INS(i: i32, j: i32) -> bool { return NWq(i, j) != 0u; }
50+fn insCell(ci: i32, cj: i32) -> bool { return INS(ci, cj) && INS(ci + 1, cj) && INS(ci + 1, cj + 1) && INS(ci, cj + 1); }
4051 fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }
4152 fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
4253 fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
@@ -44,24 +55,22 @@ fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }
4455
4556 // Weighted bending "moments" at nodes; zero where the second difference is undefined.
4657 fn Ma(i: i32, j: i32) -> f32 {
47- let nx = i32(P.nx); let ny = i32(P.ny);
48- if (i < 1 || i > nx - 2 || j < 0 || j > ny - 1) { return 0.0; }
58+ if (!(INS(i - 1, j) && INS(i, j) && INS(i + 1, j))) { return 0.0; }
4959 let a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2;
5060 var b = 0.0;
51- if (j >= 1 && j <= ny - 2) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
52- return trapWeight(i, j) * (a + P.nu * b);
61+ if (INS(i, j - 1) && INS(i, j + 1)) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
62+ return NW(i, j) * (a + P.nu * b);
5363 }
5464 fn Mb(i: i32, j: i32) -> f32 {
55- let nx = i32(P.nx); let ny = i32(P.ny);
56- if (j < 1 || j > ny - 2 || i < 0 || i > nx - 1) { return 0.0; }
65+ if (!(INS(i, j - 1) && INS(i, j) && INS(i, j + 1))) { return 0.0; }
5766 let b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2;
5867 var a = 0.0;
59- if (i >= 1 && i <= nx - 2) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
60- return trapWeight(i, j) * (b + P.nu * a);
68+ if (INS(i - 1, j) && INS(i + 1, j)) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
69+ return NW(i, j) * (b + P.nu * a);
6170 }
62-// Mixed second difference at a cell center; zero outside the plate.
71+// Mixed second difference at a cell center; zero outside the sheet.
6372 fn Cc(ci: i32, cj: i32) -> f32 {
64- if (ci < 0 || ci > i32(P.nx) - 2 || cj < 0 || cj > i32(P.ny) - 2) { return 0.0; }
73+ if (!insCell(ci, cj)) { return 0.0; }
6574 return (W(ci + 1, cj + 1) - W(ci + 1, cj) - W(ci, cj + 1) + W(ci, cj)) * ihxy;
6675 }
6776
@@ -69,7 +78,7 @@ fn Cc(ci: i32, cj: i32) -> f32 {
6978 // written to cellF at local cell index lc.
7079 fn cellForces(ci: i32, cj: i32, lc: u32) {
7180 let base = lc * 12u;
72- if (ci < 0 || ci > i32(P.nx) - 2 || cj < 0 || cj > i32(P.ny) - 2) {
81+ if (!insCell(ci, cj)) {
7382 for (var q = 0u; q < 12u; q++) { cellF[base + q] = 0.0; }
7483 return;
7584 }
@@ -121,7 +130,7 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
121130 ox = i32(wid.x) * TILE - HALO;
122131 oy = i32(wid.y) * TILE - HALO;
123132 ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
124- // Cooperative load of the w tile (400 entries, 256 threads).
133+ // Cooperative load of the w and weight tiles (400 entries each, 256 threads).
125134 let lt = lid.y * 16u + lid.x;
126135 for (var e = lt; e < 400u; e += 256u) {
127136 let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
@@ -129,6 +138,17 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
129138 if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { ww = posIn[2u * n + id(gi, gj)]; }
130139 tw[e] = ww;
131140 }
141+ if (lt < 100u) {
142+ var word = 0u;
143+ for (var q = 0u; q < 4u; q++) {
144+ let e = lt * 4u + q;
145+ let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
146+ var wq = 0u;
147+ if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { wq = u32(round(4.0 * nw[id(gi, gj)])); }
148+ word |= wq << (q * 8u);
149+ }
150+ tn[lt] = word;
151+ }
132152 workgroupBarrier();
133153 // Membrane forces of the 17x17 cells adjacent to the tile's nodes.
134154 for (var e = lt; e < 289u; e += 256u) {
@@ -140,7 +160,14 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
140160 let i = i32(gid.x); let j = i32(gid.y);
141161 if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
142162 let k = id(i, j);
143-
163+ let weight = NW(i, j);
164+ if (weight == 0.0) {
165+ // Outside the sheet: the node stays at rest.
166+ posOut[k] = 0.0; posOut[n + k] = 0.0; posOut[2u * n + k] = 0.0; posOut[3u * n + k] = 0.0;
167+ velOut[k] = 0.0; velOut[n + k] = 0.0; velOut[2u * n + k] = 0.0;
168+ if ((SP.flags & 2u) != 0u) { react[k] = 0.0; }
169+ return;
170+ }
144171
145172 let bend = (Ma(i - 1, j) - 2.0 * Ma(i, j) + Ma(i + 1, j)) * ihx2
146173 + (Mb(i, j - 1) - 2.0 * Mb(i, j) + Mb(i, j + 1)) * ihy2
@@ -150,16 +177,16 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
150177 F += cellForce(i - 1, j, 1u);
151178 F += cellForce(i - 1, j - 1, 2u);
152179 F += cellForce(i, j - 1, 3u);
153- let mass = P.rhoh * P.hx * P.hy * trapWeight(i, j);
180+ let mass = P.rhoh * P.hx * P.hy * weight;
154181 var acc = F * (1.0 / mass);
155182
156183 // Loss: frequency-independent on all components, frequency-dependent (graph Laplacian) on w.
157184 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
158185 var lap = 0.0;
159- if (i > 0) { lap += (DW(i - 1, j) - vel.z) * ihx2; }
160- if (i < i32(P.nx) - 1) { lap += (DW(i + 1, j) - vel.z) * ihx2; }
161- if (j > 0) { lap += (DW(i, j - 1) - vel.z) * ihy2; }
162- if (j < i32(P.ny) - 1) { lap += (DW(i, j + 1) - vel.z) * ihy2; }
186+ if (INS(i - 1, j)) { lap += (DW(i - 1, j) - vel.z) * ihx2; }
187+ if (INS(i + 1, j)) { lap += (DW(i + 1, j) - vel.z) * ihx2; }
188+ if (INS(i, j - 1)) { lap += (DW(i, j - 1) - vel.z) * ihy2; }
189+ if (INS(i, j + 1)) { lap += (DW(i, j + 1) - vel.z) * ihy2; }
163190 acc -= 2.0 * SP.sigma0 * vel;
164191 acc.z += 2.0 * SP.sigma1 * lap;
165192 // In-plane (stretching) waves lie at several kHz; they get their own loss.
src/sim/simulation.tsmodified+24−9View file
@@ -14,9 +14,11 @@ import {
1414 makeMaterial,
1515 initialShape,
1616 strikePoint,
17+ pushPoint,
1718 nodeX,
1819 nodeY,
19- trapWeight,
20+ nodeWeight,
21+ nodeWeights,
2022 HIST_LEN,
2123 RB_INTERVAL,
2224 RB_INTERVAL_HELD,
@@ -89,6 +91,8 @@ export class Simulation {
8991 private react: GPUBuffer;
9092 /** Stress-free initial shape (f32 per node), also used by the renderer. */
9193 readonly w0: GPUBuffer;
94+ /** Node weights (f32 per node, 0 outside the sheet), also used by the renderer. */
95+ readonly nw: GPUBuffer;
9296 private lw: GPUBuffer;
9397 private audioOut: GPUBuffer;
9498 private frameBuf: GPUBuffer;
@@ -140,6 +144,8 @@ export class Simulation {
140144 this.react = f32(n + STICK_LEN);
141145 this.w0 = f32(n);
142146 device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
147+ this.nw = f32(n);
148+ device.queue.writeBuffer(this.nw, 0, nodeWeights(g));
143149 this.lw = f32(2 * n);
144150 this.audioOut = f32(this.samplesPerBatch);
145151 this.frameBuf = f32(n);
@@ -154,17 +160,18 @@ export class Simulation {
154160 let mass = 0, Ix = 0, Iy = 0;
155161 for (let j = 0; j < g.ny; j++) {
156162 for (let i = 0; i < g.nx; i++) {
157- const mk = m.rhoh * g.hx * g.hy * trapWeight(g, i, j);
163+ const mk = m.rhoh * g.hx * g.hy * nodeWeight(g, i, j);
158164 mass += mk;
159165 Ix += mk * nodeY(g, j) ** 2;
160166 Iy += mk * nodeX(g, i) ** 2;
161167 }
162168 }
163- // Drumstick contact patch.
169+ // Drumstick contact patch: the nearest node of the sheet and the sheet nodes within the tip radius.
164170 const [sx, sy] = strikePoint(params);
165171 let strikeNode = 0, best = Infinity, patchMass = 0;
166172 for (let j = 0; j < g.ny; j++) {
167173 for (let i = 0; i < g.nx; i++) {
174+ if (nodeWeight(g, i, j) === 0) continue;
168175 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
169176 if (d2 < best) { best = d2; strikeNode = j * g.nx + i; }
170177 }
@@ -172,9 +179,10 @@ export class Simulation {
172179 for (let j = 0; j < g.ny; j++) {
173180 for (let i = 0; i < g.nx; i++) {
174181 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
175- if (d2 <= params.contactRadius ** 2 || j * g.nx + i === strikeNode) patchMass += m.rhoh * g.hx * g.hy * trapWeight(g, i, j);
182+ if (d2 <= params.contactRadius ** 2 || j * g.nx + i === strikeNode) patchMass += m.rhoh * g.hx * g.hy * nodeWeight(g, i, j);
176183 }
177184 }
185+ const [px, py] = pushPoint(params);
178186 const stick = new Float32Array(STICK_LEN);
179187 stick[0] = 0; stick[1] = -params.strikeVelocity; stick[4] = 0; stick[5] = -params.strikeVelocity;
180188 device.queue.writeBuffer(this.react, n * 4, stick);
@@ -184,7 +192,7 @@ export class Simulation {
184192 pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = HIST_LEN;
185193 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
186194 pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
187- pf[12] = params.Lx / 2; pf[13] = -params.Ly / 2; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
195+ pf[12] = px; pf[13] = py; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
188196 pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = params.fs;
189197 pf[20] = sx; pf[21] = sy; pf[22] = params.contactRadius ** 2; pf[23] = patchMass;
190198 pu[24] = strikeNode; pf[25] = params.contactStiffness; pf[26] = params.contactExponent; pf[27] = params.stickMass;
@@ -217,6 +225,7 @@ export class Simulation {
217225 entry(6, "storage"),
218226 entry(7, "storage"),
219227 entry(8, "read-only-storage"),
228+ entry(9, "read-only-storage"),
220229 ],
221230 });
222231 this.stepPipe = mk(stepWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [stepLayout] }));
@@ -230,6 +239,7 @@ export class Simulation {
230239 { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
231240 { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
232241 { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
242+ { binding: 9, visibility: GPUShaderStage.COMPUTE, buffer: { type: "read-only-storage" } },
233243 ],
234244 });
235245 const rbPipeLayout = device.createPipelineLayout({ bindGroupLayouts: [rbLayout] });
@@ -239,6 +249,7 @@ export class Simulation {
239249 entries: [
240250 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },
241251 { binding: 4, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
252+ { binding: 9, visibility: GPUShaderStage.COMPUTE, buffer: { type: "read-only-storage" } },
242253 ],
243254 });
244255 this.projectPipe = device.createComputePipeline({
@@ -262,6 +273,7 @@ export class Simulation {
262273 { binding: 6, resource: { buffer: this.hist } },
263274 { binding: 7, resource: { buffer: this.react } },
264275 { binding: 8, resource: { buffer: this.w0 } },
276+ { binding: 9, resource: { buffer: this.nw } },
265277 ],
266278 }),
267279 );
@@ -273,6 +285,7 @@ export class Simulation {
273285 { binding: 1, resource: { buffer: this.pos[p] } },
274286 { binding: 2, resource: { buffer: this.vel[p] } },
275287 { binding: 3, resource: { buffer: this.statsBuf } },
288+ { binding: 9, resource: { buffer: this.nw } },
276289 ],
277290 }),
278291 );
@@ -281,6 +294,7 @@ export class Simulation {
281294 entries: [
282295 { binding: 0, resource: { buffer: this.paramsBuf } },
283296 { binding: 4, resource: { buffer: this.react } },
297+ { binding: 9, resource: { buffer: this.nw } },
284298 ],
285299 });
286300 this.audioBG = device.createBindGroup({
@@ -327,7 +341,7 @@ export class Simulation {
327341 const k = j * g.nx + i;
328342 const rx = L[0] - nodeX(g, i), ry = L[1] - nodeY(g, j), rz = L[2];
329343 const r = Math.sqrt(rx * rx + ry * ry + rz * rz);
330- const dA = g.hx * g.hy * trapWeight(g, i, j);
344+ const dA = g.hx * g.hy * nodeWeight(g, i, j);
331345 let wgt = (RHO_AIR * dA) / (2 * Math.PI * r);
332346 if (radiation === "dipole") wgt *= Math.abs(rz) / r;
333347 data[2 * k] = wgt;
@@ -424,9 +438,10 @@ export class Simulation {
424438 this.processed = b;
425439 }
426440 this.audio.set(raw, this.samples);
427- // Dipole option: first-order high-pass at the baffle cutoff c / (2 sqrt(Lx Ly)).
441+ // Dipole option: first-order high-pass at the baffle cutoff c / (2 sqrt(area)).
428442 const dipole = this.params.radiation === "dipole";
429- const fc = SPEED_OF_SOUND / (2 * Math.sqrt(this.params.Lx * this.params.Ly));
443+ const area = this.params.Lx * this.params.Ly * (this.params.shape === "ellipse" ? Math.PI / 4 : 1);
444+ const fc = SPEED_OF_SOUND / (2 * Math.sqrt(area));
430445 const alpha = 1 / (1 + (2 * Math.PI * fc) / this.params.fs);
431446 for (let k = 0; k < raw.length; k++) {
432447 const x = raw[k];
@@ -628,7 +643,7 @@ export class Simulation {
628643 }
629644
630645 destroy(): void {
631- for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
646+ for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.nw, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
632647 }
633648
634649 get fs(): number {
tests/cpuRef.test.tsmodified+100−11View file
@@ -1,7 +1,7 @@
11 import { test } from "node:test";
22 import assert from "node:assert/strict";
33 import { accelerations, potentialEnergy, maxAngularFrequency, type Fields } from "../src/sim/cpuRef.ts";
4-import { type Grid, type Material, trapWeight } from "../src/params.ts";
4+import { type Grid, type Material, countInside, nodeWeight, withEllipseMask } from "../src/params.ts";
55
66 const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
77
@@ -41,23 +41,25 @@ function symEig(A: Float64Array, n: number): number[] {
4141 return ev.sort((x, y) => x - y);
4242 }
4343
44-/** Mass-symmetrized linear stiffness of one block ("w" or "uv"). */
44+/** Mass-symmetrized linear stiffness of one block ("w" or "uv") over the nodes of the sheet. */
4545 function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
46- const N = g.n;
46+ const nodes: number[] = [];
47+ for (let k = 0; k < g.n; k++) if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) > 0) nodes.push(k);
48+ const N = nodes.length;
4749 const n = block === "w" ? N : 2 * N;
4850 const K = new Float64Array(n * n);
49- const x = fields(N), y = fields(N);
51+ const x = fields(g.n), y = fields(g.n);
5052 const eps = 1e-9;
51- const mass = (k: number) => mat.rhoh * g.hx * g.hy * trapWeight(g, k % g.nx, Math.floor(k / g.nx));
53+ const mass = (d: number) => mat.rhoh * g.hx * g.hy * nodeWeight(g, nodes[d] % g.nx, Math.floor(nodes[d] / g.nx));
5254 for (let col = 0; col < n; col++) {
5355 x.u.fill(0); x.v.fill(0); x.w.fill(0);
54- if (block === "w") x.w[col] = eps;
55- else if (col < N) x.u[col] = eps;
56- else x.v[col - N] = eps;
56+ if (block === "w") x.w[nodes[col]] = eps;
57+ else if (col < N) x.u[nodes[col]] = eps;
58+ else x.v[nodes[col - N]] = eps;
5759 accelerations(g, mat, x, y);
5860 const mc = mass(col % N);
5961 for (let row = 0; row < n; row++) {
60- const acc = block === "w" ? y.w[row] : row < N ? y.u[row] : y.v[row - N];
62+ const acc = block === "w" ? y.w[nodes[row]] : row < N ? y.u[nodes[row]] : y.v[nodes[row - N]];
6163 const mr = mass(row % N);
6264 // K_sym = M^{1/2} (M^{-1} K) M^{-1/2}
6365 K[row * n + col] = (-acc / eps) * Math.sqrt(mr / mc);
@@ -66,6 +68,43 @@ function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
6668 return { K, n };
6769 }
6870
71+function checkSymmetricWithThreeRigidModes(K: Float64Array, n: number): number[] {
72+ let asym = 0, scale = 0;
73+ for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
74+ asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
75+ scale = Math.max(scale, Math.abs(K[r * n + c]));
76+ }
77+ assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
78+ const ev = symEig(K, n);
79+ assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
80+ assert.ok(ev[3] > 1e-4 * ev[n - 1], `fourth mode is not rigid (${ev[3]} vs ${ev[n - 1]})`);
81+ return ev;
82+}
83+
84+/** Relative gradient mismatch between the forces and the finite-difference gradient of the energy. */
85+function gradientMismatch(g: Grid, f: Fields, w0?: Float64Array): number {
86+ const acc = fields(g.n);
87+ accelerations(g, mat, f, acc, w0);
88+ const delta = 1e-7;
89+ let worst = 0;
90+ for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
91+ for (let k = 0; k < g.n; k += 5) {
92+ const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
93+ if (mass === 0) continue;
94+ const orig = arr[k];
95+ arr[k] = orig + delta;
96+ const Ep = potentialEnergy(g, mat, f, w0);
97+ arr[k] = orig - delta;
98+ const Em = potentialEnergy(g, mat, f, w0);
99+ arr[k] = orig;
100+ const Ffd = -(Ep - Em) / (2 * delta);
101+ const Fan = accArr[k] * mass;
102+ worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
103+ }
104+ }
105+ return worst;
106+}
107+
69108 test("bending stiffness is symmetric and has exactly three rigid modes", () => {
70109 const nx = 9, ny = 7;
71110 const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
@@ -126,7 +165,7 @@ test("forces are the negative gradient of the discrete energy (nonlinear state)"
126165 let worst = 0;
127166 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
128167 for (let k = 0; k < g.n; k += 5) {
129- const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));
168+ const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % nx, Math.floor(k / nx));
130169 const orig = arr[k];
131170 arr[k] = orig + delta;
132171 const Ep = potentialEnergy(g, mat, f);
@@ -168,7 +207,7 @@ test("forces are the energy gradient with an initial-shape imperfection", () =>
168207 let worst = 0;
169208 for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
170209 for (let k = 0; k < g.n; k += 5) {
171- const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));
210+ const mass = mat.rhoh * g.hx * g.hy * nodeWeight(g, k % nx, Math.floor(k / nx));
172211 const orig = arr[k];
173212 arr[k] = orig + delta;
174213 const Ep = potentialEnergy(g, mat, f, w0);
@@ -185,3 +224,53 @@ test("forces are the energy gradient with an initial-shape imperfection", () =>
185224 const rest = fields(g.n);
186225 assert.ok(potentialEnergy(g, mat, rest, w0) < 1e-12, "w0 is stress free");
187226 });
227+
228+// ---------- Masked (elliptical) grids ----------
229+
230+test("elliptical sheet: bending and membrane stiffness are symmetric with exactly three rigid modes", () => {
231+ const nx = 13, ny = 10;
232+ const g = withEllipseMask({ nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny });
233+ assert.ok(countInside(g) < g.n && countInside(g) > g.n / 2, "mask removes some nodes");
234+ for (let k = 0; k < g.n; k++) {
235+ // Every node of the sheet carries mass (no whiskers left after the mask is trimmed).
236+ if (g.mask![k]) assert.ok(nodeWeight(g, k % nx, Math.floor(k / nx)) > 0, `node ${k} has no mass`);
237+ }
238+ const b = stiffness(g, "w");
239+ checkSymmetricWithThreeRigidModes(b.K, b.n);
240+ const m = stiffness(g, "uv");
241+ checkSymmetricWithThreeRigidModes(m.K, m.n);
242+});
243+
244+test("elliptical sheet: forces are the energy gradient, with and without an initial shape", () => {
245+ const nx = 12, ny = 9;
246+ const g = withEllipseMask({ nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny });
247+ const f = fields(g.n);
248+ const w0 = new Float64Array(g.n);
249+ let seed = 5;
250+ const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
251+ for (let k = 0; k < g.n; k++) {
252+ const x = (k % nx) / (nx - 1) - 0.5;
253+ w0[k] = 0.003 * (1 - 4 * x * x);
254+ f.w[k] = 0.01 * rnd();
255+ f.u[k] = 1e-4 * rnd();
256+ f.v[k] = 1e-4 * rnd();
257+ }
258+ assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
259+ assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
260+});
261+
262+test("free circular plate frequencies approach the Leissa values", () => {
263+ const nx = 25, a = 0.15; // radius; 24 cells across the diameter
264+ const g = withEllipseMask({ nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx });
265+ const { K, n } = stiffness(g, "w");
266+ const ev = symEig(K, n);
267+ const lam2 = ev.slice(3, 12).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
268+ // Completely free circular plate, nu = 0.3 (Leissa, Table 2.4): modes (2,0) x2, (0,1),
269+ // (3,0) x2, (1,1) x2, (4,0) x2 by nodal diameters and circles. The staircase edge is
270+ // first-order accurate: about 2 percent here, under 1 percent at 60 cells across.
271+ const ref = [5.358, 5.358, 9.003, 12.439, 12.439, 20.475, 20.475, 21.835, 21.835];
272+ for (let k = 0; k < ref.length; k++) {
273+ const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
274+ assert.ok(relErr < 0.04, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
275+ }
276+});