Plate: remove the curvature (domed rest shape) and its cross-section, to keep the plate simple; Parameters stays in the right column
5 changed files+58−235
src/plate/README.mdmodified+3−3View file
@@ -1,6 +1,6 @@
11 # plate
22
3-Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices); step 5 (2026-10-09), clamped edges; step 6 (2026-10-09), Lyapunov exponent; step 7 (2026-10-09), curvature.
3+Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices); step 5 (2026-10-09), clamped edges; step 6 (2026-10-09), Lyapunov exponent.
44
55 Current state: N × N masses (N per side, 2 to 32, default 8) on vertical rails at the nodes of a square grid, spacing a, moving up and down only. Bending κ on the five-point Laplacian at each node, energy κ (L u)² / 2, force −κ L(L u). The edge nodes are supports at u = 0 with no bending term (mirror ghost), which makes the edge simply supported. Stretching (von Kármán), below, set by the thickness h (default 5 cm; 0 = linear). Damping as in the beam: drag (τ, all modes alike) and a frequency-dependent loss (dashpots alongside the grid springs, force −c L v, set by τ top at the top of the band). Parameters are given for the 8 × 8 plate (m 1 kg, κ 20 N/m, h 0.05 m, a 0.5 m, so the side is 4.5 m; τ 540 s, τ top 200 s); Masses per side refines the same plate (`refine`: m with the cell area, κ inversely, the loss coefficient c unchanged). Strike: a kick (pyramid velocity bump of half-width 0.5 m) or the shared Hertz mallet (`../common/mallet.ts`, M 1 kg, K 1000) over the same bump, at (x, y) m from the lower left corner (default (1.3, 1.1)); clicking the plate moves the strike point. The card shows the plate from above as a height map (red up, blue down, scaled to the recent largest |u|), with the strike point (cross) and the listening point (circle). The sound panel is the beam's, reduced: velocity, displacement or acceleration at a listening point (default (3.1, 1.7), bilinear between masses), pitch of the lowest mode, the whole decay as a note through the shared renderer and anti-alias filter, spectrogram with dashed lines at the distinct mode pitches, energy chart, Replay and Save WAV. Not yet: ghost, mode bars and histories, cross-section.
66
@@ -28,6 +28,6 @@ Cost issue, fixed (2026-10-09, `../common/note-render.ts`): the 4x oversampled n
2828
2929 Lyapunov exponent (step 6, 2026-10-09): as the beam's (two-trajectory method, `../common/lyapunov.ts`), with the plate's shadow in `lyapunov.ts`: a second plate D0 = 1e-8 away, distances with displacements in units of A/√(N²) (A the lowest mode's amplitude at which stretching equals bending, `nonlinearAmplitude`; 1 m for a linear plate) and momenta in m ω₁ times that; the initial offset is a fixed pseudo-random pattern (free: its rigid part removed). Layout now as the beam's: Parameters folded under the plate, a Real time card on the right with Lyapunov (λ 10 s and since the strike, first 5 s left out of the y range) and Energy error charts. Sound panel: Lyapunov chart and λ vs E scatter, both off at first (the shadow doubles the time to compute a note: default note 4.3 s without, about 12 s with, in the browser). Exponents (2026-10-09, N = 8, kick at (1.3, 1.1), undamped, 300 s, λ since the strike, 1/s): linear plate 0.007 at every energy (finite-time; the last 100 s, −0.001: zero); simply supported, h 5 cm: 0.008, 0.065, 0.26, 0.53, 1.01 at 0.001, 0.1, 1, 10, 100 J; free, h 5 cm: 0.006, 0.015, 0.16, 1.54, 8.4. The free plate's exponent grows much faster with energy.
3030
31-Curvature (step 7, 2026-10-09; Parameters, rise, default 0): the plate rests, free of stress, in a shallow dome u0 of height H (held edges: H sin(πx/L) sin(πy/L), zero at the supports; free: the paraboloid H(1 − r²/c²) about the centre, a spherical cap), as the beam's curved rest shape. Marguerre shallow shell: bending unchanged (linear, from the rest shape), stretching from the rest shape, s = l(u0 + u, u0 + u) − l(u0, u0) = l(u, u) + 2 l(u0, u), the force its exact gradient (the second differences of the total shape in the scatter). The 2 l(u0, u) term adds linear stiffness, which lifts the low modes most (the cymbal effect), and a quadratic coupling. Acts only through the stretching (h > 0). With a rise the modes are computed (symmetry blocks; the rest shapes are mirror-symmetric), also simply supported, from the linearized forces (`forces(..., linear)` now includes the curvature stiffness; `linearEnergy` = bending + that quadratic part). Nonlinear scale and amplitude generalized: the lowest mode's amplitude at which the nonlinear stretching (cubic and quartic) equals the linear energy, by bisection (flat: unchanged, 0.131 J). Checks: force against central differences 2e-10 to 5e-10 (simply supported, clamped, free, H 0.1 m, N 6); linearized force against the force at 1e-7 amplitude 6e-8; free plate still invariant under rigid motion (2e-15); flat results unchanged. Continuum check (simply supported, H/h = 0.5, 1, 2, 4): the (1, 1) shape's Rayleigh quotient ω²/ω₀² is 1 + 0.836 (H/h)² in the continuum (single mode, exact stress solution, the same 0.209 coefficient as the flat check above, times 4); grid 1.2102, 1.8410, 4.3639, 14.456 at N = 32 against 1.2090, 1.8360, 4.3440, 14.376 (second order). Lowest modes, simply supported, N = 16, relative to the flat (1, 1): H = 0.05 m (= h): 1.35, 2.64, 2.64, 4.04; H = 0.1: 1.94, 3.03, 3.03, 4.24; H = 0.2: 2.90, 4.05, 4.05, 4.89 (flat: 1, 2.48, 2.48, 3.97). Energy error with a rise of 0.1 m at 100 J: 1e-3 (simply supported), 6e-4 (free).
31+Layout (2026-10-09): Parameters back in the right column, above Real time.
3232
33-Layout and cross-section (2026-10-09): Parameters back in the right column (above Real time). Under the height map, a cross-section through the middle of the plate (y = side / 2, dashed line on the map) at true scale (the same metres per pixel across and up): the rest shape faint, the plate now (rest shape plus displacement) dark, each a band of the plate's thickness h, supports marked for held edges; the label gives the side, h, and the rise as multiples of h and fractions of the side.
33+Curvature (2026-10-09, e7f0b54, removed the same day at Jeremy's request, to keep the plate simple): an optional domed rest shape (rise; Marguerre shallow shell, stretching measured from the rest shape) with a true-scale cross-section under the height map. It worked (the (1, 1) stiffening matched the continuum 1 + 0.836 (rise/h)^2 at second order), and on the cymbal-matched square a 3 cm dome removed the flat plate's low rumble but rang at the uniform dome's membrane frequency (830 Hz). See git history to restore it.
src/plate/index.htmlmodified+2−4View file
@@ -11,7 +11,6 @@
1111 <div class="grid">
1212 <section class="card">
1313 <canvas id="plate"></canvas>
14- <canvas id="section" title="Cross-section through the middle of the plate (dashed line above), at true scale"></canvas>
1514 <div id="controls">
1615 <div class="row">
1716 <button id="strike" class="primary">Strike</button>
@@ -59,7 +58,6 @@
5958 <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
6059 <label>κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
6160 <label>h <input id="p-thick" type="number" step="any" /> <span class="unit">m (thickness, sets the stretching; 0 = linear)</span></label>
62- <label>rise <input id="p-rise" type="number" step="any" /> <span class="unit">m (height of a domed rest shape; 0 = flat)</span></label>
6361 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
6462 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
6563 <label>τ top <input id="p-tautop" type="number" step="any" /> <span class="unit">s (decay of the highest modes from the loss, 0 = off)</span></label>
@@ -121,10 +119,10 @@
121119 </template>
122120 <template id="help-params">
123121 <h2>Parameters</h2>
124- <p>Mass, bending κ, the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the rise of a domed rest shape (0 for flat; it acts through the stretching, so it needs h > 0), the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
122+ <p>Mass, bending κ, the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
125123 <p>The values are for the plate with 8 masses per side. Masses per side refines the same plate: same size, strike energies and decay, with each mass scaled with the area of its cell and κ inversely (the thickness stays), so the low modes stay put and more masses add higher ones.</p>
126124 <p>The stretching follows the von Kármán plate equations: the in-plane stress comes from a stress function, solved on the grid at every step. A simply supported edge may move outward but not along itself; a free edge is free in its own plane too. The nonlinear scale below is the energy of the lowest mode at which its stretching energy equals its bending energy.</p>
127- <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. Clamped edges hold the slope at zero too; the modes are computed (no (m, n)), and the plate is stiffer than a simply supported one, with overtones about as widely spread (plate theory: 1, 2.04, 3.01, 3.66, ...). A rise curves the plate at rest into a shallow dome (a shell, as a cymbal is): the plate rests in it without stress, and is drawn as its displacement from it. Bending is unchanged, but bending the dome stretches it in its own plane, which stiffens the plate: the lowest modes rise most (by about √(1 + 0.84 (rise / h)²) for the simply supported plate's lowest), so the low end of the spectrum crowds together, and the stretching gains a quadratic coupling besides the cubic one. With a rise the modes are computed. For a flat simply supported plate, the normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
125+ <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. Clamped edges hold the slope at zero too; the modes are computed (no (m, n)), and the plate is stiffer than a simply supported one, with overtones about as widely spread (plate theory: 1, 2.04, 3.01, 3.66, ...). For a simply supported plate, the normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
128126 </template>
129127 <template id="help-realtime">
130128 <h2>Real time</h2>
src/plate/main.tsmodified+7−79View file
@@ -2,7 +2,7 @@ import "../common/style.css";
22 import "./style.css";
33 import {
44 Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, nodePosition, nonlinearScale, plateSide, refine,
5- masses, removeRigid, restShape, stepFor, strikeAt, valueAt, type Edges, type Params, type State,
5+ masses, removeRigid, stepFor, strikeAt, valueAt, type Edges, type Params, type State,
66 } from "./plate";
77 import { fit } from "../common/canvas";
88 import { newCursor } from "../common/cursor";
@@ -42,7 +42,7 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
4242 // half an octave, 100 J is a broadband wash). base: the page's
4343 // parameters, for the plate at REF_N masses per side; P: the same plate
4444 // refined to the chosen N (see refine).
45-let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, rise: 0, a: 0.5, tau: 540, tauTop: 200 };
45+let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
4646 let P: Params = base;
4747 // The mallet: a physical object, not refined with N.
4848 let mallet: MalletParams = { M: 1, K: 1000 };
@@ -185,7 +185,6 @@ const paramInputs = {
185185 m: document.getElementById("p-m") as HTMLInputElement,
186186 kappa: document.getElementById("p-kappa") as HTMLInputElement,
187187 thick: document.getElementById("p-thick") as HTMLInputElement,
188- rise: document.getElementById("p-rise") as HTMLInputElement,
189188 a: document.getElementById("p-a") as HTMLInputElement,
190189 tau: document.getElementById("p-tau") as HTMLInputElement,
191190 tauTop: document.getElementById("p-tautop") as HTMLInputElement,
@@ -217,7 +216,7 @@ function distinctModes(): { w: number; count: number; k: number }[] {
217216 // and the decay times.
218217 function showDerived() {
219218 const w1 = modes.omega[0], hz = (w: number) => w / (2 * Math.PI);
220- const free = modes.S !== null; // computed modes: no (m, n)
219+ const free = P.edges !== "simply"; // no (m, n)
221220 const rows = distinctModes().slice(0, 8).map(({ w, count, k }) => {
222221 const m = modes.m[k], n = modes.n[k];
223222 return (
@@ -236,8 +235,8 @@ function showDerived() {
236235 (P.N === REF_N ? "" : `at ${P.N} × ${P.N}: m ${g(P.m)} kg, κ ${g(P.kappa)} N/m, a ${g(P.a)} m\n`) +
237236 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
238237 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
239- `lowest modes${P.edges === "free" ? " (free edges, rigid motions removed)" : free ? " (computed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
240- `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: nonlinear stretching = linear energy)\n` +
238+ `lowest modes${P.edges === "free" ? " (free edges, rigid motions removed)" : P.edges === "clamped" ? " (clamped, computed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
239+ `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: stretching = bending energy)\n` +
241240 `decay time ${decay}\n` +
242241 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
243242 `${(1000 * stepFor(P, modes, 0)).toPrecision(3)} ms when weak, ` +
@@ -442,10 +441,10 @@ edgesSelect.addEventListener("change", () => setParams({ ...base, edges: edgesSe
442441 document.getElementById("param-apply")!.addEventListener("click", () => {
443442 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
444443 const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
445- const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick", "rise"].includes(key) && x === 0))) &&
444+ const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick"].includes(key) && x === 0))) &&
446445 Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
447446 if (!ok) {
448- paramStatus.textContent = "parameters must be positive (h, rise, τ, τ top may be 0)";
447+ paramStatus.textContent = "parameters must be positive (h, τ, τ top may be 0)";
449448 return;
450449 }
451450 paramStatus.textContent = "";
@@ -515,15 +514,6 @@ function drawPlate() {
515514 ctx.setLineDash(free ? [4, 4] : []);
516515 ctx.strokeRect(g.x0, g.y0, g.size, g.size);
517516 ctx.setLineDash([]);
518- // Where the cross-section below is taken.
519- ctx.strokeStyle = "rgba(0, 0, 0, 0.35)";
520- ctx.lineWidth = 1;
521- ctx.setLineDash([3, 4]);
522- ctx.beginPath();
523- ctx.moveTo(px(0), py(L / 2));
524- ctx.lineTo(px(L), py(L / 2));
525- ctx.stroke();
526- ctx.setLineDash([]);
527517 // Strike point and listening point.
528518 const sx = px(coord(strikeXInput)), sy = py(coord(strikeYInput));
529519 ctx.strokeStyle = TEXT;
@@ -550,67 +540,6 @@ function drawPlate() {
550540 ctx.textAlign = "left";
551541 }
552542
553-// Cross-section through the middle of the plate (y = side / 2) at true
554-// scale, the same metres per pixel across and up, so the dome's rise can be
555-// seen against the plate's size and thickness: the rest shape (faint) and
556-// the plate now (rest shape plus displacement), each drawn as a band of the
557-// plate's thickness h (at least a hairline), the edge supports as small
558-// marks. Shapes taller than the strip are clipped.
559-const sectionCanvas = document.getElementById("section") as HTMLCanvasElement;
560-function drawSection() {
561- const ctx = fit(sectionCanvas);
562- const W = sectionCanvas.clientWidth, Hc = sectionCanvas.clientHeight;
563- ctx.clearRect(0, 0, W, Hc);
564- const g = geometry(), L = plateSide(P), sc = g.size / L; // px per m
565- const x0 = g.x0, base = Hc - 26; // z = 0
566- const u0 = restShape(P), y = L / 2, n = 240;
567- const zRest = (x: number) => (u0 ? valueAt(P, (k) => u0[k], x, y) : 0);
568- const band = (z: (x: number) => number, fill: string) => {
569- const half = Math.max(0.5, (P.thick * sc) / 2);
570- ctx.fillStyle = fill;
571- ctx.beginPath();
572- for (let i = 0; i <= n; i++) {
573- const x = (L * i) / n;
574- ctx.lineTo(x0 + x * sc, base - z(x) * sc - half);
575- }
576- for (let i = n; i >= 0; i--) {
577- const x = (L * i) / n;
578- ctx.lineTo(x0 + x * sc, base - z(x) * sc + half);
579- }
580- ctx.closePath();
581- ctx.fill();
582- };
583- ctx.save();
584- ctx.beginPath();
585- ctx.rect(0, 0, W, Hc - 14);
586- ctx.clip();
587- band(zRest, "rgba(122, 125, 130, 0.35)");
588- band((x) => zRest(x) + valueAt(P, (k) => s.u[k], x, y), TEXT);
589- ctx.restore();
590- // Supports of a held plate.
591- if (P.edges !== "free") {
592- ctx.fillStyle = MUTED;
593- for (const xs of [x0, x0 + g.size]) {
594- ctx.beginPath();
595- ctx.moveTo(xs, base);
596- ctx.lineTo(xs - 4, base + 7);
597- ctx.lineTo(xs + 4, base + 7);
598- ctx.closePath();
599- ctx.fill();
600- }
601- }
602- ctx.fillStyle = MUTED;
603- ctx.font = "11px system-ui, sans-serif";
604- ctx.textAlign = "left";
605- const rel = (x: number) => (x >= 10 ? x.toFixed(0) : x.toPrecision(2));
606- ctx.fillText(
607- `true scale · side ${fmtLength(L)}` +
608- (P.thick > 0 ? ` · h ${fmtLength(P.thick)}` : "") +
609- (P.rise > 0 ? ` · rise ${fmtLength(P.rise)}${P.thick > 0 ? ` = ${rel(P.rise / P.thick)} h` : ""} = side / ${rel(L / P.rise)}` : " · flat"),
610- x0, Hc - 2,
611- );
612-}
613-
614543 // Animation: the simulation advances `speed` seconds per second of real time.
615544 let last = performance.now();
616545 let carry = 0;
@@ -637,7 +566,6 @@ function frame(now: number) {
637566 for (const u of s.u) umax = Math.max(umax, Math.abs(u));
638567 ampl = Math.max(umax, 0.998 * ampl);
639568 drawPlate();
640- drawSection();
641569 traces.draw();
642570 if (note && soundPanel.open) {
643571 const at = player.position();
src/plate/plate.tsmodified+46−148View file
@@ -38,19 +38,6 @@
3838 // would need another solver. The stretching is set by the thickness h:
3939 // ks / kappa = 12 (1 - nu^2) / h^2 (from D = Eh^3 / 12 (1 - nu^2), nu = 0.3),
4040 // so it matters once the plate moves by about h.
41-// Curvature (rise H > 0): the plate rests, free of stress, in a dome u0 (a
42-// shallow shell, Marguerre's equations): held edges H sin(pi x / L) sin(pi
43-// y / L), zero at the supports; free, the paraboloid H (1 - r^2 / c^2)
44-// about the centre, c the half side (a spherical cap). Displacements u are
45-// measured from it. Bending is linear, so it is unchanged; the stretching is
46-// measured from the rest shape: s = l(u0 + u, u0 + u) - l(u0, u0) = l(u, u) +
47-// 2 l(u0, u), as the beam's e = d^2 + 2 D0 d. The term 2 l(u0, u) is linear
48-// in u: it adds linear stiffness, raising the low modes (in a cymbal, the
49-// dome lifts the lowest modes to a few hundred Hz), and a quadratic
50-// coupling besides the cubic one. Curvature acts only through the
51-// stretching: with h = 0 (no stretching) the rise has no effect. With a rise
52-// the modes are computed (L^2 is no longer diagonal in the coupled
53-// stiffness), also when simply supported.
5441 // Clamped: the supports also hold the slope at zero, an even ghost beyond
5542 // them (u(-2, j) = u(0, j)). The support nodes then have a bending term,
5643 // (L u) = -2 u(0, j) there, weighted 1/2 (the support's half cell): the
@@ -111,7 +98,6 @@ export interface Params {
11198 m: number; // kg, each mass
11299 kappa: number; // N/m, bending
113100 thick: number; // m, thickness, sets the stretching (von Karman); 0 = linear
114- rise: number; // m, height of the curved rest shape (a dome); 0 = flat
115101 a: number; // m, rest spacing (for positions and drawing; it does not enter the equations)
116102 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
117103 tauTop: number; // s, amplitude decay time of the top of the band from the loss alone; 0 = off
@@ -427,24 +413,19 @@ function inverseClamped(M: number, s: Float64Array, phi: Float64Array): void {
427413 for (let k = 0; k < M * M; k++) phi[k] = x[k] - w[k];
428414 }
429415
430-// Second differences of a free plate at its interior nodes (index q), from u.
431-function diffsFree(N: number, u: Float64Array, xx: Float64Array, yy: Float64Array, xy: Float64Array): void {
416+// Stretching energy of a free plate at u; with F given, its force is added.
417+function freeStretching(P: Params, ks: number, u: Float64Array, F?: Float64Array): number {
418+ const N = P.N, M = N - 2, z = stretchScratch(M);
419+ const { xx, yy, xy, s, phi } = z;
432420 for (let j = 1, q = 0; j < N - 1; j++) {
433421 for (let i = 1, k = 1 + N * j; i < N - 1; i++, k++, q++) {
434422 const cc = 2 * u[k];
435423 xx[q] = u[k - 1] - cc + u[k + 1];
436424 yy[q] = u[k - N] - cc + u[k + N];
437425 xy[q] = 0.25 * (u[k + N + 1] - u[k - N + 1] - u[k + N - 1] + u[k - N - 1]);
426+ s[q] = 2 * (xx[q] * yy[q] - xy[q] * xy[q]);
438427 }
439428 }
440-}
441-
442-// Stretching energy of a free plate at u; with F given, its force is added.
443-function freeStretching(P: Params, ks: number, u: Float64Array, F: Float64Array | undefined, linear: boolean): number {
444- const N = P.N, M = N - 2, z = stretchScratch(M);
445- const { xx, yy, xy, s, phi } = z;
446- diffsFree(N, u, xx, yy, xy);
447- if (!sourceTerm(P, M * M, xx, yy, xy, s, linear)) return 0;
448429 inverseClamped(M, s, phi);
449430 let V = 0;
450431 for (let q = 0; q < M * M; q++) V += s[q] * phi[q];
@@ -470,75 +451,6 @@ function freeStretching(P: Params, ks: number, u: Float64Array, F: Float64Array
470451 return (ks / 8) * V;
471452 }
472453
473-// The rest shape u0 (m) of each node; null when flat. Cached per Params.
474-const restCache = new WeakMap<Params, Float64Array | null>();
475-export function restShape(P: Params): Float64Array | null {
476- if (restCache.has(P)) return restCache.get(P)!;
477- let u0: Float64Array | null = null;
478- if (P.rise) {
479- const N = P.N, L = plateSide(P), c = L / 2;
480- u0 = new Float64Array(N * N);
481- for (let j = 0; j < N; j++) {
482- for (let i = 0; i < N; i++) {
483- const x = nodePosition(P, i), y = nodePosition(P, j);
484- u0[i + N * j] = isFree(P)
485- ? P.rise * (1 - ((x - c) ** 2 + (y - c) ** 2) / (c * c))
486- : P.rise * Math.sin((Math.PI * x) / L) * Math.sin((Math.PI * y) / L);
487- }
488- }
489- }
490- restCache.set(P, u0);
491- return u0;
492-}
493-// Second differences of the rest shape at the stretching's nodes, and its
494-// own s0 = l(u0, u0); null when flat.
495-interface RestDiffs { xx: Float64Array; yy: Float64Array; xy: Float64Array; s0: Float64Array }
496-const restDiffCache = new WeakMap<Params, RestDiffs | null>();
497-function restDiffs(P: Params): RestDiffs | null {
498- if (restDiffCache.has(P)) return restDiffCache.get(P)!;
499- const u0 = restShape(P);
500- let r: RestDiffs | null = null;
501- if (u0) {
502- const N = P.N, n = isFree(P) ? (N - 2) * (N - 2) : N * N, f = () => new Float64Array(n);
503- r = { xx: f(), yy: f(), xy: f(), s0: f() };
504- if (isFree(P)) diffsFree(N, u0, r.xx, r.yy, r.xy);
505- else diffsHeld(N, u0, new Float64Array((N + 2) * (N + 2)), r.xx, r.yy, r.xy);
506- for (let k = 0; k < n; k++) r.s0[k] = 2 * (r.xx[k] * r.yy[k] - r.xy[k] * r.xy[k]);
507- }
508- restDiffCache.set(P, r);
509- return r;
510-}
511-// The stretching's source s from the second differences of u (n values, in
512-// place): l(u, u) when flat; with a rest shape, l(u0 + u, u0 + u) - l(u0,
513-// u0), leaving the total shape's differences in xx, yy, xy for the force;
514-// linear: only the part linear in u, 2 l(u0, u), leaving the rest shape's
515-// differences (the force of the linearized energy). Returns false when s is
516-// zero (linear and flat).
517-function sourceTerm(P: Params, n: number, xx: Float64Array, yy: Float64Array, xy: Float64Array, s: Float64Array, linear: boolean): boolean {
518- const r = restDiffs(P);
519- if (!r) {
520- if (linear) return false;
521- for (let k = 0; k < n; k++) s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
522- return true;
523- }
524- if (linear) {
525- for (let k = 0; k < n; k++) {
526- s[k] = 2 * (r.xx[k] * yy[k] + r.yy[k] * xx[k] - 2 * r.xy[k] * xy[k]);
527- xx[k] = r.xx[k];
528- yy[k] = r.yy[k];
529- xy[k] = r.xy[k];
530- }
531- return true;
532- }
533- for (let k = 0; k < n; k++) {
534- xx[k] += r.xx[k];
535- yy[k] += r.yy[k];
536- xy[k] += r.xy[k];
537- s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]) - r.s0[k];
538- }
539- return true;
540-}
541-
542454 // Scratch for the stretching terms, by N. Padded arrays have a ring of
543455 // zeros beyond the edge, (N + 2)^2 values, node (i, j) at (i + 1) + (N + 2)(j + 1).
544456 interface Stretch { pad: Float64Array; xx: Float64Array; yy: Float64Array; xy: Float64Array; s: Float64Array; phi: Float64Array; a: Float64Array; b: Float64Array; c: Float64Array }
@@ -552,32 +464,25 @@ function stretchScratch(N: number): Stretch {
552464 }
553465 return z;
554466 }
555-// Second differences of a held plate at every node (zero beyond the edge):
556-// xx, yy and the centred cross difference xy (pad: scratch, (N + 2)^2).
557-function diffsHeld(N: number, u: Float64Array, pad: Float64Array, xx: Float64Array, yy: Float64Array, xy: Float64Array): void {
558- const R = N + 2;
467+// Stretching energy (J) at u; with F given, its force is added to F.
468+function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
469+ const ks = stretchCoefficient(P);
470+ if (ks === 0) return 0;
471+ if (isFree(P)) return freeStretching(P, ks, u, F);
472+ const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
473+ const { pad, xx, yy, xy, s, phi, a, b, c } = z;
559474 for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) pad[i + 1 + R * (j + 1)] = u[i + N * j];
475+ // Second differences at every node (zero beyond the edge): xx, yy and the
476+ // centred cross difference xy; s = l(u, u).
560477 for (let j = 0, k = 0; j < N; j++) {
561478 for (let i = 0, K = 1 + R * (j + 1); i < N; i++, k++, K++) {
562479 const cc = 2 * pad[K];
563480 xx[k] = pad[K - 1] - cc + pad[K + 1];
564481 yy[k] = pad[K - R] - cc + pad[K + R];
565482 xy[k] = 0.25 * (pad[K + R + 1] - pad[K - R + 1] - pad[K + R - 1] + pad[K - R - 1]);
483+ s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
566484 }
567485 }
568-}
569-
570-// Stretching energy (J) at u; with F given, its force is added to F.
571-// linear: only the part quadratic in u (from 2 l(u0, u); zero when flat),
572-// and its (linear) force.
573-function stretching(P: Params, u: Float64Array, F?: Float64Array, linear = false): number {
574- const ks = stretchCoefficient(P);
575- if (ks === 0) return 0;
576- if (isFree(P)) return freeStretching(P, ks, u, F, linear);
577- const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
578- const { pad, xx, yy, xy, s, phi, a, b, c } = z;
579- diffsHeld(N, u, pad, xx, yy, xy);
580- if (!sourceTerm(P, n, xx, yy, xy, s, linear)) return 0;
581486 inverseL2(N, s, phi);
582487 let V = 0;
583488 for (let k = 0; k < n; k++) V += s[k] * phi[k];
@@ -606,12 +511,11 @@ function stretching(P: Params, u: Float64Array, F?: Float64Array, linear = false
606511 }
607512
608513 // Forces: bending, F = -kappa L (L u) (free: freeBending), plus stretching;
609-// c (if given) is scratch. linear: the linearized forces (bending, and the
610-// curvature's linear stretching stiffness).
514+// c (if given) is scratch. linear: bending only.
611515 export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length), linear = false): void {
612516 if (isFree(P)) {
613517 freeBending(P, u, F);
614- stretching(P, u, F, linear);
518+ if (!linear) stretching(P, u, F);
615519 return;
616520 }
617521 laplacian(P.N, u, c);
@@ -621,7 +525,7 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Ar
621525 const N = P.N;
622526 for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) F[k] -= 2 * P.kappa * touched(N, i, j) * u[k];
623527 }
624- stretching(P, u, F, linear);
528+ if (!linear) stretching(P, u, F);
625529 }
626530
627531 export function kinetic(P: Params, s: State): number {
@@ -630,8 +534,8 @@ export function kinetic(P: Params, s: State): number {
630534 for (let k = 0; k < s.p.length; k++) T += (s.p[k] * s.p[k]) / M[k];
631535 return T / 2;
632536 }
633-// Bending energy at u.
634-export function bendingEnergy(P: Params, u: Float64Array): number {
537+// Linear energy at u: bending.
538+export function linearEnergy(P: Params, u: Float64Array): number {
635539 if (isFree(P)) return freeBending(P, u);
636540 const N = P.N, c = new Float64Array(u.length);
637541 laplacian(N, u, c);
@@ -640,11 +544,8 @@ export function bendingEnergy(P: Params, u: Float64Array): number {
640544 if (isClamped(P)) for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) V += 2 * touched(N, i, j) * u[k] * u[k];
641545 return (P.kappa * V) / 2;
642546 }
643-// Linear energy at u (the quadratic part of the potential): bending, plus
644-// the curvature's linear stretching stiffness.
645-export const linearEnergy = (P: Params, u: Float64Array) => bendingEnergy(P, u) + stretching(P, u, undefined, true);
646547 export const stretchEnergy = (P: Params, u: Float64Array) => stretching(P, u);
647-export const potential = (P: Params, s: State) => bendingEnergy(P, s.u) + stretching(P, s.u);
548+export const potential = (P: Params, s: State) => linearEnergy(P, s.u) + stretching(P, s.u);
648549 export const energy = (P: Params, s: State) => kinetic(P, s) + potential(P, s);
649550
650551 // Velocity Verlet, with drag and loss in half steps around it (Strang), as
@@ -741,7 +642,7 @@ export class Modes {
741642 constructor(P: Params) {
742643 const N = P.N;
743644 this.N = N;
744- if (isFree(P) || isClamped(P) || restDiffs(P) && stretchCoefficient(P) > 0) {
645+ if (isFree(P) || isClamped(P)) {
745646 const { lambda, S } = computedModes(P);
746647 this.count = lambda.length;
747648 this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0)));
@@ -838,36 +739,33 @@ export function decayTimes(P: Params, modes: Modes): Float64Array {
838739 });
839740 }
840741
841-// Nonlinear amplitude (m) and scale (J): the lowest mode's shape u1 (unit
842-// length, sum of u^2 = 1 m^2) at amplitude q has linear energy V2(q) (bending
843-// and the curvature's stiffness, ~ q^2) and nonlinear stretching energy
844-// Vn(q) (the rest, ~ q^3 and q^4); the amplitude is the q at which |Vn| =
845-// V2, the scale V2 there. Flat: q^2 = V2(1) / Vn(1). Infinity when linear.
846-function nonlinearPoint(P: Params, modes: Modes): { q: number; E: number } {
847- if (stretchCoefficient(P) === 0) return { q: Infinity, E: Infinity };
742+// Nonlinear scale (J): the energy of the lowest mode ((1, 1) when simply
743+// supported) at which its stretching energy equals its linear energy
744+// (Infinity when linear). Bending goes as the amplitude squared, stretching
745+// as its fourth power.
746+export function nonlinearScale(P: Params, modes: Modes): number {
747+ const { V2, V4 } = lowestModeEnergies(P, modes);
748+ return V4 > 0 ? (V2 * V2) / V4 : Infinity;
749+}
750+// Bending (V2) and stretching (V4) energy of the lowest mode's shape scaled
751+// to unit length (sum of u^2 = 1 m^2); 0 stretching when linear.
752+function lowestModeEnergies(P: Params, modes: Modes): { V2: number; V4: number } {
848753 const N = P.N;
849- let u1 = new Float64Array(N * N);
850- if (modes.S) u1 = Float64Array.from(modes.S[0]);
851- else for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u1[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
754+ let u = new Float64Array(N * N);
755+ if (modes.S) u = Float64Array.from(modes.S[0]);
756+ else for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
852757 let n2 = 0;
853- for (const x of u1) n2 += x * x;
854- for (let k = 0; k < u1.length; k++) u1[k] /= Math.sqrt(n2);
855- const at = (q: number) => {
856- const u = u1.map((x) => q * x), V2 = linearEnergy(P, u);
857- return { V2, ratio: Math.abs(stretching(P, u) - stretching(P, u, undefined, true)) / V2 };
858- };
859- // Geometric bisection on q between 1e-9 and 1e3 m.
860- let lo = 1e-9, hi = 1e3;
861- if (at(hi).ratio < 1) return { q: Infinity, E: Infinity };
862- for (let it = 0; it < 80; it++) {
863- const mid = Math.sqrt(lo * hi);
864- if (at(mid).ratio < 1) lo = mid;
865- else hi = mid;
866- }
867- return { q: hi, E: at(hi).V2 };
758+ for (const x of u) n2 += x * x;
759+ for (let k = 0; k < u.length; k++) u[k] /= Math.sqrt(n2);
760+ return { V2: linearEnergy(P, u), V4: stretching(P, u) };
761+}
762+// Nonlinear amplitude (m): the length of the lowest mode's displacement
763+// vector at which its stretching energy equals its bending energy (Infinity
764+// when linear).
765+export function nonlinearAmplitude(P: Params, modes: Modes): number {
766+ const { V2, V4 } = lowestModeEnergies(P, modes);
767+ return V4 > 0 ? Math.sqrt(V2 / V4) : Infinity;
868768 }
869-export const nonlinearScale = (P: Params, modes: Modes) => nonlinearPoint(P, modes).E;
870-export const nonlinearAmplitude = (P: Params, modes: Modes) => nonlinearPoint(P, modes).q;
871769
872770 // Integrator step (s) for motion with energy E: STEP / omega, omega at the
873771 // top of the band, as in the beam. Stretching adds an in-plane tension T
src/plate/style.cssmodified+0−1View file
@@ -5,7 +5,6 @@
55 canvas { display: block; width: 100%; }
66 [hidden] { display: none !important; }
77 #plate { aspect-ratio: 1; cursor: crosshair; }
8-#section { height: 84px; margin-top: 6px; }
98 .column { display: flex; flex-direction: column; gap: 16px; min-width: 0; }
109 #param-form { margin-top: 0; }
1110 .card-head { font-size: 12px; margin-bottom: 12px; display: flex; align-items: center; }