Initial shape as a list of features, with random configurations
The stress-free shape was one of three fixed parabolic bows. It is now the
sum of a list of features, each with its own parameters: a bow across the
width or the height with an apex offset, a twist, a Gaussian dent, a crease,
and a seeded smooth random field whose realization does not depend on the
grid. The mass-weighted least-squares plane is removed from every shape,
since a tilt would add a spurious membrane coupling.
The Initial shape panel edits the list and reports the peak-to-peak height
and the largest slope, warning above 0.1, where the shallow-shell treatment
becomes doubtful. Randomize draws a combination whose heights are of the
order of a given scale, so the effect of the shape on the sound can be
surveyed. Copy link puts every parameter that differs from its default into
a URL; the shape travels as a features=bow:x,0.003,0;dent:... list, and the
earlier bowShape and bowAmplitude parameters are still read.
The view shows the initial shape before a run and after any change that
takes effect at the next Start. During a run both views show the total shape
by default, the initial shape plus the deviation, since the sheet settles to
a zero deviation and would otherwise look flat; the toolbar checkbox falls
back to the deviation alone, which keeps a small vibration on a bowed sheet
visible.
tests/shape.test.ts covers the encoding, the plane removal, the geometry of
each feature, and the random configurations; the GPU self-test gains a
randomly imperfect sheet.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
14 changed files+902−105
README.mdmodified+8−6View file
@@ -17,16 +17,17 @@ Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefo
1717
1818 ## Using the simulator
1919
20-Press Start. With the drumstick excitation (the default) a free mass strikes the sheet at a configurable point and speed through a Hertz contact, and the recording of audio and movie frames starts at the impact. The Hits selector turns the single hit into a series: a fixed number of hits at a fixed interval, all at the same speed or rising or falling between the softest and the loudest speed; the Excitation panel reports the schedule and how much ring the recorded duration leaves after the last hit. With the push excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles (its kinetic energy drops below the threshold or stops falling), and released by ramping the finger force to zero over the release time; the recording starts at the release. The view shows the transverse displacement as a color map (2D) or a lit height field (3D; drag to orbit, 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.
20+Press Start. With the drumstick excitation (the default) a free mass strikes the sheet at a configurable point and speed through a Hertz contact, and the recording of audio and movie frames starts at the impact. The Hits selector turns the single hit into a series: a fixed number of hits at a fixed interval, all at the same speed or rising or falling between the softest and the loudest speed; the Excitation panel reports the schedule and how much ring the recorded duration leaves after the last hit. With the push excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles (its kinetic energy drops below the threshold or stops falling), and released by ramping the finger force to zero over the release time; the recording starts at the release. The view shows the shape of the sheet as a color map (2D) or a lit height field (3D; drag to orbit, wheel to zoom): by default the total shape, that is the initial shape plus the vibration, and with Show initial shape unticked only the deviation from the initial shape, which makes a small vibration on a bowed sheet visible. The spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.
2121
2222 - The sheet is a rectangle by default; the Outline selector makes it the ellipse inscribed in the given width and height. Width and height are entered in inches by default, or in millimeters; everything else is metric. Internally all values are SI.
23-- The sheet has a stress-free 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.
23+- The sheet has a stress-free initial shape, built from a list of features in the Initial shape panel: parabolic bows across the width or the height (with an optional apex offset), a twist, Gaussian dents, creases (a fold along a line), and a seeded smooth random field. The default is a single 3 mm bow across the width; Clear gives the ideal flat plate. Randomize draws a new combination whose heights are of the order of the scale field, so that the effect of the shape on the sound can be surveyed; the drawn values appear in the panel and can be edited. The panel reports the peak-to-peak height and the largest slope of the shape and warns when the slope leaves the shallow-shell range. Before Start, and after any change of a parameter that takes effect at the next Start, the view shows the initial shape alone; the movie slider or Live returns to the last run. A negative push amount pushes against a bow across the width.
24+- Copy link puts a URL with every parameter that differs from its default into the clipboard and the address bar, so a configuration (a random shape in particular) can be recovered or shared.
2425 - The Excitation panel shows the parameters of the selected excitation type only. Both excitations are self-equilibrated by inertia relief (the sheet gains no net momentum, as if the thread took it up); the periodic projection of the rigid-body modes during the ring is a Numerics option that is off by default and only removes rounding drift.
2526 - Sheet, initial shape, excitation, and numerics parameters take effect at the next Start. Losses and the listener position can be changed while a simulation runs.
2627 - Pause halts the run temporarily; Stop ends it, keeps the audio and movie, and unlocks the parameter fields. To extend a finished run, raise the duration and press Continue; the state is kept on the GPU.
2728 - 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.
2829 - After the simulation completes, drag the time slider to explore the movie; Play starts the audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.
29-- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?excitation=push&shape=ellipse` or `?hitSequence=crescendo&hitCount=6&hitInterval=0.5&duration=6`.
30+- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?excitation=push&shape=ellipse` or `?hitSequence=crescendo&hitCount=6&hitInterval=0.5&duration=6`. The initial shape is the `features` parameter, a list of `type:values` entries separated by semicolons with the values in SI in the order of the panel, for example `?features=bow:x,0.003,0;twist:0.001;dent:0.3,0.6,0.03,0.0015;noise:11,0.0005,0.05` or `?features=none`; the earlier `bowShape` and `bowAmplitude` parameters are still read.
3031
3132 The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interactive but lowers the bending partials by 3 percent near 1 kHz and by 10 percent near 3.6 kHz through grid dispersion; the Sheet panel reports these two frequencies for the chosen spacing. The stable time step is set by the in-plane wave speed, about 1.4 microseconds at 8 mm, so 4 s of sound take about 2.9 million steps. On an Intel Iris Xe integrated GPU this runs at about 57k steps per second, or roughly 50 s for 4 s of sound plus a few seconds for the push and hold; a discrete GPU is considerably faster. On grids of this size the cost of a step is the latency of one workgroup rather than the arithmetic, so a finer grid costs less than the node count suggests: halving the spacing to 4 mm costs about three times as much (about 2.5 minutes for 4 s on the same GPU) and moves the 3 percent point to about 4 kHz. The Numerics panel also offers an in-plane mass factor that lengthens the time step by its square root at the expense of lowering the in-plane modes of the sheet; see THEORY.md, Section 5.
3233
@@ -34,12 +35,13 @@ The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interacti
3435
3536 `npm test` runs two groups of tests with Node's built-in runner.
3637
37-- `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, forces as the exact gradient of the discrete energy (flat and bowed, rectangle and ellipse), dissipativity of the loss terms, and convergence of the free square and circular plate frequencies to the reference values of Gorman and Leissa.
38-- `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, and elliptical sheets (`?selftest=1` in the browser shows the same numbers), and runs a short push, a short strike, and a short series of strikes to completion. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first group.
38+- `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, forces as the exact gradient of the discrete energy (flat and bowed, rectangle and ellipse), dissipativity of the loss terms, and convergence of the free square and circular plate frequencies to the reference values of Gorman and Leissa. `tests/shape.test.ts` checks the initial-shape features: the text encoding, the removal of the mean plane, the geometry of each feature, and the random configurations.
39+- `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, elliptical, and randomly imperfect sheets (`?selftest=1` in the browser shows the same numbers), and runs a short push, a short strike, and a short series of strikes to completion. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first group.
3940
4041 ## Layout
4142
4243 - `src/params.ts`: parameters, grid and mask, material constants, geometry of the excitation, grid dispersion estimate.
44+- `src/shape.ts`: the initial-shape features, their evaluation on the grid, their text encoding, and the random configurations.
4345 - `src/sim/cpuRef.ts`: float64 reference of the spatial discretization; used for the time-step estimate and the tests.
4446 - `src/sim/shaders/`: WGSL kernels (time step, rigid-body projection, listener sum, display extraction).
4547 - `src/sim/simulation.ts`: WebGPU driver (buffers, batch encoding, phases, readback).
@@ -47,4 +49,4 @@ The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interacti
4749 - `src/render/`: 2D and 3D WebGPU rendering of the displacement field.
4850 - `src/audio/`: FFT, spectrogram, WAV encoder, playback, level mapping.
4951 - `src/ui/form.ts` and `src/main.ts`: parameter form and the rest of the interface.
50-- `tests/`: the two test files above.
52+- `tests/`: the three test files above.
THEORY.mdmodified+3−3View file
@@ -28,7 +28,7 @@ $$
2828
2929 with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. Bending acts on the deviation $w$ from the initial shape; the membrane forces act through the total slope. For $w_0 = 0$ these are the von Karman equations. The quadratic terms in the strains couple bending to stretching and produce the gong-like behavior: the pitch glide of a strongly excited plate and the transfer of energy from the low modes toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ are usually dropped in favor of an Airy stress function, since in-plane waves are much faster than bending waves. We keep them because they make the whole system explicit and local, which suits the GPU, at the cost of a time step set by the in-plane wave speed (Section 5).
3030
31-**Initial shape.** Real sheet metal is never flat; rolling and handling leave a bow of a few millimeters. A bowed sheet is a shallow shell, and its behavior differs qualitatively from a flat plate: the flat plate has a symmetric, purely hardening nonlinearity, whereas a shell has an asymmetric one, may snap between two curved states, and oscillates more slowly at large amplitude as it passes through the soft flat configuration. This is the mechanism behind oil-canning, and it appears to be what produces the slow wobble of a real sheet. Three parabolic shapes are offered: curved across the width, curved across the height, and a dish with the same curvature both ways. The sag is prescribed at the center of the bounding box; for the two cylindrical shapes it is measured relative to the straight edges, and for the dish relative to the corners, with the edge midpoints at half the sag. The default is a 3 mm bow across the width. A positive push moves the rim point in the direction of the bow and a negative push against it; the two cases sound different, as they do on a real sheet.
31+**Initial shape.** Real sheet metal is never flat; rolling and handling leave a bow of a few millimeters. A bowed sheet is a shallow shell, and its behavior differs qualitatively from a flat plate: the flat plate has a symmetric, purely hardening nonlinearity, whereas a shell has an asymmetric one, may snap between two curved states, and oscillates more slowly at large amplitude as it passes through the soft flat configuration. This is the mechanism behind oil-canning, and it appears to be what produces the slow wobble of a real sheet. The shape $w_0$ is the sum of a list of features, each of a few parameters: a parabolic bow across the width or the height, $a(1 - s^2)(1 + cs)$ with $s$ the normalized coordinate in $[-1, 1]$, whose apex offset $c$ moves the crest toward one edge; a twist $a\,XY$ with the height $a$ at the corners; a Gaussian dent of given center, $1/e$ radius, and height; a crease, which tilts the two sides of a line by half the fold angle each, $-\tan(\phi/2)(\sqrt{d^2 + s^2} - s)$ with $d$ the distance from the line and $s$ a rounding width; and a random field, a sum of 64 plane waves with Gaussian-distributed wavevectors, which has the covariance $\exp(-r^2/L^2)$ for a correlation length $L$, is scaled to a given RMS over the sheet, and is determined by its seed independently of the grid. The default is a single 3 mm bow across the width, and the interface can draw random combinations of the features with heights of a given scale. Only the slopes of $w_0$ enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation (Section 2), and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every shape. Note that bending acts on the deviation $w$ only, so a crease or a dent in $w_0$ carries no residual stress; this is the idealization of a permanent set, whereas a real crease is plastically worked and does carry some. The features should be broad compared with the grid spacing (a dent of 30 mm radius spans about eight cells at 8 mm) and their slopes small compared with one for the shallow-shell treatment to hold; the interface reports the largest slope and warns above 0.1. A positive push moves the rim point in the direction of a bow across the width and a negative push against it; the two cases sound different, as they do on a real sheet.
3232
3333 **Losses.** The two-parameter loss model is the standard one in plate sound synthesis. A bending mode of angular frequency $\omega$ and wavenumber $k$ decays at the rate $\sigma_0 + \sigma_1 k^2$, and for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster in proportion to frequency. The defaults, $\sigma_0 = 1.15\ \mathrm{s^{-1}}$ and $\sigma_1 = 5\times10^{-4}\ \mathrm{m^2/s}$, give 60 dB decay times of about 4.9 s at 100 Hz, 1.8 s at 1 kHz, and 0.5 s at 5 kHz. The in-plane motion carries an additional frequency-independent loss $\sigma_{\mathrm{in}}$ (default $20\ \mathrm{s^{-1}}$). The stretching modes of the sheet lie at several kilohertz, couple linearly to the transverse motion once the sheet is curved, and would otherwise ring for the whole low-frequency decay time; a material loss factor of about $10^{-3}$ at 5 kHz gives a rate of this order. None of these values are measured; they are chosen to sound plausible, and all three are exposed in the interface.
3434
@@ -60,7 +60,7 @@ The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $
6060
6161 **Accuracy.** Two errors limit the accuracy of the spatial scheme, and it is worth separating them. The first is the treatment of the edge. For the rectangle the scheme converges at second order: the first seven frequencies of a completely free square plate are low by about 3.5 percent at 10 cells per side, 0.8 percent at 20, and 0.3 percent at 30, compared with the reference values of Gorman. The staircase edge of the ellipse is only first-order accurate; the frequencies of a free circular plate are within about 2 percent of Leissa's values 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. The second error is grid dispersion, which affects every partial regardless of the outline. A plane bending wave of wavenumber $k$ along a grid axis has the discrete frequency ratio $(2 - 2\cos kh)/(kh)^2$, so a partial comes out low by 1 percent at $kh \approx 0.35$, by 3 percent at $kh \approx 0.60$, and by 10 percent at $kh \approx 1.12$. For the default sheet this places the 3 percent point near 1.1 kHz on the 8 mm grid, 1.9 kHz on a 6 mm grid, and 4.2 kHz on a 4 mm grid, and the 10 percent point near 3.6, 6.5, and 15 kHz respectively. The high partials of a gong sound are therefore audibly flattened at 8 mm, and the interface reports these two frequencies for the chosen spacing. Since the cost scales with the inverse cube of the spacing (Section 5), the default is a compromise in favor of interactive use, and a finer grid is a matter of patience rather than of the method.
6262
63-**Validation.** The 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 with and without a bow, that the loss terms only remove energy, that the free-plate frequencies converge as stated, and that the same holds on an elliptical mask. A further test runs one time step of the GPU kernel on a smooth random state of about 2.5 thicknesses amplitude, flat, bowed, and elliptical, and compares it with the float64 reference: the accelerations agree to about $10^{-6}$ relative, which is single-precision rounding.
63+**Validation.** The 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 with and without a bow, that the loss terms only remove energy, that the free-plate frequencies converge as stated, and that the same holds on an elliptical mask. A further test runs one time step of the GPU kernel on a smooth random state of about 2.5 thicknesses amplitude, on a flat, a bowed, an elliptical, and a randomly imperfect sheet (twist, dent, crease, and random field), and compares it with the float64 reference: the accelerations agree to about $10^{-6}$ relative, which is single-precision rounding.
6464
6565 ## 5. Time integration
6666
@@ -86,4 +86,4 @@ The plate is unbaffled in reality, so the sound from the two faces partly cancel
8686
8787 ## 7. Limitations
8888
89-The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a 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 a shallow-shell approximation, adequate for a bow of a few millimeters, 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. Grid dispersion lowers the high partials by the amounts given in Section 4, and the staircase edge of an elliptical sheet places the low frequencies within a few percent and splits modes that would be degenerate on a true ellipse. The suspension is idealized as inertia relief, so the swinging of the sheet on its thread is absent. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included, and 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.
89+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a 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 a shallow-shell approximation, adequate for a bow of a few millimeters, and the real shape of a given sheet is of course unknown; the offered features are idealizations, and the random field is a convenient generic imperfection rather than a measured one. The loss model is empirical and its parameters are not calibrated against a real sheet. Grid dispersion lowers the high partials by the amounts given in Section 4, and the staircase edge of an elliptical sheet places the low frequencies within a few percent and splits modes that would be degenerate on a true ellipse. The suspension is idealized as inertia relief, so the swinging of the sheet on its thread is absent. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included, and 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.
package.jsonmodified+1−1View file
@@ -9,7 +9,7 @@
99 "build": "tsc --noEmit && vite build",
1010 "preview": "vite preview",
1111 "test": "node --test tests/*.test.ts",
12- "test:cpu": "node --test tests/cpuRef.test.ts",
12+ "test:cpu": "node --test tests/cpuRef.test.ts tests/shape.test.ts",
1313 "test:gpu": "node --test tests/gpu.test.ts"
1414 },
1515 "devDependencies": {
src/main.tsmodified+116−17View file
@@ -1,7 +1,8 @@
11 // Application: parameter form, run loop, live and movie views, spectrogram,
22 // playback, WAV download, and the test hooks on window.steelsim.
33
4-import { defaultParams, type PhysParams } from "./params.ts";
4+import { defaultParams, type Grid, makeGrid, nodeWeights, type PhysParams } from "./params.ts";
5+import { encodeFeatures, initialShape } from "./shape.ts";
56 import { Simulation, type SimStatus } from "./sim/simulation.ts";
67 import { Renderer, type ViewMode } from "./render/renderer.ts";
78 import { Spectrogram } from "./audio/spectrogram.ts";
@@ -26,6 +27,9 @@ let live = true;
2627 let scrubFrame = 0;
2728 let movieBuf: GPUBuffer | null = null;
2829 let movieBufFrame = -1;
30+/** Show the stress-free initial shape of the current parameters instead of the simulation. */
31+let previewMode = true;
32+let preview: { grid: Grid; w0: GPUBuffer; zero: GPUBuffer; nw: GPUBuffer; wmax: number; stale: boolean } | null = null;
2933 const spec = new Spectrogram({ windowSize: 2048, hop: 512, fmax: 5000, dbFloor: 80, logFrequency: false });
3034 const player = new Player();
3135 let lastStatsTime = 0, lastStatsSteps = 0, stepsPerSec = 0;
@@ -37,6 +41,7 @@ const btnStart = mkButton("Start", "primary");
3741 const btnPause = mkButton("Pause");
3842 const btnStop = mkButton("Stop");
3943 const btnContinue = mkButton("Continue");
44+const btnLink = mkButton("Copy link");
4045 const durInput = document.createElement("input");
4146
4247 function mkButton(text: string, cls = ""): HTMLButtonElement {
@@ -71,8 +76,10 @@ function buildControls(): void {
7176 u.className = "unit";
7277 u.textContent = "s";
7378 durRow.append(l, durInput, u);
74- c.append(durRow, btnStart, btnPause, btnStop, btnContinue);
79+ c.append(durRow, btnStart, btnPause, btnStop, btnContinue, btnLink);
7580 btnStart.addEventListener("click", () => void start());
81+ btnLink.title = "Copy a link that opens the simulator with the current parameters";
82+ btnLink.addEventListener("click", () => void copyLink());
7683 btnPause.addEventListener("click", () => {
7784 paused = !paused;
7885 btnPause.textContent = paused ? "Resume" : "Pause";
@@ -104,6 +111,29 @@ function updateButtons(): void {
104111 btnContinue.disabled = !sim || running || !sim.canContinue();
105112 }
106113
114+/** The URL of the current parameters (those that differ from the defaults, in SI). */
115+function currentLink(): string {
116+ const q = form.toQuery();
117+ return `${location.origin}${location.pathname}${q ? "?" + q : ""}`;
118+}
119+
120+/** Copies the current link to the clipboard and puts it in the address bar. */
121+async function copyLink(): Promise<void> {
122+ const url = currentLink();
123+ history.replaceState(null, "", url);
124+ let ok = false;
125+ try {
126+ await navigator.clipboard.writeText(url);
127+ ok = true;
128+ } catch {
129+ window.prompt("Copy this link:", url);
130+ }
131+ if (ok) {
132+ btnLink.textContent = "Copied";
133+ setTimeout(() => (btnLink.textContent = "Copy link"), 1500);
134+ }
135+}
136+
107137 // ---------- View toolbars ----------
108138
109139 const viewMode = { value: "2d" as ViewMode };
@@ -150,11 +180,13 @@ function buildViewToolbar(): void {
150180 rangeInput = numInput(0, 0.5);
151181 shapeCheck = document.createElement("input");
152182 shapeCheck.type = "checkbox";
183+ shapeCheck.checked = true;
184+ shapeCheck.title = "Show the total shape (initial shape plus deviation); unticked, only the deviation from the initial shape";
153185 const hint = document.createElement("span");
154186 hint.className = "sub";
155187 hint.style.margin = "0";
156188 hint.textContent = "3D: drag to orbit, wheel to zoom";
157- tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), labeled("3D: add initial shape", shapeCheck), hint);
189+ tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), labeled("Show initial shape", shapeCheck), hint);
158190 }
159191
160192 function buildMovieToolbar(): void {
@@ -174,17 +206,19 @@ function buildMovieToolbar(): void {
174206 scrub.addEventListener("input", () => {
175207 live = false;
176208 liveBeforePlay = false;
209+ previewMode = false;
177210 scrubFrame = parseInt(scrub.value);
178211 if (player.playing) player.stop();
179212 updateLiveButton();
180213 });
181- btnLive.addEventListener("click", () => { live = true; liveBeforePlay = true; updateLiveButton(); });
214+ btnLive.addEventListener("click", () => { live = true; liveBeforePlay = true; previewMode = false; updateLiveButton(); });
182215 btnPlay.addEventListener("click", () => {
183216 if (player.playing) { player.stop(); endPlayback(); return; }
184217 if (!sim || sim.samples === 0) return;
185218 liveBeforePlay = live;
186219 const offset = live ? 0 : scrubFrame / sim.fps;
187220 live = false;
221+ previewMode = false;
188222 player.play(sim.recordedAudio(), sim.fs, offset, levelMode(), gainDb());
189223 btnPlay.textContent = "Stop";
190224 updateLiveButton();
@@ -275,7 +309,8 @@ async function start(): Promise<void> {
275309 btnPause.textContent = "Pause";
276310 params.duration = duration;
277311 sim = new Simulation(device, params);
278- renderer = new Renderer(device, $("view"), sim.grid, params, sim.w0, sim.nw);
312+ renderer?.setField(sim.grid, sim.w0, sim.nw, sim.w0Max);
313+ previewMode = false;
279314 movieBuf?.destroy();
280315 movieBuf = device.createBuffer({ size: Math.max(16, sim.grid.n * 4), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
281316 movieBufFrame = -1;
@@ -328,11 +363,32 @@ function frame(): void {
328363 }
329364
330365 function frameBody(): void {
331- if (!device || !sim || !renderer) return;
366+ if (!device || !renderer) return;
332367 renderer.mode = viewMode.value;
333368 renderer.exaggeration = parseFloat(exagInput.value) || 1;
334369 renderer.colorRange = (parseFloat(rangeInput.value) || 0) * 1e-3;
335370 renderer.showShape = shapeCheck.checked;
371+ const canvas = $<HTMLCanvasElement>("view");
372+ const cam = renderer.camera;
373+ const viewKey = [
374+ renderer.mode, renderer.exaggeration, renderer.colorRange, renderer.showShape,
375+ cam.yaw.toFixed(4), cam.pitch.toFixed(4), cam.dist.toFixed(4),
376+ canvas.clientWidth, canvas.clientHeight, window.devicePixelRatio,
377+ ];
378+ if (previewMode || !sim) {
379+ // The resting shape of the sheet that the next Start will use.
380+ ensurePreview(device);
381+ const renderKey = ["preview", ...viewKey, params.excitation, params.strikeFracX, params.strikeFracY, params.holeInset].join("|");
382+ if (renderKey !== lastRenderKey) {
383+ lastRenderKey = renderKey;
384+ renderPreview(device);
385+ }
386+ if (sim) {
387+ spec.draw($("spec"), sim.fs, Math.max(duration, sim.samples / sim.fs), null);
388+ updateInfo(sim.status());
389+ }
390+ return;
391+ }
336392 const pos = player.position();
337393 if (pos !== null) {
338394 scrubFrame = Math.min(sim.frames.length - 1, Math.floor(pos * sim.fps));
@@ -340,14 +396,7 @@ function frameBody(): void {
340396 }
341397 const st = sim.status();
342398 // Render only when something visible changed, so an idle or paused page does no work.
343- const canvas = $<HTMLCanvasElement>("view");
344- const cam = renderer.camera;
345- const renderKey = [
346- renderer.mode, renderer.exaggeration, renderer.colorRange, renderer.showShape,
347- cam.yaw.toFixed(4), cam.pitch.toFixed(4), cam.dist.toFixed(4),
348- live, scrubFrame, sim.frames.length, st.totalSteps, st.phase,
349- canvas.clientWidth, canvas.clientHeight, window.devicePixelRatio,
350- ].join("|");
399+ const renderKey = [...viewKey, live, scrubFrame, sim.frames.length, st.totalSteps, st.phase].join("|");
351400 if (renderKey !== lastRenderKey) {
352401 lastRenderKey = renderKey;
353402 renderView(device, st);
@@ -356,8 +405,55 @@ function frameBody(): void {
356405 updateInfo(st);
357406 }
358407
408+/** Builds the preview buffers from the current parameters when they are missing or stale. */
409+function ensurePreview(dev: GPUDevice): void {
410+ if (preview && !preview.stale) return;
411+ if (preview) for (const b of [preview.w0, preview.zero, preview.nw]) b.destroy();
412+ const grid = makeGrid(params);
413+ const shape = initialShape(params, grid);
414+ let wmax = 0;
415+ for (let k = 0; k < shape.length; k++) wmax = Math.max(wmax, Math.abs(shape[k]));
416+ const mk = (data: Float32Array<ArrayBuffer>) => {
417+ const b = dev.createBuffer({ size: Math.max(16, data.byteLength), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
418+ dev.queue.writeBuffer(b, 0, data);
419+ return b;
420+ };
421+ preview = { grid, w0: mk(shape), zero: mk(new Float32Array(grid.n)), nw: mk(nodeWeights(grid)), wmax, stale: false };
422+ renderer?.setField(grid, preview.zero, preview.nw);
423+ lastRenderKey = "";
424+}
425+
426+/** Marks the preview stale (parameters changed) and shows it; the next frame rebuilds it. */
427+function showPreview(): void {
428+ if (preview) preview.stale = true;
429+ previewMode = true;
430+ if (!sim) updateIdleInfo();
431+}
432+
433+function renderPreview(dev: GPUDevice): void {
434+ if (!renderer || !preview) return;
435+ // The renderer's "initial shape" buffer is zero here, and the shape itself is the displayed field.
436+ if (renderer.fieldGrid !== preview.grid) renderer.setField(preview.grid, preview.zero, preview.nw);
437+ const enc = dev.createCommandEncoder();
438+ renderer.render(enc, preview.w0, preview.wmax);
439+ timeLabel.textContent = "initial shape";
440+ dev.pushErrorScope("validation");
441+ dev.queue.submit([enc.finish()]);
442+ void dev.popErrorScope().then((e) => { if (e && !frameErrorShown) { frameErrorShown = true; console.error("frame validation error:", e.message); } });
443+}
444+
445+function updateIdleInfo(): void {
446+ const { grid, inside } = form.gridInfo();
447+ const html = `<div class="row"><span>Grid</span><span>${grid.nx} × ${grid.ny} nodes${grid.mask ? ` (${inside} in the sheet)` : ""}</span></div><p class="sub">The view shows the stress-free initial shape. Press Start to run.</p>`;
448+ if (html !== lastInfoHtml) {
449+ lastInfoHtml = html;
450+ $("info").innerHTML = html;
451+ }
452+}
453+
359454 function renderView(dev: GPUDevice, st: SimStatus): void {
360455 if (!sim || !renderer) return;
456+ if (renderer.fieldGrid !== sim.grid) renderer.setField(sim.grid, sim.w0, sim.nw, sim.w0Max);
361457 const enc = dev.createCommandEncoder();
362458 if (live || sim.frames.length === 0) {
363459 sim.encodeDisplay(enc);
@@ -435,7 +531,8 @@ async function selfTest(): Promise<void> {
435531 async function init(): Promise<void> {
436532 form = new ParamForm($("params"), params);
437533 form.onChange = (key, isLive) => {
438- if (!isLive || !sim) return;
534+ if (!isLive) { showPreview(); return; }
535+ if (!sim) return;
439536 if (key === "sigma0" || key === "sigma1" || key === "sigmaIn") sim.setDamping(params.sigma0, params.sigma1, params.sigmaIn);
440537 if (key.startsWith("listener")) sim.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation);
441538 };
@@ -462,8 +559,8 @@ async function init(): Promise<void> {
462559 device = await adapter.requestDevice();
463560 device.lost.then((info) => showNoGpu(`WebGPU device lost: ${info.message}`));
464561 device.addEventListener("uncapturederror", (ev) => console.error("WebGPU error:", (ev as GPUUncapturedErrorEvent).error.message));
465- const { grid } = form.gridInfo();
466- $("info").innerHTML = `<div class="row"><span>Grid</span><span>${grid.nx} × ${grid.ny} nodes</span></div><p class="sub">Press Start to run.</p>`;
562+ renderer = new Renderer(device, $("view"), params);
563+ updateIdleInfo();
467564 requestAnimationFrame(frame);
468565 if (q.get("selftest") === "1") await selfTest();
469566 else if (q.get("autostart") === "1") void start();
@@ -475,6 +572,8 @@ function exposeHooks(): void {
475572 status: () => (sim ? { ...sim.status(), dt: sim.dt, stepsPerSample: sim.stepsPerSample, nx: sim.grid.nx, ny: sim.grid.ny, stepsPerSec } : null),
476573 audio: () => (sim ? Array.from(sim.recordedAudio()) : null),
477574 frames: () => sim?.frames.length ?? 0,
575+ features: () => encodeFeatures(params.features),
576+ link: () => currentLink(),
478577 selftest: () => selfTestResult,
479578 noGpu: () => { const e = $("nogpu"); return e.hidden ? null : e.textContent; },
480579 };
src/params.tsmodified+5−30View file
@@ -2,8 +2,9 @@
22 // quantities (grid, material constants, geometry of the excitation).
33 // All values are SI.
44
5+import type { ShapeFeature } from "./shape.ts";
6+
57 export type PlateShape = "rectangle" | "ellipse";
6-export type BowShape = "none" | "width" | "height" | "dish";
78 export type Excitation = "push" | "strike";
89 export type HitSequence = "single" | "steady" | "crescendo" | "decrescendo";
910
@@ -17,9 +18,8 @@ export interface PhysParams {
1718 nu: number; // Poisson ratio
1819 rho: number; // density [kg/m^3]
1920 holeInset: number; // distance of the thread hole from the edges near the top-left corner [m]
20- // Stress-free initial shape
21- bowShape: BowShape;
22- bowAmplitude: number; // sag of the initial shape at the center [m]
21+ // Stress-free initial shape: the sum of these features (see shape.ts)
22+ features: ShapeFeature[];
2323 // Losses
2424 sigma0: number; // frequency-independent loss of all motion [1/s]
2525 sigma1: number; // loss proportional to the squared wavenumber, transverse motion [m^2/s]
@@ -66,8 +66,7 @@ export const defaultParams: PhysParams = {
6666 nu: 0.29,
6767 rho: 7900,
6868 holeInset: 0.012,
69- bowShape: "width",
70- bowAmplitude: 3e-3,
69+ features: [{ type: "bow", axis: "x", sag: 3e-3, offset: 0 }],
7170 sigma0: 1.15,
7271 sigma1: 5e-4,
7372 sigmaIn: 20,
@@ -324,30 +323,6 @@ export function pushPoint(p: PhysParams): [number, number] {
324323 return [p.Lx / 2, -p.Ly / 2];
325324 }
326325
327-/**
328- * Stress-free initial shape w0 at every node. All three shapes are parabolic
329- * with the given sag at the center of the bounding box: "width" and "height"
330- * are cylindrical with zero sag along the two opposite edges, and "dish" has
331- * the same curvature in both directions, so its sag is measured relative to
332- * the corners (the edge midpoints sit at half the sag).
333- */
334-export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
335- const w0 = new Float32Array(g.n);
336- if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;
337- for (let j = 0; j < g.ny; j++) {
338- for (let i = 0; i < g.nx; i++) {
339- const X = (2 * nodeX(g, i)) / p.Lx; // -1..1
340- const Y = (2 * nodeY(g, j)) / p.Ly;
341- let f = 0;
342- if (p.bowShape === "width") f = 1 - X * X;
343- else if (p.bowShape === "height") f = 1 - Y * Y;
344- else f = 1 - 0.5 * (X * X + Y * Y);
345- w0[j * g.nx + i] = p.bowAmplitude * f;
346- }
347- }
348- return w0;
349-}
350-
351326 /**
352327 * Nodes of the drumstick contact patch: the sheet node nearest the strike
353328 * point, plus every sheet node within the contact radius. The contact force
src/render/renderer.tsmodified+38−16View file
@@ -12,8 +12,8 @@ export class Renderer {
1212 exaggeration = 4;
1313 /** Color range |w| [m]; 0 means automatic (from the current maximum). */
1414 colorRange = 0;
15- /** Add the stress-free initial shape to the 3D height field. */
16- showShape = false;
15+ /** Show the total shape (initial shape plus deviation) rather than the deviation alone, in both views. */
16+ showShape = true;
1717 camera = { yaw: 0.6, pitch: 0.5, dist: 1.0 };
1818
1919 private device: GPUDevice;
@@ -24,22 +24,26 @@ export class Renderer {
2424 private pipe2d: GPURenderPipeline;
2525 private pipe3d: GPURenderPipeline;
2626 private bgLayout: GPUBindGroupLayout;
27- private index: GPUBuffer;
28- private indexCount: number;
27+ private index: GPUBuffer | null = null;
28+ private indexCount = 0;
2929 private depth: GPUTexture | null = null;
30- private grid: Grid;
30+ private grid: Grid | null = null;
3131 private params: PhysParams;
3232 private bindGroups = new Map<GPUBuffer, GPUBindGroup>();
3333
34- private w0: GPUBuffer;
35- private nw: GPUBuffer;
34+ private w0: GPUBuffer | null = null;
35+ private w0Max = 0;
36+ private nw: GPUBuffer | null = null;
3637
37- constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams, w0: GPUBuffer, nw: GPUBuffer) {
38+ /** The grid of the field currently set, or null before setField. */
39+ get fieldGrid(): Grid | null {
40+ return this.grid;
41+ }
42+
43+ /** The renderer draws nothing until setField gives it a grid; `params` is read at every render for the markers. */
44+ constructor(device: GPUDevice, canvas: HTMLCanvasElement, params: PhysParams) {
3845 this.device = device;
39- this.w0 = w0;
40- this.nw = nw;
4146 this.canvas = canvas;
42- this.grid = grid;
4347 this.params = params;
4448 this.ctx = canvas.getContext("webgpu") as GPUCanvasContext;
4549 this.format = navigator.gpu.getPreferredCanvasFormat();
@@ -68,6 +72,22 @@ export class Renderer {
6872 primitive: { topology: "triangle-list", cullMode: "none" },
6973 depthStencil: { format: "depth24plus", depthWriteEnabled: true, depthCompare: "less" },
7074 });
75+ this.installControls();
76+ }
77+
78+ /**
79+ * Sets the grid and the per-node buffers the views read: the stress-free
80+ * shape (added to the displayed field when showShape is on; w0Max is its
81+ * largest magnitude, for the automatic color range) and the node weights
82+ * (zero outside the sheet). Called at every Start and for the preview of
83+ * the initial shape; the camera and the canvas controls persist.
84+ */
85+ setField(grid: Grid, w0: GPUBuffer, nw: GPUBuffer, w0Max = 0): void {
86+ this.grid = grid;
87+ this.w0 = w0;
88+ this.w0Max = w0Max;
89+ this.nw = nw;
90+ this.bindGroups.clear();
7191 // Triangle indices of the height field (cells of the sheet only).
7292 const { nx, ny } = grid;
7393 const idx = new Uint32Array((nx - 1) * (ny - 1) * 6);
@@ -81,9 +101,9 @@ export class Renderer {
81101 }
82102 }
83103 this.indexCount = k;
84- this.index = device.createBuffer({ size: Math.max(16, idx.byteLength), usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
85- device.queue.writeBuffer(this.index, 0, idx);
86- this.installControls();
104+ this.index?.destroy();
105+ this.index = this.device.createBuffer({ size: Math.max(16, idx.byteLength), usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
106+ this.device.queue.writeBuffer(this.index, 0, idx);
87107 }
88108
89109 private installControls(): void {
@@ -119,11 +139,13 @@ export class Renderer {
119139 }
120140 }
121141
122- /** Renders the w field held in `wBuffer` (f32 per node). `wmax` drives the automatic color range. */
142+ /** Renders the w field held in `wBuffer` (f32 per node). `wmax`, the largest |w|, drives the automatic color range. */
123143 render(encoder: GPUCommandEncoder, wBuffer: GPUBuffer, wmax: number): void {
144+ if (!this.grid || !this.w0 || !this.nw || !this.index) return;
124145 this.resize();
125146 const aspect = this.canvas.width / this.canvas.height;
126- const range = this.colorRange > 0 ? this.colorRange : Math.max(wmax, 1e-6);
147+ const auto = wmax + (this.showShape ? this.w0Max : 0);
148+ const range = this.colorRange > 0 ? this.colorRange : Math.max(auto, 1e-6);
127149 const u = new Float32Array(36);
128150 const ui = new Uint32Array(u.buffer);
129151 const eye = [
src/render/shaders.wgslmodified+7−10View file
@@ -16,7 +16,7 @@ struct RU {
1616 lightX: f32,
1717 lightY: f32,
1818 lightZ: f32,
19- showShape: u32, // 3D: add the stress-free initial shape to the height field
19+ showShape: u32, // add the stress-free initial shape to the displayed field (both views)
2020 shape: u32, // 0: rectangle, 1: ellipse inscribed in the bounding box
2121 pad1: u32,
2222 pad2: u32,
@@ -36,23 +36,20 @@ fn colormap(t: f32) -> vec3f {
3636 return mix(mid, red, s);
3737 }
3838
39+// Displayed field: the deviation w, plus the initial shape when requested (the total shape).
3940 fn wAt(i: i32, j: i32) -> f32 {
4041 let ii = clamp(i, 0, i32(R.nx) - 1);
4142 let jj = clamp(j, 0, i32(R.ny) - 1);
42- return w[u32(jj) * R.nx + u32(ii)];
43+ var z = w[u32(jj) * R.nx + u32(ii)];
44+ if (R.showShape != 0u) { z += w0[u32(jj) * R.nx + u32(ii)]; }
45+ return z;
4346 }
4447 fn nwAt(i: i32, j: i32) -> f32 {
4548 if (i < 0 || i >= i32(R.nx) || j < 0 || j >= i32(R.ny)) { return 0.0; }
4649 return nw[u32(j) * R.nx + u32(i)];
4750 }
48-// Height used by the 3D view: w, plus the initial shape when requested.
49-fn zAt(i: i32, j: i32) -> f32 {
50- let ii = clamp(i, 0, i32(R.nx) - 1);
51- let jj = clamp(j, 0, i32(R.ny) - 1);
52- var z = w[u32(jj) * R.nx + u32(ii)];
53- if (R.showShape != 0u) { z += w0[u32(jj) * R.nx + u32(ii)]; }
54- return z;
55-}
51+// Height used by the 3D view: the same field as the color.
52+fn zAt(i: i32, j: i32) -> f32 { return wAt(i, j); }
5653 // Slope along one axis, one-sided where a neighbor lies outside the sheet.
5754 fn slope(i: i32, j: i32, di: i32, dj: i32, h: f32) -> f32 {
5855 let hasP = nwAt(i + di, j + dj) > 0.0;
src/selftest.tsmodified+7−5View file
@@ -2,7 +2,8 @@
22 // strongly nonlinear random state is compared with the float64 reference.
33 // Runs in the browser (?selftest=1) and is driven by tests/gpu.test.ts.
44
5-import { defaultParams, initialShape, makeMaterial, nodeWeight, type PhysParams } from "./params.ts";
5+import { defaultParams, makeMaterial, nodeWeight, type PhysParams } from "./params.ts";
6+import { initialShape, parseFeatures } from "./shape.ts";
67 import { accelerations, addDamping, fields } from "./sim/cpuRef.ts";
78 import { Simulation } from "./sim/simulation.ts";
89
@@ -20,10 +21,11 @@ export const ACC_TOLERANCE = 1e-3;
2021 export const POSITION_TOLERANCE = 1e-2;
2122
2223 const cases: { name: string; overrides: Partial<PhysParams> }[] = [
23- { name: "rectangle, flat", overrides: { shape: "rectangle", bowShape: "none" } },
24- { name: "rectangle, bowed across the width", overrides: { shape: "rectangle", bowShape: "width" } },
25- { name: "ellipse, dish", overrides: { shape: "ellipse", bowShape: "dish" } },
26- { name: "rectangle, bowed, in-plane mass x4", overrides: { shape: "rectangle", bowShape: "width", inPlaneMass: 4 } },
24+ { name: "rectangle, flat", overrides: { shape: "rectangle", features: [] } },
25+ { name: "rectangle, bowed across the width", overrides: { shape: "rectangle" } },
26+ { name: "ellipse, dish", overrides: { shape: "ellipse", features: parseFeatures("bow:x,0.0015,0;bow:y,0.0015,0")! } },
27+ { name: "rectangle, twist, dent, crease, and random field", overrides: { shape: "rectangle", features: parseFeatures("bow:x,0.002,0.5;twist:0.001;dent:0.3,0.6,0.03,0.0015;crease:0.6,0.4,25,-1.5,0.015;noise:11,0.0005,0.05")! } },
28+ { name: "rectangle, bowed, in-plane mass x4", overrides: { shape: "rectangle", inPlaneMass: 4 } },
2729 ];
2830
2931 /** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
src/shape.tsadded+331−0View file
@@ -0,0 +1,331 @@
1+// Stress-free initial shape of the sheet. The shape is a list of features
2+// (bows, a twist, dents, creases, a random field) summed at the nodes, with
3+// the mean plane removed. Only the slopes of the shape enter the physics, so
4+// the mean is immaterial; the tilt is not, because the von Karman strains are
5+// not invariant under rotation, and a tilted stress-free shape would add a
6+// spurious membrane coupling. Every shape is therefore made plane-free here.
7+
8+import { type Grid, type PhysParams, insideNode, nodeWeight, nodeX, nodeY } from "./params.ts";
9+
10+export type FeatureType = "bow" | "twist" | "dent" | "crease" | "noise";
11+
12+/** Parabolic bow across the width (axis x) or the height (axis y); the offset in -1..1 moves the crest toward one edge. */
13+export interface BowFeature { type: "bow"; axis: "x" | "y"; sag: number; offset: number }
14+/** Saddle x*y with the given height at the corners (alternating in sign). */
15+export interface TwistFeature { type: "twist"; height: number }
16+/** Gaussian bump at (x, y) given as fractions of the width and the height, with the 1/e radius and the height. */
17+export interface DentFeature { type: "dent"; x: number; y: number; radius: number; height: number }
18+/** Fold along the line through (x, y) at `angle` degrees from the x axis, kinked by `fold` degrees (positive: ridge), rounded over `width`. */
19+export interface CreaseFeature { type: "crease"; x: number; y: number; angle: number; fold: number; width: number }
20+/** Seeded random field with Gaussian correlation of length `length` and the given RMS over the sheet. */
21+export interface NoiseFeature { type: "noise"; seed: number; rms: number; length: number }
22+export type ShapeFeature = BowFeature | TwistFeature | DentFeature | CreaseFeature | NoiseFeature;
23+
24+export interface FeatureField {
25+ key: string;
26+ label: string;
27+ unit: string;
28+ scale: number; // form value = SI value / scale
29+ step: number;
30+ options?: [string, string][]; // a choice rather than a number
31+}
32+
33+export interface FeatureSpec {
34+ label: string;
35+ description: string;
36+ fields: FeatureField[];
37+}
38+
39+/** Presentation and encoding order of the fields of each feature type. */
40+export const FEATURE_SPECS: Record<FeatureType, FeatureSpec> = {
41+ bow: {
42+ label: "Bow",
43+ description: "Parabolic bow with zero sag along two opposite edges. The apex offset (-1 to 1) moves the crest toward one edge.",
44+ fields: [
45+ { key: "axis", label: "Across", unit: "", scale: 1, step: 1, options: [["x", "the width"], ["y", "the height"]] },
46+ { key: "sag", label: "Sag", unit: "mm", scale: 1e-3, step: 0.5 },
47+ { key: "offset", label: "Apex offset", unit: "", scale: 1, step: 0.1 },
48+ ],
49+ },
50+ twist: {
51+ label: "Twist",
52+ description: "Saddle shape: two opposite corners up by the height, the other two down.",
53+ fields: [{ key: "height", label: "Height at the corners", unit: "mm", scale: 1e-3, step: 0.5 }],
54+ },
55+ dent: {
56+ label: "Dent",
57+ description: "Gaussian bump; the radius is where the height has fallen to 1/e. A negative height is a dent toward the back.",
58+ fields: [
59+ { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
60+ { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
61+ { key: "radius", label: "Radius", unit: "mm", scale: 1e-3, step: 5 },
62+ { key: "height", label: "Height", unit: "mm", scale: 1e-3, step: 0.5 },
63+ ],
64+ },
65+ crease: {
66+ label: "Crease",
67+ description: "Fold along a line: the two sides tilt by half the fold angle each. A positive fold is a ridge, a negative one a valley; the width rounds the kink.",
68+ fields: [
69+ { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
70+ { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
71+ { key: "angle", label: "Line angle from x axis", unit: "deg", scale: 1, step: 15 },
72+ { key: "fold", label: "Fold angle", unit: "deg", scale: 1, step: 0.5 },
73+ { key: "width", label: "Rounding width", unit: "mm", scale: 1e-3, step: 5 },
74+ ],
75+ },
76+ noise: {
77+ label: "Random field",
78+ description: "Smooth random undulation with Gaussian correlation; the same seed gives the same field on any grid.",
79+ fields: [
80+ { key: "seed", label: "Seed", unit: "", scale: 1, step: 1 },
81+ { key: "rms", label: "RMS amplitude", unit: "mm", scale: 1e-3, step: 0.1 },
82+ { key: "length", label: "Correlation length", unit: "mm", scale: 1e-3, step: 10 },
83+ ],
84+ },
85+};
86+
87+export const FEATURE_TYPES = Object.keys(FEATURE_SPECS) as FeatureType[];
88+
89+/** A feature of the given type with sensible default values. */
90+export function newFeature(type: FeatureType): ShapeFeature {
91+ switch (type) {
92+ case "bow": return { type, axis: "x", sag: 3e-3, offset: 0 };
93+ case "twist": return { type, height: 1e-3 };
94+ case "dent": return { type, x: 0.5, y: 0.5, radius: 0.03, height: 1e-3 };
95+ case "crease": return { type, x: 0.5, y: 0.5, angle: 30, fold: 1, width: 0.015 };
96+ case "noise": return { type, seed: 1, rms: 0.5e-3, length: 0.06 };
97+ }
98+}
99+
100+export function cloneFeatures(features: ShapeFeature[]): ShapeFeature[] {
101+ return features.map((f) => ({ ...f }));
102+}
103+
104+/** The shapes of the earlier interface: a parabolic bow across the width, the height, or both (dish). */
105+export function legacyFeatures(bowShape: string, amplitude: number): ShapeFeature[] {
106+ if (bowShape === "none" || amplitude === 0) return [];
107+ if (bowShape === "height") return [{ type: "bow", axis: "y", sag: amplitude, offset: 0 }];
108+ if (bowShape === "dish") return [{ type: "bow", axis: "x", sag: amplitude / 2, offset: 0 }, { type: "bow", axis: "y", sag: amplitude / 2, offset: 0 }];
109+ return [{ type: "bow", axis: "x", sag: amplitude, offset: 0 }];
110+}
111+
112+// ---------- Text encoding (URL) ----------
113+
114+function fmt(v: number): string {
115+ return v === 0 ? "0" : String(+v.toPrecision(8));
116+}
117+
118+/** Compact text form: `type:v1,v2,...;type:...` with the values in SI, or `none` for a flat sheet. */
119+export function encodeFeatures(features: ShapeFeature[]): string {
120+ if (features.length === 0) return "none";
121+ return features
122+ .map((f) => `${f.type}:${FEATURE_SPECS[f.type].fields.map((fld) => {
123+ const v = (f as unknown as Record<string, number | string>)[fld.key];
124+ return typeof v === "number" ? fmt(v) : String(v);
125+ }).join(",")}`)
126+ .join(";");
127+}
128+
129+/** Inverse of encodeFeatures. Returns null when nothing in the text can be read; entries that cannot be read are skipped. */
130+export function parseFeatures(text: string): ShapeFeature[] | null {
131+ const s = text.trim();
132+ if (s === "" || s === "none" || s === "flat") return [];
133+ const out: ShapeFeature[] = [];
134+ let any = false;
135+ for (const item of s.split(";")) {
136+ const colon = item.indexOf(":");
137+ if (colon < 0) continue;
138+ const type = item.slice(0, colon).trim() as FeatureType;
139+ const spec = FEATURE_SPECS[type];
140+ if (!spec) continue;
141+ const vals = item.slice(colon + 1).split(",").map((v) => v.trim());
142+ if (vals.length !== spec.fields.length) continue;
143+ const f = newFeature(type) as unknown as Record<string, number | string>;
144+ let ok = true;
145+ spec.fields.forEach((fld, i) => {
146+ if (fld.options) {
147+ if (!fld.options.some(([v]) => v === vals[i])) ok = false;
148+ else f[fld.key] = vals[i];
149+ } else {
150+ const num = parseFloat(vals[i]);
151+ if (!isFinite(num)) ok = false;
152+ else f[fld.key] = num;
153+ }
154+ });
155+ if (!ok) continue;
156+ out.push(f as unknown as ShapeFeature);
157+ any = true;
158+ }
159+ return any ? out : null;
160+}
161+
162+// ---------- Evaluation ----------
163+
164+/** Deterministic uniform generator on [0, 1) (mulberry32). */
165+function mulberry32(seed: number): () => number {
166+ let a = Math.floor(seed) >>> 0;
167+ return () => {
168+ a = (a + 0x6d2b79f5) >>> 0;
169+ let t = a;
170+ t = Math.imul(t ^ (t >>> 15), t | 1);
171+ t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
172+ return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
173+ };
174+}
175+
176+const NOISE_TERMS = 64;
177+
178+/**
179+ * Adds the field of one feature to f (float64 per node, in meters). Nodes
180+ * outside the sheet are filled too and zeroed at the end.
181+ */
182+function addFeature(f: Float64Array, feat: ShapeFeature, p: PhysParams, g: Grid): void {
183+ const { nx, ny } = g;
184+ if (feat.type === "noise") {
185+ // Sum of random plane waves whose wavevectors are Gaussian with standard
186+ // deviation sqrt(2)/L per component, which gives the covariance exp(-r^2/L^2).
187+ const rnd = mulberry32(feat.seed);
188+ const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());
189+ const L = Math.max(feat.length, 1e-4);
190+ const terms = Array.from({ length: NOISE_TERMS }, () => ({ kx: (gauss() * Math.SQRT2) / L, ky: (gauss() * Math.SQRT2) / L, ph: 2 * Math.PI * rnd() }));
191+ const raw = new Float64Array(g.n);
192+ let sw = 0, s1 = 0, s2 = 0;
193+ for (let j = 0; j < ny; j++) {
194+ for (let i = 0; i < nx; i++) {
195+ const x = nodeX(g, i), y = nodeY(g, j);
196+ let v = 0;
197+ for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);
198+ raw[j * nx + i] = v;
199+ const wgt = nodeWeight(g, i, j);
200+ sw += wgt; s1 += wgt * v; s2 += wgt * v * v;
201+ }
202+ }
203+ const mean = s1 / sw, rms = Math.sqrt(Math.max(0, s2 / sw - mean * mean));
204+ const scale = rms > 0 ? feat.rms / rms : 0;
205+ for (let k = 0; k < g.n; k++) f[k] += scale * (raw[k] - mean);
206+ return;
207+ }
208+ for (let j = 0; j < ny; j++) {
209+ for (let i = 0; i < nx; i++) {
210+ const x = nodeX(g, i), y = nodeY(g, j);
211+ const X = (2 * x) / p.Lx, Y = (2 * y) / p.Ly; // -1..1 over the bounding box
212+ let v = 0;
213+ switch (feat.type) {
214+ case "bow": {
215+ const s = feat.axis === "x" ? X : Y;
216+ v = feat.sag * (1 - s * s) * (1 + feat.offset * s);
217+ break;
218+ }
219+ case "twist":
220+ v = feat.height * X * Y;
221+ break;
222+ case "dent": {
223+ const dx = x - (feat.x - 0.5) * p.Lx, dy = y - (feat.y - 0.5) * p.Ly;
224+ const r = Math.max(feat.radius, 1e-4);
225+ v = feat.height * Math.exp(-(dx * dx + dy * dy) / (r * r));
226+ break;
227+ }
228+ case "crease": {
229+ const th = (feat.angle * Math.PI) / 180;
230+ const d = -(x - (feat.x - 0.5) * p.Lx) * Math.sin(th) + (y - (feat.y - 0.5) * p.Ly) * Math.cos(th);
231+ const s = Math.max(feat.width, 0);
232+ v = -Math.tan((feat.fold * Math.PI) / 360) * (Math.sqrt(d * d + s * s) - s);
233+ break;
234+ }
235+ }
236+ f[j * nx + i] += v;
237+ }
238+ }
239+}
240+
241+/** Subtracts the mass-weighted least-squares plane a + b x + c y from f over the sheet nodes. */
242+export function removePlane(f: Float64Array, g: Grid): void {
243+ let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
244+ for (let j = 0; j < g.ny; j++) {
245+ for (let i = 0; i < g.nx; i++) {
246+ const wgt = nodeWeight(g, i, j);
247+ if (wgt === 0) continue;
248+ const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
249+ sw += wgt; sx += wgt * x; sy += wgt * y;
250+ sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
251+ sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
252+ }
253+ }
254+ if (sw === 0) return;
255+ // Normal equations of the 3-parameter fit, solved by Cramer's rule.
256+ const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
257+ if (Math.abs(det) < 1e-300) return;
258+ const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
259+ const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
260+ const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
261+ 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);
262+}
263+
264+/** Stress-free initial shape w0 at every node: the sum of the features with the mean plane removed; zero outside the sheet. */
265+export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
266+ const w0 = new Float32Array(g.n);
267+ if (p.features.length === 0) return w0;
268+ const f = new Float64Array(g.n);
269+ for (const feat of p.features) addFeature(f, feat, p, g);
270+ removePlane(f, g);
271+ 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] = f[j * g.nx + i];
272+ return w0;
273+}
274+
275+/** Peak-to-peak height and the largest slope of a shape over the sheet nodes. */
276+export function shapeStats(w0: ArrayLike<number>, g: Grid): { peakToPeak: number; maxSlope: number } {
277+ let lo = Infinity, hi = -Infinity, maxSlope = 0;
278+ const at = (i: number, j: number) => w0[j * g.nx + i];
279+ const slope = (i: number, j: number, di: number, dj: number, h: number) => {
280+ const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
281+ if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * h);
282+ if (hasP) return (at(i + di, j + dj) - at(i, j)) / h;
283+ if (hasM) return (at(i, j) - at(i - di, j - dj)) / h;
284+ return 0;
285+ };
286+ for (let j = 0; j < g.ny; j++) {
287+ for (let i = 0; i < g.nx; i++) {
288+ if (!insideNode(g, i, j)) continue;
289+ const v = at(i, j);
290+ lo = Math.min(lo, v); hi = Math.max(hi, v);
291+ maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
292+ }
293+ }
294+ return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
295+}
296+
297+// ---------- Random configurations ----------
298+
299+/**
300+ * A random feature list whose heights are of the order of `scale` (m): always
301+ * a bow across one axis, sometimes one across the other, a twist, up to two
302+ * dents, a crease, and a random field. `rnd` is a uniform generator on [0, 1).
303+ */
304+export function randomFeatures(p: PhysParams, scale: number, rnd: () => number = Math.random): ShapeFeature[] {
305+ const U = (a: number, b: number) => a + (b - a) * rnd();
306+ const sign = () => (rnd() < 0.5 ? -1 : 1);
307+ const Lmin = Math.min(p.Lx, p.Ly);
308+ const out: ShapeFeature[] = [];
309+ const first: "x" | "y" = rnd() < 0.5 ? "x" : "y";
310+ out.push({ type: "bow", axis: first, sag: sign() * U(0.3, 1) * scale, offset: U(-0.6, 0.6) });
311+ if (rnd() < 0.5) out.push({ type: "bow", axis: first === "x" ? "y" : "x", sag: sign() * U(0.1, 0.5) * scale, offset: U(-0.6, 0.6) });
312+ if (rnd() < 0.5) out.push({ type: "twist", height: sign() * U(0.1, 0.5) * scale });
313+ const dents = rnd() < 0.3 ? 0 : rnd() < 0.6 ? 1 : 2;
314+ for (let k = 0; k < dents; k++) out.push({ type: "dent", x: U(0.15, 0.85), y: U(0.15, 0.85), radius: U(0.1, 0.2) * Lmin, height: sign() * U(0.2, 0.7) * scale });
315+ if (rnd() < 0.3) {
316+ // A rise of a fraction of the scale over a quarter of the smaller side.
317+ const slope = (U(0.2, 0.6) * scale) / (Lmin / 4);
318+ out.push({ type: "crease", x: U(0.3, 0.7), y: U(0.3, 0.7), angle: U(0, 180), fold: sign() * ((2 * Math.atan(slope) * 180) / Math.PI), width: 0.05 * Lmin });
319+ }
320+ if (rnd() < 0.7) out.push({ type: "noise", seed: 1 + Math.floor(rnd() * 99999), rms: U(0.05, 0.25) * scale, length: U(0.1, 0.35) * Lmin });
321+ // Round the values to what the form shows, so the list reads the same in the form and in a link.
322+ for (const f of out) {
323+ const rec = f as unknown as Record<string, number | string>;
324+ for (const fld of FEATURE_SPECS[f.type].fields) {
325+ if (fld.options || fld.key === "seed") continue;
326+ const v = rec[fld.key] as number;
327+ rec[fld.key] = +((v / fld.scale).toPrecision(4)) * fld.scale;
328+ }
329+ }
330+ return out;
331+}
src/sim/simulation.tsmodified+6−2View file
@@ -16,7 +16,6 @@ import {
1616 type Material,
1717 makeGrid,
1818 makeMaterial,
19- initialShape,
2019 strikePoint,
2120 pushPoint,
2221 contactPatch,
@@ -35,6 +34,7 @@ import {
3534 SETTLE_WINDOW,
3635 SPEED_OF_SOUND,
3736 } from "../params.ts";
37+import { initialShape } from "../shape.ts";
3838 import { maxAngularFrequency } from "./cpuRef.ts";
3939 import commonWgsl from "./shaders/common.wgsl?raw";
4040 import stepWgsl from "./shaders/step.wgsl?raw";
@@ -141,6 +141,8 @@ export class Simulation {
141141 private react: GPUBuffer;
142142 /** Stress-free initial shape (f32 per node), also used by the renderer. */
143143 readonly w0: GPUBuffer;
144+ /** Largest magnitude of the initial shape [m]. */
145+ readonly w0Max: number;
144146 /** Node weights (f32 per node, 0 outside the sheet), also used by the renderer. */
145147 readonly nw: GPUBuffer;
146148 private lw: GPUBuffer;
@@ -194,7 +196,9 @@ export class Simulation {
194196 this.hist = f32(this.histLen * n);
195197 this.react = f32(n + STICK_LEN);
196198 this.w0 = f32(n);
197- device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
199+ const shape = initialShape(params, g);
200+ device.queue.writeBuffer(this.w0, 0, shape);
201+ this.w0Max = shape.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
198202 this.nw = f32(n);
199203 device.queue.writeBuffer(this.nw, 0, nodeWeights(g));
200204 this.lw = f32(2 * n);
src/style.cssmodified+7−0View file
@@ -60,3 +60,10 @@ button.primary:hover { background: #f8c453; }
6060 .field select.unit-select { padding: 3px 2px; font-size: 11px; }
6161 .cacherow { display: flex; justify-content: space-between; align-items: center; gap: 6px; margin: 3px 0; font-size: 12px; color: var(--muted); }
6262 .cacherow span { overflow-wrap: anywhere; }
63+.feature { border: 1px solid var(--border); border-radius: 5px; padding: 3px 6px 4px; margin: 5px 0; }
64+.feature-head { display: flex; justify-content: space-between; align-items: center; margin-bottom: 2px; color: var(--muted); font-size: 12px; }
65+.feature-row { display: flex; align-items: center; gap: 6px; margin: 5px 0; }
66+.feature-row select { width: auto; flex: 1; }
67+.feature-row input[type=number] { width: 64px; }
68+.feature-row label { color: var(--muted); }
69+button.small { padding: 1px 6px; font-size: 11px; }
src/ui/form.tsmodified+229−12View file
@@ -1,14 +1,15 @@
1-// Parameter form: numeric fields in convenient units, a few selects, and the
2-// derived-quantity notes. Everything is stored in SI in `params`.
1+// Parameter form: numeric fields in convenient units, a few selects, the
2+// editable list of initial-shape features, and the derived-quantity notes.
3+// Everything is stored in SI in `params`.
34
45 import {
5- type BowShape,
66 type Excitation,
77 type Grid,
88 hitSchedule,
99 type PhysParams,
1010 type PlateShape,
1111 countInside,
12+ defaultParams,
1213 dispersionFrequency,
1314 makeGrid,
1415 makeMaterial,
@@ -16,6 +17,19 @@ import {
1617 CFL,
1718 } from "../params.ts";
1819 import { maxAngularFrequency } from "../sim/cpuRef.ts";
20+import {
21+ FEATURE_SPECS,
22+ FEATURE_TYPES,
23+ type FeatureType,
24+ type ShapeFeature,
25+ encodeFeatures,
26+ initialShape,
27+ legacyFeatures,
28+ newFeature,
29+ parseFeatures,
30+ randomFeatures,
31+ shapeStats,
32+} from "../shape.ts";
1933
2034 export type SizeUnit = "in" | "mm";
2135
@@ -41,7 +55,8 @@ interface Group {
4155 open: boolean;
4256 selects?: SelectSpec[]; // shown at the top of the group
4357 fields: FieldSpec[];
44- note?: string | "grid" | "damping" | "numerics" | "hits";
58+ features?: boolean; // the group holds the initial-shape feature list
59+ note?: string | "grid" | "damping" | "numerics" | "hits" | "shape";
4560 }
4661
4762 const push = (p: PhysParams) => p.excitation === "push";
@@ -49,6 +64,9 @@ const strike = (p: PhysParams) => p.excitation === "strike";
4964 const series = (p: PhysParams) => strike(p) && p.hitSequence !== "single";
5065 const graded = (p: PhysParams) => strike(p) && (p.hitSequence === "crescendo" || p.hitSequence === "decrescendo");
5166
67+/** Slope above which the shallow-shell treatment of the initial shape is doubtful. */
68+const SLOPE_LIMIT = 0.1;
69+
5270 const groups: Group[] = [
5371 {
5472 title: "Sheet",
@@ -68,9 +86,9 @@ const groups: Group[] = [
6886 {
6987 title: "Initial shape",
7088 open: true,
71- selects: [{ key: "bowShape", label: "Bow shape", options: [["width", "Curved across the width"], ["height", "Curved across the height"], ["dish", "Dish (curved both ways)"], ["none", "Flat"]] }],
72- fields: [{ key: "bowAmplitude", label: "Bow (center sag)", unit: "mm", scale: 1e-3, step: 0.5 }],
73- note: "Real sheet metal is never flat. A bow of a few millimeters makes the sheet a shallow shell with an asymmetric stiffness; a negative push amount pushes against the bow.",
89+ fields: [],
90+ features: true,
91+ note: "shape",
7492 },
7593 {
7694 title: "Excitation",
@@ -134,10 +152,9 @@ const groups: Group[] = [
134152 ];
135153
136154 const SIZE_SCALE: Record<SizeUnit, number> = { in: 0.0254, mm: 1e-3 };
137-const selectKeys: (keyof PhysParams)[] = ["shape", "bowShape", "excitation", "hitSequence", "rigidProjection"];
155+const selectKeys: (keyof PhysParams)[] = ["shape", "excitation", "hitSequence", "rigidProjection"];
138156 const selectValues: Record<string, string[]> = {
139157 shape: ["rectangle", "ellipse"],
140- bowShape: ["none", "width", "height", "dish"],
141158 excitation: ["push", "strike"],
142159 hitSequence: ["single", "steady", "crescendo", "decrescendo"],
143160 rigidProjection: ["false", "true"],
@@ -155,9 +172,16 @@ export function formatNumber(v: number): string {
155172 return String(+v.toPrecision(6));
156173 }
157174
175+/** A number in SI for a link: shortest exact-enough decimal. */
176+function formatSI(v: number): string {
177+ return v === 0 ? "0" : String(+v.toPrecision(10));
178+}
179+
158180 export class ParamForm {
159181 readonly params: PhysParams;
160182 sizeUnit: SizeUnit = "in";
183+ /** Height scale of the random shapes drawn by Randomize [m]. */
184+ randomScale = 3e-3;
161185 /** Called after any change; `live` tells whether the field may change during a run. */
162186 onChange: (key: keyof PhysParams, live: boolean) => void = () => {};
163187
@@ -169,8 +193,12 @@ export class ParamForm {
169193 private hitsNote: HTMLElement | null = null;
170194 private dampingNote: HTMLElement | null = null;
171195 private numericsNote: HTMLElement | null = null;
196+ private shapeNote: HTMLElement | null = null;
197+ private featureList: HTMLElement | null = null;
198+ private featureControls: (HTMLInputElement | HTMLSelectElement | HTMLButtonElement)[] = [];
172199 private omegaCache: { key: string; omega: number } | null = null;
173200 private liveKeys = new Set<keyof PhysParams>();
201+ private enabled = true;
174202
175203 constructor(root: HTMLElement, params: PhysParams) {
176204 this.params = params;
@@ -191,6 +219,7 @@ export class ParamForm {
191219 if (f.live) this.liveKeys.add(f.key);
192220 det.appendChild(row);
193221 }
222+ if (g.features) det.appendChild(this.featureSection());
194223 if (g.note) {
195224 const note = document.createElement("p");
196225 note.className = "sub";
@@ -198,6 +227,7 @@ export class ParamForm {
198227 else if (g.note === "hits") this.hitsNote = note;
199228 else if (g.note === "damping") this.dampingNote = note;
200229 else if (g.note === "numerics") this.numericsNote = note;
230+ else if (g.note === "shape") this.shapeNote = note;
201231 else note.textContent = g.note;
202232 det.appendChild(note);
203233 }
@@ -282,6 +312,151 @@ export class ParamForm {
282312 return f.unit === "size" ? SIZE_SCALE[this.sizeUnit] : f.scale;
283313 }
284314
315+ // ---------- Initial-shape feature list ----------
316+
317+ /** The list of features, the add row, and the randomize row. */
318+ private featureSection(): HTMLElement {
319+ const wrap = document.createElement("div");
320+ this.featureList = document.createElement("div");
321+ wrap.appendChild(this.featureList);
322+
323+ const addRow = document.createElement("div");
324+ addRow.className = "feature-row";
325+ const typeSel = document.createElement("select");
326+ for (const t of FEATURE_TYPES) {
327+ const o = document.createElement("option");
328+ o.value = t; o.textContent = FEATURE_SPECS[t].label;
329+ typeSel.appendChild(o);
330+ }
331+ const btnAdd = document.createElement("button");
332+ btnAdd.textContent = "Add";
333+ btnAdd.addEventListener("click", () => {
334+ this.params.features.push(newFeature(typeSel.value as FeatureType));
335+ this.featuresChanged();
336+ });
337+ const btnFlat = document.createElement("button");
338+ btnFlat.textContent = "Clear (flat)";
339+ btnFlat.addEventListener("click", () => {
340+ this.params.features = [];
341+ this.featuresChanged();
342+ });
343+ addRow.append(typeSel, btnAdd, btnFlat);
344+ wrap.appendChild(addRow);
345+
346+ const rndRow = document.createElement("div");
347+ rndRow.className = "feature-row";
348+ const btnRnd = document.createElement("button");
349+ btnRnd.textContent = "Randomize";
350+ btnRnd.title = "Draw a new random combination of features with heights of the order of the scale";
351+ const scaleLabel = document.createElement("label");
352+ scaleLabel.textContent = "scale";
353+ const scaleInp = document.createElement("input");
354+ scaleInp.type = "number";
355+ scaleInp.step = "0.5";
356+ scaleInp.value = formatNumber(this.randomScale * 1e3);
357+ scaleInp.addEventListener("change", () => {
358+ const v = parseFloat(scaleInp.value);
359+ if (isFinite(v) && v > 0) this.randomScale = v * 1e-3;
360+ });
361+ const scaleUnit = document.createElement("span");
362+ scaleUnit.className = "unit";
363+ scaleUnit.textContent = "mm";
364+ btnRnd.addEventListener("click", () => {
365+ this.params.features = randomFeatures(this.params, this.randomScale);
366+ this.featuresChanged();
367+ });
368+ rndRow.append(btnRnd, scaleLabel, scaleInp, scaleUnit);
369+ wrap.appendChild(rndRow);
370+ this.persistentFeatureControls = [typeSel, btnAdd, btnFlat, btnRnd, scaleInp];
371+ this.renderFeatures();
372+ return wrap;
373+ }
374+
375+ private persistentFeatureControls: (HTMLInputElement | HTMLSelectElement | HTMLButtonElement)[] = [];
376+
377+ /** Rebuilds the feature rows from params.features. */
378+ private renderFeatures(): void {
379+ const list = this.featureList;
380+ if (!list) return;
381+ list.replaceChildren();
382+ this.featureControls = [...this.persistentFeatureControls];
383+ const features = this.params.features;
384+ if (features.length === 0) {
385+ const p = document.createElement("p");
386+ p.className = "sub";
387+ p.textContent = "Flat sheet (no features).";
388+ list.appendChild(p);
389+ }
390+ features.forEach((feat, index) => {
391+ const spec = FEATURE_SPECS[feat.type];
392+ const box = document.createElement("div");
393+ box.className = "feature";
394+ box.title = spec.description;
395+ const head = document.createElement("div");
396+ head.className = "feature-head";
397+ const title = document.createElement("span");
398+ title.textContent = `${index + 1}. ${spec.label}`;
399+ const btnRemove = document.createElement("button");
400+ btnRemove.className = "small";
401+ btnRemove.textContent = "Remove";
402+ btnRemove.addEventListener("click", () => {
403+ this.params.features.splice(index, 1);
404+ this.featuresChanged();
405+ });
406+ head.append(title, btnRemove);
407+ box.appendChild(head);
408+ this.featureControls.push(btnRemove);
409+ const rec = feat as unknown as Record<string, number | string>;
410+ for (const fld of spec.fields) {
411+ const row = document.createElement("div");
412+ row.className = "field";
413+ const label = document.createElement("label");
414+ label.textContent = fld.label;
415+ const unit = document.createElement("span");
416+ unit.className = "unit";
417+ unit.textContent = fld.unit;
418+ if (fld.options) {
419+ const sel = document.createElement("select");
420+ for (const [v, t] of fld.options) {
421+ const o = document.createElement("option");
422+ o.value = v; o.textContent = t;
423+ sel.appendChild(o);
424+ }
425+ sel.value = String(rec[fld.key]);
426+ sel.addEventListener("change", () => {
427+ rec[fld.key] = sel.value;
428+ this.featuresChanged(false);
429+ });
430+ row.append(label, sel, unit);
431+ this.featureControls.push(sel);
432+ } else {
433+ const inp = document.createElement("input");
434+ inp.type = "number";
435+ inp.step = String(fld.step);
436+ inp.value = formatNumber((rec[fld.key] as number) / fld.scale);
437+ inp.addEventListener("change", () => {
438+ const v = parseFloat(inp.value);
439+ if (!isFinite(v)) return;
440+ rec[fld.key] = v * fld.scale;
441+ this.featuresChanged(false);
442+ });
443+ row.append(label, inp, unit);
444+ this.featureControls.push(inp);
445+ }
446+ box.appendChild(row);
447+ }
448+ list.appendChild(box);
449+ });
450+ for (const el of this.featureControls) el.disabled = !this.enabled;
451+ }
452+
453+ /** After a change of the feature list: rebuild the rows (unless only a value changed), refresh the note, and notify. */
454+ private featuresChanged(rebuild = true): void {
455+ if (rebuild) this.renderFeatures();
456+ this.updateNotes();
457+ this.onChange("features", false);
458+ }
459+
285460 /** Shows the fields that apply to the selected excitation and hit sequence only. */
286461 private updateVisibility(): void {
287462 for (const c of this.conditional) c.el.hidden = !c.show(this.params);
@@ -296,19 +471,31 @@ export class ParamForm {
296471 }
297472 }
298473 for (const k of selectKeys) (this.inputs.get(k) as HTMLSelectElement).value = String(this.params[k]);
474+ this.renderFeatures();
299475 this.updateVisibility();
300476 this.updateNotes();
301477 }
302478
303479 /** Enables or disables the fields that cannot change while a simulation runs. */
304480 setEnabled(enabled: boolean): void {
481+ this.enabled = enabled;
305482 for (const [key, el] of this.inputs) if (!this.liveKeys.has(key)) el.disabled = !enabled;
483+ for (const el of this.featureControls) el.disabled = !enabled;
306484 }
307485
308- /** Applies URL parameters: any PhysParams key in SI units (booleans as 0/1 or false/true). Returns the keys that were set. */
486+ /**
487+ * Applies URL parameters: any PhysParams key in SI units (booleans as 0/1 or
488+ * false/true), `features` in the text form of shape.ts, and the earlier
489+ * `bowShape` and `bowAmplitude` keys. Returns the keys that were set.
490+ */
309491 applyUrl(q: URLSearchParams): (keyof PhysParams)[] {
310492 const set: (keyof PhysParams)[] = [];
311493 for (const [k, v] of q) {
494+ if (k === "features") {
495+ const f = parseFeatures(v);
496+ if (f) { this.params.features = f; set.push("features"); }
497+ continue;
498+ }
312499 if (!(k in this.params)) continue;
313500 const key = k as keyof PhysParams;
314501 if (selectKeys.includes(key)) {
@@ -325,10 +512,31 @@ export class ParamForm {
325512 }
326513 }
327514 }
515+ if (!set.includes("features") && (q.has("bowShape") || q.has("bowAmplitude"))) {
516+ const amp = parseFloat(q.get("bowAmplitude") ?? "");
517+ this.params.features = legacyFeatures(q.get("bowShape") ?? "width", isFinite(amp) ? amp : 3e-3);
518+ set.push("features");
519+ }
328520 this.refresh();
329521 return set;
330522 }
331523
524+ /** Query string (without the leading ?) of every parameter that differs from its default, in the form applyUrl reads. */
525+ toQuery(): string {
526+ const parts: string[] = [];
527+ const enc = (s: string) => encodeURIComponent(s).replace(/%3A/gi, ":").replace(/%2C/gi, ",").replace(/%3B/gi, ";");
528+ for (const k of Object.keys(defaultParams) as (keyof PhysParams)[]) {
529+ if (k === "features") {
530+ const e = encodeFeatures(this.params.features);
531+ if (e !== encodeFeatures(defaultParams.features)) parts.push(`features=${enc(e)}`);
532+ } else if (this.params[k] !== defaultParams[k]) {
533+ const v = this.params[k];
534+ parts.push(`${k}=${enc(typeof v === "number" ? formatSI(v) : String(v))}`);
535+ }
536+ }
537+ return parts.join("&");
538+ }
539+
332540 gridInfo(): { grid: Grid; inside: number } {
333541 const grid = makeGrid(this.params);
334542 return { grid, inside: countInside(grid) };
@@ -361,6 +569,15 @@ export class ParamForm {
361569 const outline = p.shape === "ellipse" ? "ellipse with a stepped edge" : "rectangle";
362570 this.gridNote.textContent = `${outline}; ${grid.nx} × ${grid.ny} grid, ${inside} nodes in the sheet. Bending partials come out low by 3% near ${(f3 / 1e3).toFixed(1)} kHz and by 10% near ${(f10 / 1e3).toFixed(1)} kHz (grid dispersion).`;
363571 }
572+ if (this.shapeNote) {
573+ const { grid } = this.gridInfo();
574+ const { peakToPeak, maxSlope } = shapeStats(initialShape(p, grid), grid);
575+ const base = "Real sheet metal is never flat. The features are summed and the mean plane is removed; only the slopes matter to the physics. A negative push amount pushes against a bow across the width.";
576+ const stats = p.features.length === 0 ? "" : ` This shape spans ${(peakToPeak * 1e3).toFixed(2)} mm peak to peak with a largest slope of ${maxSlope.toFixed(3)}`;
577+ const warn = maxSlope > SLOPE_LIMIT ? `, beyond the shallow-shell range of about ${SLOPE_LIMIT}; reduce the heights or widen the features.` : stats ? "." : "";
578+ this.shapeNote.textContent = base + stats + warn;
579+ this.shapeNote.className = maxSlope > SLOPE_LIMIT ? "warn" : "sub";
580+ }
364581 if (this.dampingNote) {
365582 const m = makeMaterial(p);
366583 const f = (hz: number) => t60(m, p, 2 * Math.PI * hz).toFixed(2);
@@ -389,7 +606,7 @@ export class ParamForm {
389606 get shape(): PlateShape {
390607 return this.params.shape;
391608 }
392- get bowShape(): BowShape {
393- return this.params.bowShape;
609+ get features(): ShapeFeature[] {
610+ return this.params.features;
394611 }
395612 }
tests/gpu.test.tsmodified+5−3View file
@@ -34,6 +34,7 @@ interface Status {
3434 nx: number;
3535 ny: number;
3636 stick: { maxForce: number } | null;
37+ features?: string;
3738 hits: { done: number; total: number } | null;
3839 }
3940
@@ -82,7 +83,7 @@ async function runToCompletion(query: string): Promise<Status> {
8283 await withPage(`autostart=1&${query}`, async (page, errors) => {
8384 for (let i = 0; i < 600; i++) {
8485 await page.waitForTimeout(250);
85- st = await page.evaluate<Status | null>("window.steelsim ? window.steelsim.status() : null");
86+ st = await page.evaluate<Status | null>("window.steelsim ? { ...window.steelsim.status(), features: window.steelsim.features() } : null");
8687 const noGpu = await page.evaluate<string | null>("window.steelsim.noGpu()");
8788 assert.equal(noGpu, null, `no WebGPU: ${noGpu}`);
8889 if (st && st.phase === "done") break;
@@ -102,8 +103,9 @@ test("a short push-and-release run completes with finite output", { timeout: 180
102103 assert.ok(st.peakPressure > 1e-3 && isFinite(st.peakPressure), "the sheet radiates");
103104 });
104105
105-test("a short drumstick strike on an elliptical sheet completes with finite output", { timeout: 180_000 }, async () => {
106- const st = await runToCompletion("spacing=0.012&duration=0.2&excitation=strike&shape=ellipse");
106+test("a short drumstick strike on an elliptical sheet with dents and a random field completes with finite output", { timeout: 180_000 }, async () => {
107+ const st = await runToCompletion("spacing=0.012&duration=0.2&excitation=strike&shape=ellipse&features=bow:x,0.003,0.3;dent:0.5,0.5,0.03,0.002;noise:7,0.0005,0.05");
108+ assert.equal(st.features, "bow:x,0.003,0.3;dent:0.5,0.5,0.03,0.002;noise:7,0.0005,0.05", "the features were read from the URL");
107109 console.log(` strike: ${st.nx}x${st.ny}, ${(st.stepsPerSec / 1e3).toFixed(0)} k steps/s, peak ${st.peakPressure.toFixed(3)} Pa, max contact force ${st.stick?.maxForce.toFixed(1)} N`);
108110 assert.equal(st.phase, "done");
109111 assert.ok(isFinite(st.wmax) && st.wmax < 0.05, "displacement stays bounded");
tests/shape.test.tsadded+139−0View file
@@ -0,0 +1,139 @@
1+// Tests of the initial-shape features (no GPU needed): the text encoding, the
2+// removal of the mean plane, the geometry of the individual features, and
3+// the random configurations.
4+
5+import { test } from "node:test";
6+import assert from "node:assert/strict";
7+import { defaultParams, makeGrid, nodeWeight, nodeX, nodeY, type Grid, type PhysParams } from "../src/params.ts";
8+import { encodeFeatures, initialShape, parseFeatures, randomFeatures, shapeStats, type ShapeFeature } from "../src/shape.ts";
9+
10+function withFeatures(features: ShapeFeature[], extra: Partial<PhysParams> = {}): PhysParams {
11+ return { ...defaultParams, ...extra, features };
12+}
13+
14+/** Mass-weighted mean, x-moment, and y-moment of a field over the sheet, normalized by the total weight. */
15+function planeMoments(w0: ArrayLike<number>, g: Grid): [number, number, number] {
16+ let sw = 0, s = 0, sx = 0, sy = 0;
17+ for (let j = 0; j < g.ny; j++) {
18+ for (let i = 0; i < g.nx; i++) {
19+ const wgt = nodeWeight(g, i, j);
20+ const v = w0[j * g.nx + i];
21+ sw += wgt; s += wgt * v; sx += wgt * nodeX(g, i) * v; sy += wgt * nodeY(g, j) * v;
22+ }
23+ }
24+ return [s / sw, sx / sw, sy / sw];
25+}
26+
27+const everything = "bow:x,0.003,0.4;bow:y,-0.001,0;twist:0.0008;dent:0.3,0.7,0.03,0.0015;crease:0.6,0.4,25,-1.5,0.015;noise:11,0.0005,0.05";
28+
29+test("features encode to text and back", () => {
30+ const f = parseFeatures(everything);
31+ assert.ok(f);
32+ assert.equal(f.length, 6);
33+ assert.equal(encodeFeatures(f), everything);
34+ assert.deepEqual(parseFeatures("none"), []);
35+ assert.deepEqual(parseFeatures(""), []);
36+ assert.equal(parseFeatures("nonsense"), null);
37+ // An unreadable entry is skipped, the rest is kept.
38+ assert.equal(encodeFeatures(parseFeatures("bow:z,1,0;twist:0.001")!), "twist:0.001");
39+ assert.equal(encodeFeatures(defaultParams.features), "bow:x,0.003,0");
40+});
41+
42+test("the mean plane is removed on the rectangle and on the ellipse", () => {
43+ for (const shape of ["rectangle", "ellipse"] as const) {
44+ const p = withFeatures(parseFeatures(everything)!, { shape });
45+ const g = makeGrid(p);
46+ const w0 = initialShape(p, g);
47+ const [m, mx, my] = planeMoments(w0, g);
48+ const { peakToPeak } = shapeStats(w0, g);
49+ assert.ok(peakToPeak > 3e-3, `${shape}: the shape is not flat (${peakToPeak})`);
50+ assert.ok(Math.abs(m) < 1e-6 * peakToPeak, `${shape}: mean ${m}`);
51+ assert.ok(Math.abs(mx) < 1e-6 * peakToPeak * p.Lx, `${shape}: x tilt ${mx}`);
52+ assert.ok(Math.abs(my) < 1e-6 * peakToPeak * p.Ly, `${shape}: y tilt ${my}`);
53+ // Nodes outside the sheet carry no shape.
54+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (nodeWeight(g, i, j) === 0) assert.equal(w0[j * g.nx + i], 0);
55+ }
56+});
57+
58+test("a bow across the width sags by the given amount and is flat along the height", () => {
59+ const p = withFeatures([{ type: "bow", axis: "x", sag: 3e-3, offset: 0 }]);
60+ const g = makeGrid(p);
61+ const w0 = initialShape(p, g);
62+ const mid = Math.floor(g.ny / 2), ci = (g.nx - 1) / 2;
63+ const center = w0[mid * g.nx + Math.floor(ci)], edge = w0[mid * g.nx];
64+ assert.ok(Math.abs(center - edge - 3e-3 * (1 - ((Math.floor(ci) - ci) * 2 / (g.nx - 1)) ** 2)) < 1e-9, "sag");
65+ for (let j = 0; j < g.ny; j++) assert.ok(Math.abs(w0[j * g.nx + 3] - w0[3]) < 1e-12, "cylindrical");
66+ // The apex offset moves the crest toward the right edge.
67+ const q = withFeatures([{ type: "bow", axis: "x", sag: 3e-3, offset: 0.5 }]);
68+ const v = initialShape(q, g);
69+ let best = 0;
70+ for (let i = 0; i < g.nx; i++) if (v[mid * g.nx + i] > v[mid * g.nx + best]) best = i;
71+ assert.ok(best > ci, "crest right of center");
72+});
73+
74+test("a crease tilts the two sides by half the fold angle", () => {
75+ const fold = 2;
76+ const p = withFeatures([{ type: "crease", x: 0.5, y: 0.5, angle: 90, fold, width: 0 }]);
77+ const g = makeGrid(p);
78+ const w0 = initialShape(p, g);
79+ const { maxSlope } = shapeStats(w0, g);
80+ const expected = Math.tan((fold * Math.PI) / 360);
81+ assert.ok(Math.abs(maxSlope - expected) < 0.05 * expected, `slope ${maxSlope} vs ${expected}`);
82+ // A positive fold is a ridge: the line is the highest part after the plane is removed.
83+ const mid = Math.floor(g.ny / 2);
84+ const onLine = w0[mid * g.nx + Math.floor((g.nx - 1) / 2)];
85+ assert.ok(onLine > w0[mid * g.nx] && onLine > w0[mid * g.nx + g.nx - 1], "ridge");
86+});
87+
88+test("a dent has the requested height and radius", () => {
89+ const p = withFeatures([{ type: "dent", x: 0.5, y: 0.5, radius: 0.04, height: 2e-3 }], { spacing: 0.004 });
90+ const g = makeGrid(p);
91+ const w0 = initialShape(p, g);
92+ const ci = Math.floor((g.nx - 1) / 2), cj = Math.floor((g.ny - 1) / 2);
93+ const center = w0[cj * g.nx + ci];
94+ const far = w0[cj * g.nx];
95+ assert.ok(Math.abs(center - far - 2e-3) < 1e-4, `height ${center - far}`);
96+ const atRadius = w0[cj * g.nx + ci + Math.round(0.04 / g.hx)];
97+ assert.ok(Math.abs((atRadius - far) / (center - far) - Math.exp(-1)) < 0.03, "1/e at the radius");
98+});
99+
100+test("the random field has the requested RMS and depends on the seed but not on the grid", () => {
101+ const p = withFeatures([{ type: "noise", seed: 5, rms: 1e-3, length: 0.06 }]);
102+ const g = makeGrid(p);
103+ const w0 = initialShape(p, g);
104+ let sw = 0, s2 = 0;
105+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) { const wgt = nodeWeight(g, i, j); sw += wgt; s2 += wgt * w0[j * g.nx + i] ** 2; }
106+ const rms = Math.sqrt(s2 / sw);
107+ // The plane removal takes a little of the variance away.
108+ assert.ok(rms > 0.8e-3 && rms <= 1.0e-3 + 1e-9, `rms ${rms}`);
109+ const other = initialShape(withFeatures([{ type: "noise", seed: 6, rms: 1e-3, length: 0.06 }]), g);
110+ let diff = 0;
111+ for (let k = 0; k < g.n; k++) diff = Math.max(diff, Math.abs(other[k] - w0[k]));
112+ assert.ok(diff > 0.5e-3, "another seed gives another field");
113+ // The same seed on a finer grid gives the same field at the shared nodes.
114+ const fine = makeGrid({ ...p, spacing: 0.004 });
115+ const wf = initialShape(p, fine);
116+ assert.equal(fine.nx, 2 * g.nx - 1);
117+ let err = 0;
118+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) err = Math.max(err, Math.abs(wf[2 * j * fine.nx + 2 * i] - w0[j * g.nx + i]));
119+ assert.ok(err < 2e-5, `grid dependence ${err}`);
120+});
121+
122+test("random configurations are readable and stay within the shallow-shell range", () => {
123+ let seed = 12345;
124+ const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
125+ const types = new Set<string>();
126+ for (let trial = 0; trial < 40; trial++) {
127+ const features = randomFeatures(defaultParams, 3e-3, rnd);
128+ for (const f of features) types.add(f.type);
129+ assert.equal(features[0].type, "bow");
130+ const text = encodeFeatures(features);
131+ assert.equal(encodeFeatures(parseFeatures(text)!), text);
132+ const p = withFeatures(features);
133+ const g = makeGrid(p);
134+ const { peakToPeak, maxSlope } = shapeStats(initialShape(p, g), g);
135+ assert.ok(peakToPeak > 0.5e-3 && peakToPeak < 15e-3, `peak to peak ${peakToPeak}`);
136+ assert.ok(maxSlope < 0.1, `slope ${maxSlope}`);
137+ }
138+ for (const t of ["bow", "twist", "dent", "crease", "noise"]) assert.ok(types.has(t), `${t} was drawn`);
139+});