Free edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examples
16 changed files+387−34
README.mdmodified+6−3View file
@@ -17,6 +17,7 @@ npm test # solver tests (Node)
1717 Each bead has a mass and a position in 3D. The script defines:
1818
1919 - *Springs* between pairs of beads, with stiffness k, rest length L0 and optional cubic stiffness k3: the force along the spring is k e + k3 e^3, with e the current length minus L0. A spring whose rest length is shorter than the distance between its beads pulls with a force at rest, and this pre-tension gives a flat net its transverse stiffness, as tension does for a drum head. Because the springs are real springs in 3D, the geometric nonlinearity (a transverse displacement stretches the springs) is exact rather than a truncated expansion.
20+- *Bending triplets* on three beads in a row, i, j, k, with energy kb/2 |x_i - 2 x_j + x_k|^2: the middle bead resists leaving the straight line through its neighbours. The force is linear and the energy is unchanged by translations and rotations. This gives plates and bars their stiffness without tension, so their edges can be free. The rest shape must be straight and evenly spaced along each triplet (flat plates and straight bars; a domed shape will need a different element).
2021 - *Dashpots* between pairs of beads, with force c (v_j - v_i) on bead i. Short-wavelength motion has larger relative velocities between neighbours, so this damping grows with frequency.
2122 - *Friction* on each bead, force -c0 v, which damps all frequencies alike.
2223 - *Pinned* beads, held fixed.
@@ -25,9 +26,9 @@ Each bead has a mass and a position in 3D. The script defines:
2526
2627 Damping is easiest to set by its effect. `[sys.friction, sys.dashpots] = bp_decay(pos, mass, springs, [f1 T1; f2 T2])` chooses them so that every transverse mode of a flat net falls 60 dB in time T(f) following the law sigma(f) = sigma0 + sigma1 (2 pi f)^2 for the amplitude decay rate, through the two given points (one point gives the same decay time at every frequency). It does this by giving each bead friction 2 sigma0 m and each spring a parallel dashpot 2 sigma1 times its tension per length, which makes the damping matrix of the transverse motion 2 sigma0 M + 2 sigma1 K; the modes then stay uncoupled and the law holds exactly at small amplitude. Large motion and the in-plane modes are damped by the same elements, but not by this law. The page plots the small-amplitude decay time of every mode for whatever damping the script sets. A script that sets neither `sys.friction` nor `sys.dashpots` gets a default under which every mode falls 60 dB in 2.5 s; `sys.friction = 0` turns damping off.
2728
28-Bending stiffness (for plates and cymbals) is not in yet.
29+A free object needs support somewhere, or it can drift and spin: pin a small cluster of beads (a single pinned bead lets a plate tilt about it), or hold beads softly with `bp_anchor`, which ties each to four fixed anchors in its own plane with pre-tensioned springs of a chosen out-of-plane stiffness, like a cord or a felt. The run warns when the linear analysis finds modes at 0 Hz.
2930
30-The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads; `bp_chain`: a line of beads), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`), set damping from decay times (`bp_decay`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
31+The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads and bending triplets; `bp_chain`: a line of beads, with its triplets), add soft supports (`bp_anchor`), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`), set damping from decay times (`bp_decay`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
3132
3233 ## Time integration
3334
@@ -58,13 +59,15 @@ The "perturbed twin" option runs the system again with every hit changed by one
5859 | `drum` | A 15 cm drum head (823 beads, fundamental 237 Hz) struck by a soft mallet. Nearly linear at first (2% from the small-strike reference over 10 ms); the slight pitch change at large amplitude accumulates to 25% by 0.3 s. About 7x real time. |
5960 | `cascade` | The same net with a 0.5% pre-stretch, weaker dashpots and a hard strike. The share of energy in modes above 1 kHz grows from 4% (small strike) to about 25% within 30 ms. About 10x real time. |
6061 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
62+| `plate` | A flat plate of 1015 beads, 40 cm across, with a free rim and its center clamped (seven pinned beads), struck sharply near the rim: a flat cymbal or gong. Stiffness from bending triplets only; lowest modes 39 and 58 Hz, highest about 22 kHz. Mildly nonlinear at the given force (it is "thicker" than a real cymbal: nonlinear once the motion reaches about 4 mm). About 20x real time. |
63+| `bar` | A free bar of 41 beads on two soft supports at the nodes of its lowest mode, tuned to 880 Hz (it comes out at 838 Hz). Its overtones are at 2.75, 5.39 and 8.88 times the fundamental, against 2.76, 5.40 and 8.93 for an ideal free bar. |
6164 | `chaos` | The minimal chaotic example: two beads on a short string, free to move in 3D (6 degrees of freedom), struck hard and off-axis, with every mode decaying 60 dB in 2.5 s. With the perturbed twin, the separation e-folds every 2 ms and the runs are unrelated within 0.15 s; at 0.025 of the strike the motion is regular and the separation stays below 1e-8 for 3 s. With half the strike, the damping pulls the motion out of the chaotic range before the runs fully separate (transient chaos). Runs in a fraction of a second. |
6265
6366 Timings are for a laptop (Chrome, single thread).
6467
6568 ## Limitations and next steps
6669
67-- No bending stiffness yet, so no plates or cymbals; the plan is hinge (dihedral angle) elements on pairs of triangles, with a domed rest shape.
70+- Bending works only for straight or flat rest shapes; a domed cymbal needs hinge (dihedral angle) elements on pairs of triangles.
6871 - The solver is plain TypeScript in one Web Worker. Cost scales with the number of springs times the number of steps; nets of several thousand beads will need WASM or WebGPU.
6972 - A bead net has an upper frequency set by its spacing (about 4 kHz for the 1 cm drum example), so a dense high-frequency spectrum needs many beads.
7073 - Linear modes are computed only for flat nets.
matlab/examples/bar.madded+42−0View file
@@ -0,0 +1,42 @@
1+% A bar with free ends, like a glockenspiel or xylophone bar, held softly
2+% (like cords) at the nodes of its lowest mode, 0.224 of the length from
3+% each end. Its stiffness comes from bending triplets, so its overtones are
4+% not harmonic: for an ideal free bar they sit at 2.76 and 5.40 times the
5+% fundamental. The soft supports add low "bouncing" modes of the whole bar
6+% below the fundamental. (Pinning the support beads instead holds them
7+% rigidly, which moves the overtones a long way.)
8+%
9+% The bending stiffness is set from the target pitch: an ideal free bar
10+% has f1 = 22.37/(2*pi*L^2) * sqrt(EI/rhoA), and the triplets give
11+% EI = kb*a^3 and rhoA = m/a.
12+
13+N = 41; % beads
14+L = 0.2; % length (m)
15+rhoA = 0.2; % mass per length (kg/m), e.g. 25 x 3 mm aluminium
16+f1 = 880; % target fundamental (Hz)
17+thick = 0.15; % kb/k: bending relative to stretching stiffness
18+
19+[pos, E, trip] = bp_chain(N, L);
20+a = L / (N - 1);
21+m = rhoA * a;
22+kb = m / a^4 * (2*pi*f1 * L^2 / 22.37)^2;
23+
24+sys.springs = bp_springs(pos, E, kb / thick, 0);
25+supports = [bp_near(pos, [0.224*L 0 0], 0), bp_near(pos, [0.776*L 0 0], 0)];
26+[pos, cords, anchors] = bp_anchor(pos, supports, 500, 0.015); % 500 N/m each
27+sys.pos = pos;
28+sys.mass = m;
29+sys.pinned = anchors;
30+sys.springs = [sys.springs; cords];
31+sys.bends = [trip, kb * ones(size(trip, 1), 1)];
32+[sys.friction, sys.dashpots, sys.bends] = ...
33+ bp_decay(pos, sys.mass, sys.springs, [f1 2.5; 10*f1 0.5], sys.bends);
34+
35+hit.t = 0;
36+hit.beads = bp_near(pos, [L/2 0 0], 0);
37+hit.force = [0 0 -40]; % N
38+hit.duration = 5e-4;
39+sys.hits = hit;
40+
41+sys.listen = struct('beads', N - 2, 'dir', [0 0 1]);
42+sys.duration = 2.5;
matlab/examples/plate.madded+43−0View file
@@ -0,0 +1,43 @@
1+% A flat, free-edged plate held at its center, struck near the rim: a flat
2+% cymbal or gong. Units are SI.
3+%
4+% The plate's stiffness comes from bending triplets (three beads in a row
5+% resist bending, sys.bends), not from tension: the springs carry no load
6+% at rest and only resist stretching. Stretching is what makes the plate
7+% nonlinear, once the motion is comparable to the plate's effective
8+% thickness, about a*sqrt(kb/k) (here 4 mm). The center is clamped by
9+% pinning the middle bead and its six neighbours; with a single pinned bead
10+% the plate could tilt and spin freely.
11+
12+R = 0.2; % radius (m)
13+a = 0.012; % bead spacing (m)
14+sigma = 8.9; % mass per area (kg/m^2), about 1 mm of bronze
15+kb = 7e5; % bending stiffness of each triplet (N/m); sets the pitch
16+thick = 0.1; % kb/k: larger is "thicker", nonlinear only at larger motion
17+
18+[pos, tri, rim, trip] = bp_hexdisk(R, a);
19+E = bp_edges(tri);
20+
21+sys.pos = pos;
22+sys.mass = sigma * (sqrt(3)/2) * a^2;
23+sys.pinned = bp_near(pos, [0 0 0], 1.05 * a); % the mount
24+sys.springs = bp_springs(pos, E, kb / thick, 0); % no tension at rest
25+sys.bends = [trip, kb * ones(size(trip, 1), 1)];
26+[sys.friction, sys.dashpots, sys.bends] = ...
27+ bp_decay(pos, sys.mass, sys.springs, [100 2.5; 5000 1.2], sys.bends);
28+
29+% A stick on the edge: a sharp 0.3 ms push on the beads within 2 cm of the
30+% rim point. A flat plate is only mildly nonlinear at this force (a few
31+% percent from the small-strike reference); at 4x it differs by about 10%
32+% after 0.3 s. A longer, softer push mostly rocks the plate on its mount.
33+hit.t = 0;
34+hit.beads = bp_near(pos, [0.18 0 0], 0.02);
35+hit.force = [0 0 -150]; % peak force per bead (N)
36+hit.duration = 3e-4;
37+sys.hits = hit;
38+
39+sys.listen = struct('beads', bp_near(pos, [-0.1 0.06 0], 0), 'dir', [0 0 1]);
40+sys.duration = 2.5;
41+
42+fprintf('%d beads, %d springs, %d bending triplets, %d struck\n', ...
43+ size(pos, 1), size(E, 1), size(trip, 1), numel(hit.beads));
matlab/helpers/bp_anchor.madded+32−0View file
@@ -0,0 +1,32 @@
1+function [pos, springs, anchors] = bp_anchor(pos, beads, kz, h)
2+% BP_ANCHOR Hold beads softly, like a cord or a felt, rather than pinning them.
3+% [pos, springs, anchors] = bp_anchor(pos, beads, kz, h) adds four fixed
4+% anchor beads around each listed bead, at distance h (default 0.02 m)
5+% along +x, -x, +y and -y in the bead's own plane, and returns
6+% pos the positions with the anchors appended
7+% springs rows [i j k L0] for the new springs (append to sys.springs)
8+% anchors indices of the new beads (append to sys.pinned)
9+% The springs are pre-tensioned so that each held bead feels a stiffness
10+% kz (N/m) against motion out of the plane, and 2*kz against motion in
11+% the plane. The forces balance at rest, and the anchors lie in the
12+% object's plane, so a flat object keeps its linear-mode analysis.
13+if nargin < 4
14+ h = 0.02;
15+end
16+beads = beads(:);
17+nb = numel(beads);
18+n0 = size(pos, 1);
19+dirs = [1 0 0; -1 0 0; 0 1 0; 0 -1 0];
20+f0 = kz * h / 4; % four springs of tension f0 give kz out of plane
21+newpos = zeros(4 * nb, 3);
22+springs = zeros(4 * nb, 4);
23+for q = 1:nb
24+ for d = 1:4
25+ r = 4 * (q - 1) + d;
26+ newpos(r, :) = pos(beads(q), :) + h * dirs(d, :);
27+ springs(r, :) = [beads(q), n0 + r, kz, h - f0 / kz];
28+ end
29+end
30+pos = [pos; newpos];
31+anchors = (n0 + 1:n0 + 4 * nb)';
32+end
matlab/helpers/bp_chain.mmodified+4−2View file
@@ -1,8 +1,10 @@
1-function [pos, E] = bp_chain(n, len)
1+function [pos, E, trip] = bp_chain(n, len)
22 % BP_CHAIN A straight line of n beads along x, from 0 to len.
33 % [pos, E] = bp_chain(n, len) returns positions (n-by-3) and the n-1
4-% neighbor pairs E. The end beads are 1 and n.
4+% neighbor pairs E. The end beads are 1 and n. trip lists the n-2
5+% bending triplets [i-1 i i+1] (for sys.bends).
56 x = linspace(0, len, n)';
67 pos = [x, zeros(n, 2)];
78 E = [(1:n-1)', (2:n)'];
9+trip = [(1:n-2)', (2:n-1)', (3:n)'];
810 end
matlab/helpers/bp_decay.mmodified+18−8View file
@@ -1,4 +1,4 @@
1-function [friction, dashpots, law] = bp_decay(pos, mass, springs, T60)
1+function [friction, dashpots, bends, law] = bp_decay(pos, mass, springs, T60, bends)
22 % BP_DECAY Friction and dashpots that give chosen decay times.
33 % [friction, dashpots] = bp_decay(pos, mass, springs, T60) returns values
44 % for sys.friction and sys.dashpots such that, at small amplitude, each
@@ -9,14 +9,19 @@ function [friction, dashpots, law] = bp_decay(pos, mass, springs, T60)
99 % decay time. T1 must be at least T2 when f1 < f2 (higher modes cannot
1010 % ring longer under this law).
1111 %
12-% Each bead gets friction 2*sigma0*m, and each spring gets a dashpot
13-% 2*sigma1*g in parallel, where g is the spring's tension divided by its
14-% length at rest. The damping is then 2*sigma0*M + 2*sigma1*K for the
15-% transverse motion, so the modes stay uncoupled and the law is exact for
16-% them. Larger motion, and the in-plane modes, are damped by the same
17-% friction and dashpots but not by this law.
12+% [friction, dashpots, bends] = bp_decay(pos, mass, springs, T60, bends)
13+% also damps the bending triplets (rows [i j k kb]): the returned bends
14+% get a fifth column cb, for sys.bends. Pass the bends whenever the
15+% system has them, or plate modes will not follow the law.
1816 %
19-% law = [sigma0 sigma1] is the third output.
17+% Each bead gets friction 2*sigma0*m, each spring a dashpot 2*sigma1*g in
18+% parallel (g is the spring's tension divided by its length at rest), and
19+% each bending triplet the damping cb = 2*sigma1*kb. The damping is then
20+% 2*sigma0*M + 2*sigma1*K for the transverse motion, so the modes stay
21+% uncoupled and the law is exact for them. Larger motion, and the in-plane
22+% modes, are damped by the same elements but not by this law.
23+%
24+% law = [sigma0 sigma1] is the fourth output.
2025 n = size(pos, 1);
2126 s = log(1000) ./ T60(:, 2); % amplitude decay rates (1/s)
2227 w2 = (2 * pi * T60(:, 1)).^2;
@@ -45,5 +50,10 @@ if size(springs, 2) > 4
4550 f = f + springs(:, 5) .* e.^3;
4651 end
4752 dashpots = [i, j, 2 * sigma1 * f ./ L];
53+if nargin > 4
54+ bends = [bends(:, 1:4), 2 * sigma1 * bends(:, 4)];
55+else
56+ bends = zeros(0, 5);
57+end
4858 law = [sigma0, sigma1];
4959 end
matlab/helpers/bp_hexdisk.mmodified+9−1View file
@@ -1,4 +1,4 @@
1-function [pos, tri, rim] = bp_hexdisk(R, a)
1+function [pos, tri, rim, trip] = bp_hexdisk(R, a)
22 % BP_HEXDISK Triangular lattice of beads with spacing a, clipped to a disk of radius R.
33 % [pos, tri, rim] = bp_hexdisk(R, a)
44 % pos n-by-3 bead positions in the plane z = 0, centered at the origin
@@ -6,6 +6,8 @@ function [pos, tri, rim] = bp_hexdisk(R, a)
66 % rim beads with fewer than six lattice neighbors (the edge); pin these
77 % for a drum head. Every other bead then sits at an equilibrium
88 % under uniform spring tension.
9+% trip bending triplets [i j k]: three beads in a row along one of the
10+% three lattice directions, j in the middle (for sys.bends)
911 m = ceil(R / a) + 2;
1012 [I, J] = ndgrid(-2*m:2*m, -m:m);
1113 X = a * (I + J/2);
@@ -23,4 +25,10 @@ tri = tri(all(tri > 0, 2), :);
2325 E = bp_edges(tri);
2426 deg = accumarray(E(:), 1, [size(pos, 1), 1]);
2527 rim = find(deg < 6);
28+% triplets along the three lattice directions: (1,0), (0,1), (1,-1) in (I,J)
29+t1 = [reshape(id(1:end-2, :), [], 1), reshape(id(2:end-1, :), [], 1), reshape(id(3:end, :), [], 1)];
30+t2 = [reshape(id(:, 1:end-2), [], 1), reshape(id(:, 2:end-1), [], 1), reshape(id(:, 3:end), [], 1)];
31+t3 = [reshape(id(1:end-2, 3:end), [], 1), reshape(id(2:end-1, 2:end-1), [], 1), reshape(id(3:end, 1:end-2), [], 1)];
32+trip = [t1; t2; t3];
33+trip = trip(all(trip > 0, 2), :);
2634 end
src/App.tsxmodified+8−4View file
@@ -303,7 +303,7 @@ function Report({ p }: { p: RunPayload }) {
303303 const nPinned = p.sys.pinned.reduce((a, b) => a + b, 0);
304304 const rows: [string, string][] = [
305305 ["beads", `${p.sys.n} (${nPinned} pinned)`],
306- ["springs, dashpots", `${p.sys.si.length}, ${p.sys.nDashpots}`],
306+ ["springs, bends, dashpots", `${p.sys.si.length}, ${p.sys.nBends}, ${p.sys.nDashpots}`],
307307 ["method", p.sys.method],
308308 ["dt", `${r.dt.toExponential(3)} s (${r.substeps} per audio sample)`],
309309 ["highest linear frequency", `${(r.omegaMax / (2 * Math.PI)).toFixed(0)} Hz`],
@@ -382,12 +382,16 @@ const SYS_DOC = `Required
382382 sys.listen struct: beads (index or indices, averaged), dir (default [0 0 1])
383383 Elements
384384 sys.springs rows [i j k L0] or [i j k L0 k3]; force k*e + k3*e^3, e = length - L0
385+ sys.bends rows [i j k kb] or [i j k kb cb], j the middle bead: bending,
386+ energy kb/2 |x_i - 2x_j + x_k|^2 (straight, evenly spaced rest);
387+ cb damps (v_i - 2v_j + v_k). Triplets come from bp_hexdisk, bp_chain
385388 sys.dashpots rows [i j c]; force c*(v_j - v_i) on bead i
386389 sys.friction scalar or n-by-1; force -c0*v
387390 (if neither friction nor dashpots is given: every mode falls 60 dB in 2.5 s)
388- [sys.friction, sys.dashpots] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2])
389- sets both from decay times (T60, s) at two frequencies (Hz)
390- sys.pinned bead indices (or logical mask) held fixed
391+ [sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2], sys.bends)
392+ sets them from decay times (T60, s) at two frequencies (Hz)
393+ sys.pinned bead indices (or logical mask) held fixed; for a soft hold
394+ (cord, felt) use bp_anchor
391395 Excitation
392396 sys.hits struct or struct array, each with
393397 beads, t (default 0), and either
src/sim/analysis.tsmodified+5−2View file
@@ -90,13 +90,16 @@ export function modeBands(
9090 nBands = 48
9191 ): { edges: Float64Array; e: Float32Array; nBands: number } {
9292 const nModes = freqs.length;
93- const fmin = Math.max(freqs[0], 1e-3) * 0.97;
93+ // modes at (nearly) 0 Hz, from free rigid motion, go in the lowest band
94+ const fLow = freqs.find(f => f > 1e-4 * freqs[nModes - 1]) ?? freqs[nModes - 1];
95+ const fmin = fLow * 0.97;
9496 const fmax = freqs[nModes - 1] * 1.03;
9597 const edges = new Float64Array(nBands + 1);
9698 for (let b = 0; b <= nBands; b++) edges[b] = fmin * Math.pow(fmax / fmin, b / nBands);
9799 const band = new Int32Array(nModes);
98100 for (let k = 0; k < nModes; k++) {
99- band[k] = Math.min(nBands - 1, Math.max(0, Math.floor((Math.log(freqs[k] / fmin) / Math.log(fmax / fmin)) * nBands)));
101+ const fk = Math.max(freqs[k], fmin);
102+ band[k] = Math.min(nBands - 1, Math.max(0, Math.floor((Math.log(fk / fmin) / Math.log(fmax / fmin)) * nBands)));
100103 }
101104 const e = new Float32Array(nFrames * nBands);
102105 for (let f = 0; f < nFrames; f++) {
src/sim/modes.tsmodified+13−2View file
@@ -3,8 +3,10 @@
33 // When every bead rests in a plane z = const and every spring lies in that
44 // plane, small transverse (z) motion decouples from in-plane motion, and its
55 // stiffness matrix is a graph Laplacian weighted by spring tension over length:
6-// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T. We return the eigenpairs of
7-// M^{-1/2} K M^{-1/2} over the free beads.
6+// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T, plus, for each bending
7+// triplet, kb (e_i - 2 e_j + e_k)(e_i - 2 e_j + e_k)^T (the triplet energy
8+// separates by component, so its z part decouples too). We return the
9+// eigenpairs of M^{-1/2} K M^{-1/2} over the free beads.
810
911 import type { BeadSystem, Modes } from "./types";
1012
@@ -45,6 +47,14 @@ export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: I
4547 K[b * nf + a] -= g;
4648 }
4749 }
50+ for (let t = 0; t < sys.bi.length; t++) {
51+ const ids = [freeIdx[sys.bi[t]], freeIdx[sys.bj[t]], freeIdx[sys.bk[t]]];
52+ const c = [1, -2, 1];
53+ for (let p = 0; p < 3; p++) {
54+ if (ids[p] < 0) continue;
55+ for (let q = 0; q < 3; q++) if (ids[q] >= 0) K[ids[p] * nf + ids[q]] += sys.bkb[t] * c[p] * c[q];
56+ }
57+ }
4858 for (let a = 0; a < nf; a++) {
4959 for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
5060 }
@@ -82,6 +92,7 @@ export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
8292 let c = 0;
8393 for (let r = 0; r < nf; r++) c += friction[free[r]] * u[free[r]] ** 2;
8494 for (let d = 0; d < di.length; d++) c += dc[d] * (u[di[d]] - u[dj[d]]) ** 2;
95+ for (let t = 0; t < sys.bi.length; t++) c += sys.bcb[t] * (u[sys.bi[t]] - 2 * u[sys.bj[t]] + u[sys.bk[t]]) ** 2;
8596 out[k] = c / 2;
8697 }
8798 return out;
src/sim/parse.tsmodified+29−4View file
@@ -86,7 +86,7 @@ function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeV
8686 }
8787
8888 const KNOWN = new Set([
89- "pos", "mass", "pinned", "springs", "dashpots", "friction", "hits", "listen",
89+ "pos", "mass", "pinned", "springs", "bends", "dashpots", "friction", "hits", "listen",
9090 "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
9191 ]);
9292
@@ -161,8 +161,33 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
161161 dc[d] = col(2);
162162 if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
163163 }
164+ const B = asMat(f.get("bends"), "sys.bends");
165+ const nbend = B && B.data.length ? B.rows : 0;
166+ if (B && nbend && B.cols !== 4 && B.cols !== 5) {
167+ throw new SysError("sys.bends must have 4 or 5 columns: [i j k kb] or [i j k kb cb] (j the middle bead)");
168+ }
169+ const bi = new Int32Array(nbend);
170+ const bj = new Int32Array(nbend);
171+ const bk = new Int32Array(nbend);
172+ const bkb = new Float64Array(nbend);
173+ const bcb = new Float64Array(nbend);
174+ for (let r = 0; r < nbend; r++) {
175+ const col = (c: number) => B!.data[r + c * nbend];
176+ const ids = [col(0), col(1), col(2)];
177+ for (const b of ids) {
178+ if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.bends row ${r + 1}: bead index ${b} is not in 1..${n}`);
179+ }
180+ if (new Set(ids).size !== 3) throw new SysError(`sys.bends row ${r + 1} repeats a bead`);
181+ bi[r] = ids[0] - 1;
182+ bj[r] = ids[1] - 1;
183+ bk[r] = ids[2] - 1;
184+ bkb[r] = col(3);
185+ bcb[r] = B!.cols === 5 ? col(4) : 0;
186+ if (bkb[r] < 0 || bcb[r] < 0) throw new SysError(`sys.bends row ${r + 1}: kb and cb must be nonnegative`);
187+ }
188+
164189 let friction = perBead(f.get("friction"), "sys.friction", n, 0);
165- if (!f.has("friction") && !f.has("dashpots")) {
190+ if (!f.has("friction") && !f.has("dashpots") && !(B && B.cols === 5)) {
166191 // Default damping: every mode falls 60 dB in DEFAULT_T60 seconds
167192 // (amplitude rate sigma = ln(1000)/T60, from friction 2 sigma m per bead).
168193 const sigma = Math.log(1000) / DEFAULT_T60;
@@ -202,7 +227,7 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
202227 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
203228
204229 const sys: BeadSystem = {
205- n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, hits, listenBeads,
230+ n, pos0, mass, pinned, si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, di, dj, dc, friction, hits, listenBeads,
206231 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
207232 duration: scalar(f.get("duration"), "sys.duration", 2),
208233 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
@@ -214,6 +239,6 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
214239 };
215240 if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
216241 if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
217- if (ns === 0) warnings.push("sys.springs is empty: the beads do not interact.");
242+ if (ns === 0 && nbend === 0) warnings.push("sys.springs and sys.bends are empty: the beads do not interact.");
218243 return { sys, warnings };
219244 }
src/sim/solver.tsmodified+78−4View file
@@ -84,8 +84,9 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
8484 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
8585
8686 // ---- forces ------------------------------------------------------------
87- const { si, sj, sk, sL0, sk3 } = sys;
87+ const { si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb } = sys;
8888 const ns = si.length;
89+ const nb = bi.length;
8990 const internalForces = (): number => {
9091 F.fill(0);
9192 let pe = 0;
@@ -109,6 +110,25 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
109110 F[j3 + 1] -= g * dy;
110111 F[j3 + 2] -= g * dz;
111112 }
113+ for (let b = 0; b < nb; b++) {
114+ const i3 = 3 * bi[b];
115+ const j3 = 3 * bj[b];
116+ const k3 = 3 * bk[b];
117+ const kb = bkb[b];
118+ const d0 = x[i3] - 2 * x[j3] + x[k3];
119+ const d1 = x[i3 + 1] - 2 * x[j3 + 1] + x[k3 + 1];
120+ const d2 = x[i3 + 2] - 2 * x[j3 + 2] + x[k3 + 2];
121+ pe += 0.5 * kb * (d0 * d0 + d1 * d1 + d2 * d2);
122+ F[i3] -= kb * d0;
123+ F[i3 + 1] -= kb * d1;
124+ F[i3 + 2] -= kb * d2;
125+ F[j3] += 2 * kb * d0;
126+ F[j3 + 1] += 2 * kb * d1;
127+ F[j3 + 2] += 2 * kb * d2;
128+ F[k3] -= kb * d0;
129+ F[k3 + 1] -= kb * d1;
130+ F[k3 + 2] -= kb * d2;
131+ }
112132 return pe;
113133 };
114134 const pulses = sys.hits.filter(h => h.duration > 0);
@@ -142,7 +162,17 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
142162 const e = L - sL0[s];
143163 maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
144164 }
145- if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
165+ // force scale: the larger of the rest tension and the force of a 1e-3 strain
166+ // on any spring or bending triplet
167+ let forceScale = maxTension;
168+ for (let s = 0; s < ns; s++) forceScale = Math.max(forceScale, 1e-3 * sk[s] * sL0[s]);
169+ for (let b = 0; b < nb; b++) {
170+ const i3 = 3 * bi[b];
171+ const j3 = 3 * bj[b];
172+ const len = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
173+ forceScale = Math.max(forceScale, 1e-3 * bkb[b] * len);
174+ }
175+ if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
146176 warnings.push(
147177 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
148178 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
@@ -170,7 +200,12 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
170200 const cs = stepCoefficients(sys.method);
171201
172202 // ---- damping: v <- exp(-tau A) v ------------------------------------------
173- const { di, dj, dc, friction } = sys;
203+ // dashpots with c = 0 (e.g. from bp_decay on untensioned springs) are dropped
204+ const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0);
205+ const di = Int32Array.from(keepD, d => sys.di[d]);
206+ const dj = Int32Array.from(keepD, d => sys.dj[d]);
207+ const dc = Float64Array.from(keepD, d => sys.dc[d]);
208+ const { friction } = sys;
174209 const nd = di.length;
175210 let normA = 0; // infinity-norm bound of A
176211 {
@@ -180,10 +215,15 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
180215 row[di[d]] += 2 * dc[d];
181216 row[dj[d]] += 2 * dc[d];
182217 }
218+ for (let b = 0; b < nb; b++) {
219+ row[bi[b]] += 4 * bcb[b];
220+ row[bj[b]] += 8 * bcb[b];
221+ row[bk[b]] += 4 * bcb[b];
222+ }
183223 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
184224 }
185225 const hasDamping = normA > 0;
186- const frictionOnly = dc.every(c => c === 0); // then exp(-tau A) is diagonal
226+ const frictionOnly = dc.every(c => c === 0) && bcb.every(c => c === 0); // then exp(-tau A) is diagonal
187227 const term = new Float64Array(N3);
188228 const tmp = new Float64Array(N3);
189229 const fa = new Float64Array(n);
@@ -216,6 +256,28 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
216256 out[j3 + 1] -= aj * r1;
217257 out[j3 + 2] -= aj * r2;
218258 }
259+ for (let b = 0; b < nb; b++) {
260+ const cb = bcb[b];
261+ if (cb === 0) continue;
262+ const i3 = 3 * bi[b];
263+ const j3 = 3 * bj[b];
264+ const k3 = 3 * bk[b];
265+ const ci = cb * invm[bi[b]];
266+ const cj = -2 * cb * invm[bj[b]];
267+ const ck = cb * invm[bk[b]];
268+ const r0 = w[i3] - 2 * w[j3] + w[k3];
269+ const r1 = w[i3 + 1] - 2 * w[j3 + 1] + w[k3 + 1];
270+ const r2 = w[i3 + 2] - 2 * w[j3 + 2] + w[k3 + 2];
271+ out[i3] += ci * r0;
272+ out[i3 + 1] += ci * r1;
273+ out[i3 + 2] += ci * r2;
274+ out[j3] += cj * r0;
275+ out[j3 + 1] += cj * r1;
276+ out[j3 + 2] += cj * r2;
277+ out[k3] += ck * r0;
278+ out[k3 + 1] += ck * r1;
279+ out[k3 + 2] += ck * r2;
280+ }
219281 };
220282 const kinetic = () => {
221283 let ke = 0;
@@ -519,6 +581,18 @@ export function estimateOmegaMax(sys: BeadSystem): number {
519581 y[j3 + 1] -= f1;
520582 y[j3 + 2] -= f2;
521583 }
584+ for (let b = 0; b < sys.bi.length; b++) {
585+ const i3 = 3 * sys.bi[b];
586+ const j3 = 3 * sys.bj[b];
587+ const k3 = 3 * sys.bk[b];
588+ const kb = sys.bkb[b];
589+ for (let q = 0; q < 3; q++) {
590+ const r = kb * (u[i3 + q] - 2 * u[j3 + q] + u[k3 + q]);
591+ y[i3 + q] += r;
592+ y[j3 + q] -= 2 * r;
593+ y[k3 + q] += r;
594+ }
595+ }
522596 let rq = 0;
523597 const wNext = new Float64Array(3 * n);
524598 for (let q = 0; q < 3 * n; q++) {
src/sim/types.tsmodified+7−0View file
@@ -24,6 +24,13 @@ export interface BeadSystem {
2424 sk: Float64Array;
2525 sL0: Float64Array;
2626 sk3: Float64Array;
27+ // Bending triplets (j the middle bead): energy kb/2 |x_i - 2 x_j + x_k|^2,
28+ // with damping force -cb (v_i - 2 v_j + v_k) (1, -2, 1) on (i, j, k)
29+ bi: Int32Array;
30+ bj: Int32Array;
31+ bk: Int32Array;
32+ bkb: Float64Array;
33+ bcb: Float64Array;
2734 // Dashpots on relative velocity: F_i = c (v_j - v_i)
2835 di: Int32Array;
2936 dj: Int32Array;
src/ui/Plots.tsxmodified+4−2View file
@@ -310,9 +310,11 @@ export function DecayTimes({ p }: { p: RunPayload }) {
310310 const ref = useCanvas(
311311 c => {
312312 const f = setupCanvas(c);
313- const T60 = Array.from(sig, s => (s > 0 ? Math.log(1000) / s : Infinity));
313+ const fTop = freqs[freqs.length - 1];
314+ const shown = (k: number) => freqs[k] > 1e-4 * fTop; // skip free rigid motion at 0 Hz
315+ const T60 = Array.from(sig, (s, k) => (s > 0 && shown(k) ? Math.log(1000) / s : Infinity));
314316 const finite = T60.filter(Number.isFinite);
315- const fx0 = Math.log10(freqs[0] * 0.9);
317+ const fx0 = Math.log10((freqs.find(f => f > 1e-4 * fTop) ?? fTop) * 0.9);
316318 const fx1 = Math.log10(freqs[freqs.length - 1] * 1.1);
317319 let ty0 = Math.floor(Math.log10(Math.min(...finite, 1)) * 2) / 2;
318320 let ty1 = Math.ceil(Math.log10(Math.max(...finite, 1)) * 2) / 2;
src/worker.tsmodified+11−1View file
@@ -21,6 +21,7 @@ export interface RunPayload {
2121 si: Int32Array;
2222 sj: Int32Array;
2323 nDashpots: number;
24+ nBends: number;
2425 listenBeads: Int32Array;
2526 hitBeads: Int32Array;
2627 method: string;
@@ -73,6 +74,14 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
7374 else {
7475 post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
7576 modes = computeModes(sys, numblEig);
77+ const fTop = modes.freqs[modes.freqs.length - 1];
78+ const nZero = modes.freqs.filter(f => f < 1e-4 * fTop).length;
79+ if (nZero > 0) {
80+ notes.push(
81+ `${nZero} transverse mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` +
82+ "Pin more beads (a cluster rather than one bead), or tie the object to pinned anchor beads with soft springs."
83+ );
84+ }
7685 }
7786 }
7887
@@ -144,7 +153,8 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
144153 pinned: sys.pinned,
145154 si: sys.si,
146155 sj: sys.sj,
147- nDashpots: sys.di.length,
156+ nDashpots: sys.dc.filter(c => c !== 0).length,
157+ nBends: sys.bi.length,
148158 listenBeads: sys.listenBeads,
149159 hitBeads,
150160 method: sys.method,
tests/solver.test.tsmodified+78−1View file
@@ -206,7 +206,7 @@ for i = 1:N
206206 sys.hits(i).beads = i + 1;
207207 sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
208208 end
209-[sys.friction, sys.dashpots, law] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1]);
209+[sys.friction, sys.dashpots] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1]);
210210 sys.listen.beads = 2:N+1;
211211 sys.listen.dir = [0 0 1];
212212 `).replace("sys.listen = struct('beads', 4);", "")
@@ -258,3 +258,80 @@ describe("default damping", () => {
258258 expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
259259 });
260260 });
261+
262+describe("bending", () => {
263+ // Stiff string: tension plus bending triplets. With theta = k pi/(N+1),
264+ // omega_k^2 = [g (2 - 2 cos theta) + kb (2 - 2 cos theta)^2] / m exactly.
265+ const STIFF = (extra: string) =>
266+ STRING(`
267+[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
268+sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
269+${extra}`);
270+
271+ it("gives the stiff string its exact modes", () => {
272+ const { sys } = build(STIFF(""));
273+ const modes = computeModes(sys, numblEig);
274+ const g = 10 / 0.01;
275+ for (let k = 1; k <= 15; k++) {
276+ const c = 2 - 2 * Math.cos((k * Math.PI) / 16);
277+ const fk = Math.sqrt((g * c + 50 * c * c) / 1e-3) / (2 * Math.PI);
278+ expect(modes.freqs[k - 1] / fk).toBeCloseTo(1, 10);
279+ }
280+ });
281+
282+ it("conserves energy to second order", () => {
283+ const { sys } = build(STIFF("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
284+ const d1 = simulate({ ...sys, dt: 1 / 48000 / 8, duration: 0.05 }).energyDrift;
285+ const d2 = simulate({ ...sys, dt: 1 / 48000 / 16, duration: 0.05 }).energyDrift;
286+ expect(d1).toBeLessThan(1e-3);
287+ expect(d1 / d2).toBeGreaterThan(3);
288+ expect(d1 / d2).toBeLessThan(5);
289+ });
290+
291+ it("follows the bp_decay law with bending damping", () => {
292+ const { sys } = build(STIFF("[sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1], sys.bends);"));
293+ const modes = computeModes(sys, numblEig);
294+ const md = modalDamping(sys, modes);
295+ const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
296+ const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
297+ const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
298+ const sigma0 = s[0] - sigma1 * w2[0];
299+ for (let k = 0; k < modes.freqs.length; k++) {
300+ const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
301+ expect(md[k] / expected).toBeCloseTo(1, 9);
302+ }
303+ // and the solver's damping step agrees: energy balance still exact
304+ const r1 = simulate({ ...sys, hits: [{ ...sys.hits[0], velocity: [0, 0.3, 2] }], dt: 1 / 48000 / 8, duration: 0.05 });
305+ expect(r1.energyDrift).toBeLessThan(1e-3);
306+ });
307+
308+ it("leaves a straight chain in equilibrium after a rigid rotation", () => {
309+ const { sys } = build(
310+ STIFF("R3 = [cos(0.7) -sin(0.7) 0; sin(0.7) cos(0.7) 0; 0 0 1] * [1 0 0; 0 cos(0.4) -sin(0.4); 0 sin(0.4) cos(0.4)]; sys.pos = pos * R3';")
311+ );
312+ const r = simulate({ ...sys, duration: 0.001 });
313+ expect(r.residualForce).toBeLessThan(1e-9);
314+ });
315+});
316+
317+describe("free objects", () => {
318+ it("bp_anchor holds a bead with the chosen out-of-plane stiffness, at equilibrium", () => {
319+ const { sys } = build(`
320+[pos, springs, anchors] = bp_anchor([0 0 0], 1, 400, 0.02);
321+sys.pos = pos; sys.mass = 1e-3; sys.pinned = anchors; sys.springs = springs;
322+hit.beads = 1; hit.velocity = [0 0 1e-4]; sys.hits = hit; sys.listen = struct('beads', 1);
323+sys.friction = 0; sys.duration = 0.01;`);
324+ const modes = computeModes(sys, numblEig);
325+ expect(modes.freqs[0]).toBeCloseTo(Math.sqrt(400 / 1e-3) / (2 * Math.PI), 8);
326+ expect(simulate(sys).residualForce).toBeLessThan(1e-9);
327+ });
328+
329+ it("bar example has the overtones of a free bar", () => {
330+ const { sys } = build(readFileSync("matlab/examples/bar.m", "utf8"));
331+ const f = computeModes(sys, numblEig).freqs;
332+ // two low modes are the bar bouncing on its supports
333+ expect(f[1]).toBeLessThan(100);
334+ const f1 = f[2];
335+ [2.756, 5.404, 8.933].forEach((ratio, q) => expect(Math.abs(f[3 + q] / f1 / ratio - 1)).toBeLessThan(0.01));
336+ });
337+});