Dome: hinge elements (sys.hinges, Bridson gradient, rest angle from sys.pos), bp_hinges, full 3D linear modes for curved nets, axial dashpots so bp_decay is exact for all modes, bp_decay takes and returns sys; dome example
21 changed files+671−163
README.mdmodified+8−6View file
@@ -18,17 +18,18 @@ 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.
2020 - *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).
21-- *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.
21+- *Hinges* across the shared edge j-k of two triangles (i, j, k) and (l, k, j), with energy kh/2 (theta - theta0)^2, where theta is the angle between the triangles and theta0 its value at the rest positions. Any rest shape is therefore free of stress, which is what a dome or a bell needs; `bp_hinges(tri)` lists the hinges of a triangle mesh. The gradient of theta is that of Bridson, Marino and Fedkiw (2003). At rest the stiffness is exactly kh grad(theta) grad(theta)^T, so the time step and the linear modes remain analytic.
22+- *Dashpots* between pairs of beads, with force c (v_j - v_i) on bead i, plus an optional axial part ca (d . (v_j - v_i)) d along the current direction d between the beads. Short-wavelength motion has larger relative velocities between neighbours, so this damping grows with frequency.
2223 - *Friction* on each bead, force -c0 v, which damps all frequencies alike.
2324 - *Pinned* beads, held fixed.
2425 - *Hits*: an instant velocity given to a set of beads, or a raised-cosine force pulse of given peak and duration. Any number, at any times.
2526 - A *listener*: a bead (or the average of several), a direction, and a quantity: the acceleration (the default), velocity or displacement along that direction is the audio signal. Acceleration is the closest of the three to what a small vibrating surface radiates into the air; displacement weights each frequency by 1/f^2 relative to it, so it sounds dull and is dominated by the lowest modes. The acceleration is computed from the forces on the listener bead at each sample (springs, bending, force pulses and damping), not by differencing positions.
2627
27-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.
28+Damping is easiest to set by its effect. `sys = bp_decay(sys, [f1 T1; f2 T2])`, called after the elements are set, chooses the friction, the dashpots and the damping of the bending triplets and hinges 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; each spring a parallel dashpot whose isotropic part is 2 sigma1 times the spring's tension per length g and whose axial part is 2 sigma1 (keff - g) (a spring's stiffness at rest is (keff - g) d d^T + g I); and each triplet and hinge a damping 2 sigma1 times its stiffness. The damping matrix at rest is then 2 sigma0 M + 2 sigma1 K for all the motion, so the modes stay uncoupled and the law holds exactly at small amplitude, for flat and curved nets and for in-plane as well as transverse modes. At large amplitude the same elements damp the motion, but not exactly 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.
2829
2930 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.
3031
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".
32+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`), list hinges of a mesh (`bp_hinges`), 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".
3233
3334 ## Time integration
3435
@@ -44,7 +45,7 @@ Each run reports:
4445
4546 ## Linear modes
4647
47-For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, and its stiffness matrix is a graph Laplacian weighted by spring tension over length. The app computes its eigenpairs with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example) and records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
48+The stiffness at rest is assembled exactly from the elements. For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, so its modes need only the z coordinates; the app computes them with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example). For a curved net (a dome) the modes need all three coordinates of every free bead (up to 1200 coordinates, about 12 s). The app records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
4849
4950 The "small-strike reference" option runs the same system with every hit scaled down by 1000 and the output scaled back up, which is what the system would do if it were linear. The app plays it for comparison and reports the difference from the actual run.
5051
@@ -60,6 +61,7 @@ The "perturbed twin" option runs the system again with every hit changed by one
6061 | `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. |
6162 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
6263 | `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. Bright (spectral centroid of the first 50 ms about 3.2 kHz). Mildly nonlinear at the given force, 1% from the small-strike reference over 10 ms and 14% by 0.3 s; at 4x the force, 14% and 93%. It is "thicker" than a real cymbal: nonlinear once the motion reaches about 4 mm. About 20x real time. |
64+| `dome` | A domed cymbal: 367 beads, 40 cm across, with a 2 cm bow and a 2 cm bell, free at the rim and clamped at the top of the bell, its bending stiffness from hinges. Modes from 63 Hz to about 20 kHz (full 3D analysis, about 12 s). Curvature couples bending to stretching at small amplitude: at the given force the listener signal is 2.2% from the small-strike reference within 10 ms, against 0.3% for the same plate made flat; at 4x the force, the share of energy in modes above 2 kHz rises from 9% to 16% within 10 ms, while the flat plate's stays at 8%. About 20x real time. |
6365 | `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. |
6466 | `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. |
6567
@@ -67,7 +69,7 @@ Timings are for a laptop (Chrome, single thread).
6769
6870 ## Limitations and next steps
6971
70-- Bending works only for straight or flat rest shapes; a domed cymbal needs hinge (dihedral angle) elements on pairs of triangles.
72+- Bending triplets need straight or flat rest shapes; hinges handle curved ones.
7173 - 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.
7274 - 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.
73-- Linear modes are computed only for flat nets.
75+- Linear modes of curved nets need all three coordinates of each free bead and a dense eigensolver, so they are limited to about 400 free beads (1200 coordinates, about 12 s); flat nets go to 1500 free beads.
matlab/examples/bar.mmodified+1−2View file
@@ -29,8 +29,7 @@ sys.mass = m;
2929 sys.pinned = anchors;
3030 sys.springs = [sys.springs; cords];
3131 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);
32+sys = bp_decay(sys, [f1 2.5; 10*f1 0.5]);
3433
3534 hit.t = 0;
3635 hit.beads = bp_near(pos, [L/2 0 0], 0);
matlab/examples/cascade.mmodified+1−1View file
@@ -23,7 +23,7 @@ sys.pos = pos;
2323 sys.mass = sigma * (sqrt(3)/2) * a^2;
2424 sys.pinned = rim;
2525 sys.springs = bp_springs(pos, E, k, f0);
26-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.6]);
26+sys = bp_decay(sys, [200 2.5; 2000 1.6]);
2727
2828 hit.t = 0;
2929 hit.beads = bp_near(pos, [0.05 0 0], 0.03);
matlab/examples/chaos.mmodified+1−1View file
@@ -19,7 +19,7 @@ sys.mass = m;
1919 sys.pinned = [1 4];
2020 sys.springs = bp_springs(pos, E, f0 / (strain * a), f0);
2121 % every mode falls 60 dB in 2.5 s
22-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5]);
22+sys = bp_decay(sys, [200 2.5]);
2323
2424 s = 1; % strike strength; try 0.025 for regular motion
2525 hit.beads = 2;
matlab/examples/dome.madded+50−0View file
@@ -0,0 +1,50 @@
1+% A domed cymbal: a shallow bow with a raised bell, free at the rim and
2+% clamped at the top of the bell, struck sharply near the edge. Units are SI.
3+%
4+% Its bending stiffness comes from hinges across the interior edges of the
5+% triangle mesh (sys.hinges). A hinge resists changes of the angle between
6+% its two triangles from their angle at the rest positions, so the curved
7+% shape below rests without stress. The springs carry no load at rest and
8+% only resist stretching. Curvature couples bending to stretching already at
9+% small amplitude, which a flat plate lacks: set Hbow = 0 and Hbell = 0 to
10+% compare with the same plate flat.
11+%
12+% The linear modes of a curved net need all three coordinates of every free
13+% bead (here about 1100), which takes some seconds; set sys.modes = false to
14+% skip them.
15+
16+R = 0.2; % radius (m)
17+a = 0.02; % bead spacing (m)
18+sigma = 8.9; % mass per area (kg/m^2), about 1 mm of bronze
19+Hbow = 0.02; % height of the bow (m)
20+Hbell = 0.02; % extra height of the bell (m)
21+Rbell = 0.06; % radius of the bell (m)
22+kh = 230; % hinge stiffness (N m per rad^2); sets the pitch
23+k = 8e6; % spring stiffness (N/m); resists stretching
24+
25+[pos, tri] = bp_hexdisk(R, a);
26+r = sqrt(sum(pos(:, 1:2).^2, 2));
27+pos(:, 3) = Hbow * (1 - (r / R).^2) + Hbell * max(0, 1 - (r / Rbell).^2).^2;
28+E = bp_edges(tri);
29+hh = bp_hinges(tri);
30+
31+sys.pos = pos;
32+sys.mass = sigma * (sqrt(3)/2) * a^2;
33+sys.pinned = bp_near(pos, [0 0 max(pos(:, 3))], 1.05 * a); % the mount
34+sys.springs = bp_springs(pos, E, k, 0); % no tension at rest
35+sys.hinges = [hh, kh * ones(size(hh, 1), 1)];
36+sys = bp_decay(sys, [100 2.5; 5000 1.2]);
37+
38+% A stick on the edge: a sharp 0.3 ms push, straight down, on the beads
39+% within 3 cm of a point near the rim.
40+hit.t = 0;
41+hit.beads = bp_near(pos, [0.17 0 0], 0.03);
42+hit.force = [0 0 -150]; % peak force per bead (N)
43+hit.duration = 3e-4;
44+sys.hits = hit;
45+
46+sys.listen = struct('beads', bp_near(pos, [-0.1 0.06 0], 0.01), 'dir', [0 0 1]);
47+sys.duration = 2.5;
48+
49+fprintf('%d beads, %d springs, %d hinges, %d struck\n', ...
50+ size(pos, 1), size(E, 1), size(hh, 1), numel(hit.beads));
matlab/examples/drum.mmodified+2−2View file
@@ -22,9 +22,9 @@ sys.pinned = rim;
2222 sys.springs = bp_springs(pos, E, k, f0);
2323 % Damping, set by decay times (T60: seconds to fall 60 dB) at two
2424 % frequencies: low modes ring about 2.5 s, modes near 2 kHz about 0.75 s.
25-% bp_decay turns this into sys.friction and sys.dashpots; the page plots
25+% bp_decay sets sys.friction and sys.dashpots from them; the page plots
2626 % the decay time each mode gets.
27-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 0.75]);
27+sys = bp_decay(sys, [200 2.5; 2000 0.75]);
2828
2929 % A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
3030 hit.t = 0;
matlab/examples/plate.mmodified+1−2View file
@@ -23,8 +23,7 @@ sys.mass = sigma * (sqrt(3)/2) * a^2;
2323 sys.pinned = bp_near(pos, [0 0 0], 1.05 * a); % the mount
2424 sys.springs = bp_springs(pos, E, kb / thick, 0); % no tension at rest
2525 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);
26+sys = bp_decay(sys, [100 2.5; 5000 1.2]);
2827
2928 % A stick on the edge: a sharp 0.3 ms push on the beads within 2 cm of the
3029 % rim point. The plate is mildly nonlinear at this force (14% from the
matlab/examples/string.mmodified+1−1View file
@@ -19,7 +19,7 @@ sys.pos = pos;
1919 sys.mass = mu * a;
2020 sys.pinned = [1, N + 2];
2121 sys.springs = bp_springs(pos, E, k, T);
22-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.5]); % T60 (s)
22+sys = bp_decay(sys, [200 2.5; 2000 1.5]); % T60 (s)
2323
2424 hit.t = 0;
2525 hit.beads = bp_near(pos, [0.12 0 0], 0.01);
matlab/helpers/bp_decay.mmodified+40−33View file
@@ -1,28 +1,24 @@
1-function [friction, dashpots, bends, law] = bp_decay(pos, mass, springs, T60, bends)
2-% BP_DECAY Friction and dashpots that give chosen decay times.
3-% [friction, dashpots] = bp_decay(pos, mass, springs, T60) returns values
4-% for sys.friction and sys.dashpots such that, at small amplitude, each
5-% transverse mode of a flat net decays at the amplitude rate
1+function sys = bp_decay(sys, T60)
2+% BP_DECAY Set the damping from chosen decay times.
3+% sys = bp_decay(sys, T60) sets sys.friction and sys.dashpots, and adds a
4+% damping column to sys.bends and sys.hinges if present, such that at
5+% small amplitude each mode decays at the amplitude rate
66 % sigma(f) = sigma0 + sigma1 * (2*pi*f)^2,
77 % passing through the decay times in T60 = [f1 T1; f2 T2] (Hz, seconds to
88 % fall by 60 dB). With a single row [f1 T1], every mode gets the same
99 % decay time. T1 must be at least T2 when f1 < f2 (higher modes cannot
10-% ring longer under this law).
10+% ring longer under this law). Call it after the elements are set.
1111 %
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.
12+% Each bead gets friction 2*sigma0*m; each spring a dashpot in parallel
13+% with isotropic part 2*sigma1*g and axial part 2*sigma1*(keff - g), where
14+% g is the spring's tension divided by its length at rest and keff its
15+% stiffness there (a spring's stiffness at rest is (keff - g) d d' + g I);
16+% each bending triplet the damping 2*sigma1*kb, and each hinge 2*sigma1*kh.
17+% The damping matrix at rest is then 2*sigma0*M + 2*sigma1*K for all the
18+% motion, so the modes stay uncoupled and the law holds exactly at small
19+% amplitude, for flat and curved nets alike.
1620 %
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.
25-n = size(pos, 1);
21+% The law [sigma0 sigma1] is stored in sys.decayLaw (not used by the solver).
2622 s = log(1000) ./ T60(:, 2); % amplitude decay rates (1/s)
2723 w2 = (2 * pi * T60(:, 1)).^2;
2824 if size(T60, 1) == 1
@@ -40,20 +36,31 @@ end
4036 if sigma0 < 0
4137 error('bp_decay: decay at the higher frequency is too fast for sigma0 + sigma1*w^2; raise its T60 or lower the other');
4238 end
43-friction = 2 * sigma0 * mass(:) .* ones(n, 1);
44-i = springs(:, 1);
45-j = springs(:, 2);
46-L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
47-e = L - springs(:, 4);
48-f = springs(:, 3) .* e;
49-if size(springs, 2) > 4
50- f = f + springs(:, 5) .* e.^3;
51-end
52-dashpots = [i, j, 2 * sigma1 * f ./ L];
53-if nargin > 4
54- bends = [bends(:, 1:4), 2 * sigma1 * bends(:, 4)];
39+pos = sys.pos;
40+n = size(pos, 1);
41+sys.friction = 2 * sigma0 * sys.mass(:) .* ones(n, 1);
42+if isfield(sys, 'springs') && ~isempty(sys.springs)
43+ S = sys.springs;
44+ i = S(:, 1);
45+ j = S(:, 2);
46+ L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
47+ e = L - S(:, 4);
48+ f = S(:, 3) .* e;
49+ keff = S(:, 3);
50+ if size(S, 2) > 4
51+ f = f + S(:, 5) .* e.^3;
52+ keff = keff + 3 * S(:, 5) .* e.^2;
53+ end
54+ g = f ./ L;
55+ sys.dashpots = [i, j, 2 * sigma1 * g, 2 * sigma1 * (keff - g)];
5556 else
56- bends = zeros(0, 5);
57+ sys.dashpots = zeros(0, 4);
58+end
59+if isfield(sys, 'bends') && ~isempty(sys.bends)
60+ sys.bends = [sys.bends(:, 1:4), 2 * sigma1 * sys.bends(:, 4)];
61+end
62+if isfield(sys, 'hinges') && ~isempty(sys.hinges)
63+ sys.hinges = [sys.hinges(:, 1:5), 2 * sigma1 * sys.hinges(:, 5)];
5764 end
58-law = [sigma0, sigma1];
65+sys.decayLaw = [sigma0, sigma1];
5966 end
matlab/helpers/bp_hinges.madded+15−0View file
@@ -0,0 +1,15 @@
1+function H = bp_hinges(tri)
2+% BP_HINGES Hinges across the interior edges of a triangle mesh.
3+% H = bp_hinges(tri) returns rows [i j k l] for sys.hinges (append a
4+% stiffness column kh): j-k is an edge shared by two triangles, and i and
5+% l are the third beads of those triangles. A hinge resists changes of the
6+% angle between its two triangles from their angle at the rest positions,
7+% so any rest shape (flat, domed, bell-shaped) is free of stress.
8+t = tri;
9+% each triangle edge, sorted, with the opposite bead
10+E3 = [sort(t(:, [1 2]), 2), t(:, 3); sort(t(:, [2 3]), 2), t(:, 1); sort(t(:, [3 1]), 2), t(:, 2)];
11+[~, ord] = sortrows(E3(:, 1:2));
12+E3 = E3(ord, :);
13+pair = find(all(E3(1:end-1, 1:2) == E3(2:end, 1:2), 2));
14+H = [E3(pair, 3), E3(pair, 1), E3(pair, 2), E3(pair + 1, 3)];
15+end
src/App.tsxmodified+10−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, bends, dashpots", `${p.sys.si.length}, ${p.sys.nBends}, ${p.sys.nDashpots}`],
306+ ["springs, bends, hinges, dashpots", `${p.sys.si.length}, ${p.sys.nBends}, ${p.sys.nHinges}, ${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`],
@@ -388,11 +388,17 @@ Elements
388388 sys.bends rows [i j k kb] or [i j k kb cb], j the middle bead: bending,
389389 energy kb/2 |x_i - 2x_j + x_k|^2 (straight, evenly spaced rest);
390390 cb damps (v_i - 2v_j + v_k). Triplets come from bp_hexdisk, bp_chain
391- sys.dashpots rows [i j c]; force c*(v_j - v_i) on bead i
391+ sys.hinges rows [i j k l kh] or [i j k l kh ch]: bending across the edge j-k
392+ shared by triangles (i,j,k) and (l,k,j); energy kh/2 (theta-theta0)^2,
393+ theta0 the angle at sys.pos (so any rest shape, e.g. a dome, is
394+ unstressed). Rows come from bp_hinges(tri)
395+ sys.dashpots rows [i j c] or [i j c ca]; force c*(v_j - v_i) on bead i, plus
396+ ca*(d.(v_j - v_i))*d along the direction d between the beads
392397 sys.friction scalar or n-by-1; force -c0*v
393398 (if neither friction nor dashpots is given: every mode falls 60 dB in 2.5 s)
394- [sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2], sys.bends)
395- sets them from decay times (T60, s) at two frequencies (Hz)
399+ sys = bp_decay(sys, [f1 T1; f2 T2]) sets friction, dashpots and the
400+ damping columns of bends and hinges from decay times (T60, s) at
401+ two frequencies (Hz); call it after the elements are set
396402 sys.pinned bead indices (or logical mask) held fixed; for a soft hold
397403 (cord, felt) use bp_anchor
398404 Excitation
src/sim/hinge.tsadded+67−0View file
@@ -0,0 +1,67 @@
1+// Dihedral angle across a hinge and its gradient.
2+//
3+// Two triangles (i, j, k) and (l, k, j) share the edge j-k; i and l are the
4+// opposite vertices. theta is the signed angle between the two triangle
5+// normals about the edge (0 when the triangles are coplanar). The gradient
6+// is that of Bridson, Marino and Fedkiw (2003), "Simulation of clothing
7+// with folds and wrinkles".
8+
9+// Writes d theta / d(x_i, x_j, x_k, x_l) into g[off .. off+12) and returns theta.
10+export function hingeAngle(x: Float64Array, i: number, j: number, k: number, l: number, g: Float64Array, off = 0): number {
11+ const x1 = 3 * i; // opposite vertex of the first triangle
12+ const x2 = 3 * l; // opposite vertex of the second triangle
13+ const x3 = 3 * j; // edge start
14+ const x4 = 3 * k; // edge end
15+ const e0 = x[x4] - x[x3];
16+ const e1 = x[x4 + 1] - x[x3 + 1];
17+ const e2 = x[x4 + 2] - x[x3 + 2];
18+ const eLen = Math.sqrt(e0 * e0 + e1 * e1 + e2 * e2);
19+ // n1 = (x1 - x3) x (x1 - x4), n2 = (x2 - x4) x (x2 - x3)
20+ const a0 = x[x1] - x[x3], a1 = x[x1 + 1] - x[x3 + 1], a2 = x[x1 + 2] - x[x3 + 2];
21+ const b0 = x[x1] - x[x4], b1 = x[x1 + 1] - x[x4 + 1], b2 = x[x1 + 2] - x[x4 + 2];
22+ const c0 = x[x2] - x[x4], c1 = x[x2 + 1] - x[x4 + 1], c2 = x[x2 + 2] - x[x4 + 2];
23+ const d0 = x[x2] - x[x3], d1 = x[x2 + 1] - x[x3 + 1], d2 = x[x2 + 2] - x[x3 + 2];
24+ const n10 = a1 * b2 - a2 * b1, n11 = a2 * b0 - a0 * b2, n12 = a0 * b1 - a1 * b0;
25+ const n20 = c1 * d2 - c2 * d1, n21 = c2 * d0 - c0 * d2, n22 = c0 * d1 - c1 * d0;
26+ const n1sq = n10 * n10 + n11 * n11 + n12 * n12;
27+ const n2sq = n20 * n20 + n21 * n21 + n22 * n22;
28+ // angle from atan2((n1 x n2) . e/|e|, n1 . n2)
29+ const cr0 = n11 * n22 - n12 * n21, cr1 = n12 * n20 - n10 * n22, cr2 = n10 * n21 - n11 * n20;
30+ const sinPart = (cr0 * e0 + cr1 * e1 + cr2 * e2) / eLen;
31+ const cosPart = n10 * n20 + n11 * n21 + n12 * n22;
32+ // sign chosen to match the gradient below
33+ const theta = -Math.atan2(sinPart, cosPart);
34+ // gradient
35+ const s1 = eLen / n1sq;
36+ const s2 = eLen / n2sq;
37+ // u1 = |e| n1/|n1|^2, u2 = |e| n2/|n2|^2
38+ const p1 = (a0 * e0 + a1 * e1 + a2 * e2) / eLen; // (x1 - x3).e/|e|
39+ const q1 = (b0 * e0 + b1 * e1 + b2 * e2) / eLen; // (x1 - x4).e/|e|
40+ const p2 = (d0 * e0 + d1 * e1 + d2 * e2) / eLen; // (x2 - x3).e/|e|
41+ const q2 = (c0 * e0 + c1 * e1 + c2 * e2) / eLen; // (x2 - x4).e/|e|
42+ const m1 = 1 / n1sq;
43+ const m2 = 1 / n2sq;
44+ const aj = q1 * m1, bj = q2 * m2; // x_j (edge start)
45+ const ak = -p1 * m1, bk = -p2 * m2; // x_k (edge end)
46+ g[off] = s1 * n10; // x_i
47+ g[off + 1] = s1 * n11;
48+ g[off + 2] = s1 * n12;
49+ g[off + 3] = aj * n10 + bj * n20;
50+ g[off + 4] = aj * n11 + bj * n21;
51+ g[off + 5] = aj * n12 + bj * n22;
52+ g[off + 6] = ak * n10 + bk * n20;
53+ g[off + 7] = ak * n11 + bk * n21;
54+ g[off + 8] = ak * n12 + bk * n22;
55+ g[off + 9] = s2 * n20; // x_l
56+ g[off + 10] = s2 * n21;
57+ g[off + 11] = s2 * n22;
58+ return theta;
59+}
60+
61+// theta - theta0 wrapped to (-pi, pi]
62+export function angleDiff(theta: number, theta0: number): number {
63+ let d = theta - theta0;
64+ if (d > Math.PI) d -= 2 * Math.PI;
65+ else if (d <= -Math.PI) d += 2 * Math.PI;
66+ return d;
67+}
src/sim/modes.tsmodified+123−56View file
@@ -1,13 +1,20 @@
1-// Linear transverse modes of a flat bead system.
1+// Linear modes of a bead system about its rest positions.
22 //
3-// When every bead rests in a plane z = const and every spring lies in that
4-// plane, small transverse (z) motion decouples from in-plane motion, and its
5-// 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, 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.
3+// The stiffness at rest is assembled from the elements: a spring contributes
4+// the 3x3 block keff d d^T + (f/L)(I - d d^T) (d its direction, f its tension,
5+// keff = k + 3 k3 e^2); a bending triplet kb (e_i - 2 e_j + e_k)(...)^T in each
6+// component; a hinge kh grad(theta) grad(theta)^T (exact at rest, where
7+// theta = theta0). We return eigenpairs of M^{-1/2} K M^{-1/2} over a set of
8+// degrees of freedom (dof: indices 3*bead + component).
9+//
10+// When every bead rests in a plane z = const, small z motion decouples from
11+// in-plane motion (springs in the plane have no z-to-x/y coupling, triplets
12+// separate by component, and flat hinges respond only to z), so the
13+// "transverse" modes use z alone. Otherwise all three components are needed
14+// ("full" modes), which is about 27 times the work of a transverse problem
15+// with the same number of beads.
1016
17+import { hingeAngle } from "./hinge";
1118 import type { BeadSystem, Modes } from "./types";
1219
1320 export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
@@ -25,74 +32,134 @@ export function isPlanar(sys: BeadSystem): boolean {
2532 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
2633 }
2734
28-export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
29- const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
30- const freeIdx = new Int32Array(n).fill(-1);
31- const free: number[] = [];
32- for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
33- const nf = free.length;
34- const K = new Float64Array(nf * nf);
35- for (let s = 0; s < si.length; s++) {
36- const i3 = 3 * si[s];
37- const j3 = 3 * sj[s];
38- const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
39- const e = L - sL0[s];
40- const g = (sk[s] * e + sk3[s] * e * e * e) / L;
41- const a = freeIdx[si[s]];
42- const b = freeIdx[sj[s]];
43- if (a >= 0) K[a * nf + a] += g;
44- if (b >= 0) K[b * nf + b] += g;
45- if (a >= 0 && b >= 0) {
46- K[a * nf + b] -= g;
47- K[b * nf + a] -= g;
35+export function freeBeads(sys: BeadSystem): number[] {
36+ const out: number[] = [];
37+ for (let i = 0; i < sys.n; i++) if (!sys.pinned[i]) out.push(i);
38+ return out;
39+}
40+
41+export function modeDofs(sys: BeadSystem, kind: Modes["kind"]): Int32Array {
42+ const free = freeBeads(sys);
43+ return kind === "transverse"
44+ ? Int32Array.from(free, b => 3 * b + 2)
45+ : Int32Array.from(free.flatMap(b => [3 * b, 3 * b + 1, 3 * b + 2]));
46+}
47+
48+// M^{-1/2} K M^{-1/2} restricted to dof, dense, row-major.
49+export function stiffnessMatrix(sys: BeadSystem, dof: Int32Array): Float64Array {
50+ const { n, pos0, mass } = sys;
51+ const nd = dof.length;
52+ const index = new Int32Array(3 * n).fill(-1);
53+ dof.forEach((q, r) => (index[q] = r));
54+ const K = new Float64Array(nd * nd);
55+ const add = (qa: number, qb: number, v: number) => {
56+ const a = index[qa];
57+ const b = index[qb];
58+ if (a >= 0 && b >= 0) K[a * nd + b] += v;
59+ };
60+ for (let s = 0; s < sys.si.length; s++) {
61+ const i = sys.si[s];
62+ const j = sys.sj[s];
63+ const dx = [0, 1, 2].map(c => pos0[3 * j + c] - pos0[3 * i + c]);
64+ const L = Math.hypot(dx[0], dx[1], dx[2]);
65+ const e = L - sys.sL0[s];
66+ const keff = sys.sk[s] + 3 * sys.sk3[s] * e * e;
67+ const g = (sys.sk[s] * e + sys.sk3[s] * e * e * e) / L;
68+ const u = dx.map(d => d / L);
69+ for (let a = 0; a < 3; a++) {
70+ for (let b = 0; b < 3; b++) {
71+ const B = (keff - g) * u[a] * u[b] + (a === b ? g : 0);
72+ if (B === 0) continue;
73+ add(3 * i + a, 3 * i + b, B);
74+ add(3 * j + a, 3 * j + b, B);
75+ add(3 * i + a, 3 * j + b, -B);
76+ add(3 * j + a, 3 * i + b, -B);
77+ }
4878 }
4979 }
80+ const cb = [1, -2, 1];
5081 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];
82+ const beads = [sys.bi[t], sys.bj[t], sys.bk[t]];
5383 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];
84+ for (let q = 0; q < 3; q++) {
85+ for (let c = 0; c < 3; c++) add(3 * beads[p] + c, 3 * beads[q] + c, sys.bkb[t] * cb[p] * cb[q]);
86+ }
5687 }
5788 }
58- for (let a = 0; a < nf; a++) {
59- for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
89+ const g = new Float64Array(12);
90+ for (let h = 0; h < sys.hi.length; h++) {
91+ const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
92+ hingeAngle(pos0, beads[0], beads[1], beads[2], beads[3], g);
93+ for (let p = 0; p < 12; p++) {
94+ if (g[p] === 0) continue;
95+ for (let q = 0; q < 12; q++) {
96+ add(3 * beads[(p / 3) | 0] + (p % 3), 3 * beads[(q / 3) | 0] + (q % 3), sys.hkh[h] * g[p] * g[q]);
97+ }
98+ }
99+ }
100+ for (let a = 0; a < nd; a++) {
101+ for (let b = 0; b < nd; b++) K[a * nd + b] /= Math.sqrt(mass[(dof[a] / 3) | 0] * mass[(dof[b] / 3) | 0]);
60102 }
61- return { K, free: Int32Array.from(free) };
103+ return K;
62104 }
63105
64-export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
65- const { K, free } = transverseStiffness(sys);
66- const nf = free.length;
67- const { values, vectors } = eig(K, nf);
106+export function computeModes(sys: BeadSystem, eig: SymEig, kind: Modes["kind"] = isPlanar(sys) ? "transverse" : "full"): Modes {
107+ const dof = modeDofs(sys, kind);
108+ const nd = dof.length;
109+ const { values, vectors } = eig(stiffnessMatrix(sys, dof), nd);
68110 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
69- const freqs = new Float64Array(nf);
70- const psi = new Float64Array(nf * nf);
111+ const freqs = new Float64Array(nd);
112+ const psi = new Float64Array(nd * nd);
71113 order.forEach((src, k) => {
72114 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
73- psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
115+ psi.set(vectors.subarray(src * nd, (src + 1) * nd), k * nd);
74116 });
75- return { freqs, psi, free };
117+ return { freqs, psi, dof, kind };
76118 }
77119
78-// Small-amplitude decay rate of each transverse mode (amplitude, 1/s), from
79-// the friction and dashpots: sigma_k = (1/2) u^T C u with u the mode shape
80-// normalized so u^T M u = 1. Exact when the damping matrix is a combination
81-// of M and K (as bp_decay builds it); otherwise it neglects the coupling of
82-// modes by damping, which is small for light damping.
120+// Small-amplitude decay rate of each mode (amplitude, 1/s): sigma_k =
121+// (1/2) u^T C u with u the mode shape normalized so u^T M u = 1 and C the
122+// damping matrix at rest (friction, dashpots, bending and hinge damping).
123+// Exact when C is a combination of M and K (as bp_decay builds it for flat
124+// nets); otherwise it neglects the coupling of modes by damping, which is
125+// small for light damping.
83126 export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
84- const { n, mass, friction, di, dj, dc } = sys;
85- const { free, psi, freqs } = modes;
86- const nf = free.length;
87- const u = new Float64Array(n);
127+ const { n, mass, friction, di, dj, dc, pos0 } = sys;
128+ const { dof, psi, freqs } = modes;
129+ const nd = dof.length;
130+ const u = new Float64Array(3 * n);
88131 const out = new Float64Array(freqs.length);
132+ const hgs = new Float64Array(12 * sys.hi.length);
133+ for (let h = 0; h < sys.hi.length; h++) {
134+ hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], hgs.subarray(12 * h, 12 * h + 12));
135+ }
89136 for (let k = 0; k < freqs.length; k++) {
90137 u.fill(0);
91- for (let r = 0; r < nf; r++) u[free[r]] = psi[k * nf + r] / Math.sqrt(mass[free[r]]);
138+ for (let r = 0; r < nd; r++) u[dof[r]] = psi[k * nd + r] / Math.sqrt(mass[(dof[r] / 3) | 0]);
92139 let c = 0;
93- for (let r = 0; r < nf; r++) c += friction[free[r]] * u[free[r]] ** 2;
94- 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;
140+ for (let i = 0; i < n; i++) c += friction[i] * (u[3 * i] ** 2 + u[3 * i + 1] ** 2 + u[3 * i + 2] ** 2);
141+ for (let d = 0; d < di.length; d++) {
142+ let ax = 0;
143+ let len = 0;
144+ for (let q = 0; q < 3; q++) {
145+ const r = u[3 * di[d] + q] - u[3 * dj[d] + q];
146+ c += dc[d] * r * r;
147+ const e = pos0[3 * dj[d] + q] - pos0[3 * di[d] + q];
148+ ax += e * r;
149+ len += e * e;
150+ }
151+ if (sys.dca[d] !== 0) c += (sys.dca[d] * ax * ax) / len;
152+ }
153+ for (let t = 0; t < sys.bi.length; t++) {
154+ for (let q = 0; q < 3; q++) c += sys.bcb[t] * (u[3 * sys.bi[t] + q] - 2 * u[3 * sys.bj[t] + q] + u[3 * sys.bk[t] + q]) ** 2;
155+ }
156+ for (let h = 0; h < sys.hi.length; h++) {
157+ if (sys.hch[h] === 0) continue;
158+ const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
159+ let s = 0;
160+ for (let p = 0; p < 4; p++) for (let q = 0; q < 3; q++) s += hgs[12 * h + 3 * p + q] * u[3 * beads[p] + q];
161+ c += sys.hch[h] * s * s;
162+ }
96163 out[k] = c / 2;
97164 }
98165 return out;
src/sim/parse.tsmodified+43−6View file
@@ -2,6 +2,7 @@
22 // The script uses 1-based bead indices; they become 0-based here.
33
44 import type { RuntimeValue } from "numbl";
5+import { hingeAngle } from "./hinge";
56 import type { BeadSystem, Hit, ListenQuantity, Method } from "./types";
67
78 type Mat = { rows: number; cols: number; data: Float64Array };
@@ -86,8 +87,8 @@ function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeV
8687 }
8788
8889 const KNOWN = new Set([
89- "pos", "mass", "pinned", "springs", "bends", "dashpots", "friction", "hits", "listen",
90- "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
90+ "pos", "mass", "pinned", "springs", "bends", "hinges", "dashpots", "friction", "hits", "listen",
91+ "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes", "decayLaw",
9192 ]);
9293
9394 export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
@@ -147,10 +148,11 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
147148
148149 const D = asMat(f.get("dashpots"), "sys.dashpots");
149150 const nd = D && D.data.length ? D.rows : 0;
150- if (D && nd && D.cols !== 3) throw new SysError("sys.dashpots must have 3 columns: [i j c]");
151+ if (D && nd && D.cols !== 3 && D.cols !== 4) throw new SysError("sys.dashpots must have 3 or 4 columns: [i j c] or [i j c ca]");
151152 const di = new Int32Array(nd);
152153 const dj = new Int32Array(nd);
153154 const dc = new Float64Array(nd);
155+ const dca = new Float64Array(nd);
154156 for (let d = 0; d < nd; d++) {
155157 const col = (c: number) => D!.data[d + c * nd];
156158 for (const b of [col(0), col(1)]) {
@@ -159,6 +161,8 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
159161 di[d] = col(0) - 1;
160162 dj[d] = col(1) - 1;
161163 dc[d] = col(2);
164+ dca[d] = D!.cols === 4 ? col(3) : 0;
165+ if (dca[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: ca must be nonnegative`);
162166 if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
163167 }
164168 const B = asMat(f.get("bends"), "sys.bends");
@@ -186,8 +190,41 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
186190 if (bkb[r] < 0 || bcb[r] < 0) throw new SysError(`sys.bends row ${r + 1}: kb and cb must be nonnegative`);
187191 }
188192
193+ const Hm = asMat(f.get("hinges"), "sys.hinges");
194+ const nh = Hm && Hm.data.length ? Hm.rows : 0;
195+ if (Hm && nh && Hm.cols !== 5 && Hm.cols !== 6) {
196+ throw new SysError("sys.hinges must have 5 or 6 columns: [i j k l kh] or [i j k l kh ch] (j-k the shared edge, i and l the opposite beads)");
197+ }
198+ const hi = new Int32Array(nh);
199+ const hj = new Int32Array(nh);
200+ const hk = new Int32Array(nh);
201+ const hl = new Int32Array(nh);
202+ const hkh = new Float64Array(nh);
203+ const hch = new Float64Array(nh);
204+ const hth0 = new Float64Array(nh);
205+ const gtmp = new Float64Array(12);
206+ for (let r = 0; r < nh; r++) {
207+ const col = (c: number) => Hm!.data[r + c * nh];
208+ const ids = [col(0), col(1), col(2), col(3)];
209+ for (const b of ids) {
210+ if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.hinges row ${r + 1}: bead index ${b} is not in 1..${n}`);
211+ }
212+ if (new Set(ids).size !== 4) throw new SysError(`sys.hinges row ${r + 1} repeats a bead`);
213+ hi[r] = ids[0] - 1;
214+ hj[r] = ids[1] - 1;
215+ hk[r] = ids[2] - 1;
216+ hl[r] = ids[3] - 1;
217+ hkh[r] = col(4);
218+ hch[r] = Hm!.cols === 6 ? col(5) : 0;
219+ if (hkh[r] < 0 || hch[r] < 0) throw new SysError(`sys.hinges row ${r + 1}: kh and ch must be nonnegative`);
220+ hth0[r] = hingeAngle(pos0, hi[r], hj[r], hk[r], hl[r], gtmp);
221+ if (!Number.isFinite(hth0[r]) || !gtmp.every(Number.isFinite)) {
222+ throw new SysError(`sys.hinges row ${r + 1}: a triangle has zero area at the rest positions`);
223+ }
224+ }
225+
189226 let friction = perBead(f.get("friction"), "sys.friction", n, 0);
190- if (!f.has("friction") && !f.has("dashpots") && !(B && B.cols === 5)) {
227+ if (!f.has("friction") && !f.has("dashpots") && !(B && B.cols === 5) && !(Hm && Hm.cols === 6)) {
191228 // Default damping: every mode falls 60 dB in DEFAULT_T60 seconds
192229 // (amplitude rate sigma = ln(1000)/T60, from friction 2 sigma m per bead).
193230 const sigma = Math.log(1000) / DEFAULT_T60;
@@ -231,7 +268,7 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
231268 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
232269
233270 const sys: BeadSystem = {
234- n, pos0, mass, pinned, si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, di, dj, dc, friction, hits, listenBeads,
271+ n, pos0, mass, pinned, si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, hi, hj, hk, hl, hkh, hch, hth0, di, dj, dc, dca, friction, hits, listenBeads,
235272 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
236273 listenQuantity,
237274 duration: scalar(f.get("duration"), "sys.duration", 2),
@@ -244,6 +281,6 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
244281 };
245282 if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
246283 if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
247- if (ns === 0 && nbend === 0) warnings.push("sys.springs and sys.bends are empty: the beads do not interact.");
284+ if (ns === 0 && nbend === 0 && nh === 0) warnings.push("sys.springs, sys.bends and sys.hinges are empty: the beads do not interact.");
248285 return { sys, warnings };
249286 }
src/sim/solver.tsmodified+160−28View file
@@ -9,6 +9,7 @@
99 // which is small compared with the oscillation frequencies; the conservative
1010 // part is second order (Verlet) or fourth order (Yoshida) in the step h.
1111
12+import { angleDiff, hingeAngle } from "./hinge";
1213 import type { BeadSystem, Method, Modes, SimResult } from "./types";
1314
1415 const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
@@ -84,7 +85,10 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
8485 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
8586
8687 // ---- forces ------------------------------------------------------------
87- const { si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb } = sys;
88+ const { si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, hi, hj, hk, hl, hkh, hch, hth0 } = sys;
89+ const nhg = hi.length;
90+ const hg = new Float64Array(12 * nhg); // hinge gradients, refreshed where needed
91+ const hb = (h: number) => [hi[h], hj[h], hk[h], hl[h]]; // the hinge's beads, in gradient order
8892 const ns = si.length;
8993 const nb = bi.length;
9094 const internalForces = (): number => {
@@ -129,6 +133,25 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
129133 F[k3 + 1] -= kb * d1;
130134 F[k3 + 2] -= kb * d2;
131135 }
136+ for (let h = 0; h < nhg; h++) {
137+ const o = 12 * h;
138+ const d = angleDiff(hingeAngle(x, hi[h], hj[h], hk[h], hl[h], hg, o), hth0[h]);
139+ pe += 0.5 * hkh[h] * d * d;
140+ const s = hkh[h] * d;
141+ const i3 = 3 * hi[h], j3 = 3 * hj[h], k3 = 3 * hk[h], l3 = 3 * hl[h];
142+ F[i3] -= s * hg[o];
143+ F[i3 + 1] -= s * hg[o + 1];
144+ F[i3 + 2] -= s * hg[o + 2];
145+ F[j3] -= s * hg[o + 3];
146+ F[j3 + 1] -= s * hg[o + 4];
147+ F[j3 + 2] -= s * hg[o + 5];
148+ F[k3] -= s * hg[o + 6];
149+ F[k3 + 1] -= s * hg[o + 7];
150+ F[k3 + 2] -= s * hg[o + 8];
151+ F[l3] -= s * hg[o + 9];
152+ F[l3 + 1] -= s * hg[o + 10];
153+ F[l3 + 2] -= s * hg[o + 11];
154+ }
132155 return pe;
133156 };
134157 const pulses = sys.hits.filter(h => h.duration > 0);
@@ -172,6 +195,13 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
172195 const len = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
173196 forceScale = Math.max(forceScale, 1e-3 * bkb[b] * len);
174197 }
198+ {
199+ const g = new Float64Array(12);
200+ for (let h = 0; h < nhg; h++) {
201+ hingeAngle(x, hi[h], hj[h], hk[h], hl[h], g);
202+ forceScale = Math.max(forceScale, 1e-3 * hkh[h] * Math.max(...g.map(Math.abs)));
203+ }
204+ }
175205 if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
176206 warnings.push(
177207 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
@@ -200,11 +230,14 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
200230 const cs = stepCoefficients(sys.method);
201231
202232 // ---- damping: v <- exp(-tau A) v ------------------------------------------
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);
233+ // dashpots with c = ca = 0 are dropped
234+ const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0 || sys.dca[d] !== 0);
205235 const di = Int32Array.from(keepD, d => sys.di[d]);
206236 const dj = Int32Array.from(keepD, d => sys.dj[d]);
207237 const dc = Float64Array.from(keepD, d => sys.dc[d]);
238+ const dca = Float64Array.from(keepD, d => sys.dca[d]);
239+ const axialDamped = dca.some(c => c > 0);
240+ const dd = new Float64Array(3 * dca.length); // dashpot directions, refreshed where needed
208241 const { friction } = sys;
209242 const nd = di.length;
210243 let normA = 0; // infinity-norm bound of A
@@ -214,6 +247,8 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
214247 for (let d = 0; d < nd; d++) {
215248 row[di[d]] += 2 * dc[d];
216249 row[dj[d]] += 2 * dc[d];
250+ row[di[d]] += 2 * dca[d];
251+ row[dj[d]] += 2 * dca[d];
217252 }
218253 for (let b = 0; b < nb; b++) {
219254 row[bi[b]] += 4 * bcb[b];
@@ -222,18 +257,49 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
222257 }
223258 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
224259 }
225- const hasDamping = normA > 0;
226- const frictionOnly = dc.every(c => c === 0) && bcb.every(c => c === 0); // then exp(-tau A) is diagonal
260+ const hingeDamped = hch.some(c => c > 0);
261+ const hasDamping = normA > 0 || hingeDamped;
262+ const frictionOnly = dc.every(c => c === 0) && !axialDamped && bcb.every(c => c === 0) && !hingeDamped; // then exp(-tau A) is diagonal
263+ // Hinge damping depends on the configuration through grad theta, so its
264+ // gradients (and its part of the norm bound) are refreshed at the current
265+ // positions before each damping step; positions do not change during it.
266+ let normAh = 0;
267+ const hingeRow = new Float64Array(n);
268+ const refreshHinges = () => {
269+ if (axialDamped) {
270+ for (let d = 0; d < dca.length; d++) {
271+ const i3 = 3 * di[d];
272+ const j3 = 3 * dj[d];
273+ const e0 = x[j3] - x[i3];
274+ const e1 = x[j3 + 1] - x[i3 + 1];
275+ const e2 = x[j3 + 2] - x[i3 + 2];
276+ const L = Math.sqrt(e0 * e0 + e1 * e1 + e2 * e2) || 1;
277+ dd[3 * d] = e0 / L;
278+ dd[3 * d + 1] = e1 / L;
279+ dd[3 * d + 2] = e2 / L;
280+ }
281+ }
282+ if (!hingeDamped) return;
283+ hingeRow.fill(0);
284+ const row = hingeRow;
285+ for (let h = 0; h < nhg; h++) {
286+ const o = 12 * h;
287+ hingeAngle(x, hi[h], hj[h], hk[h], hl[h], hg, o);
288+ let total = 0;
289+ for (let q = 0; q < 12; q++) total += Math.abs(hg[o + q]);
290+ const c = hch[h] * total;
291+ row[hi[h]] += c * (Math.abs(hg[o]) + Math.abs(hg[o + 1]) + Math.abs(hg[o + 2]));
292+ row[hj[h]] += c * (Math.abs(hg[o + 3]) + Math.abs(hg[o + 4]) + Math.abs(hg[o + 5]));
293+ row[hk[h]] += c * (Math.abs(hg[o + 6]) + Math.abs(hg[o + 7]) + Math.abs(hg[o + 8]));
294+ row[hl[h]] += c * (Math.abs(hg[o + 9]) + Math.abs(hg[o + 10]) + Math.abs(hg[o + 11]));
295+ }
296+ normAh = 0;
297+ for (let i = 0; i < n; i++) normAh = Math.max(normAh, row[i] * invm[i]);
298+ };
227299 const term = new Float64Array(N3);
228300 const tmp = new Float64Array(N3);
229301 const fa = new Float64Array(n);
230302 for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
231- const dai = new Float64Array(nd);
232- const daj = new Float64Array(nd);
233- for (let d = 0; d < nd; d++) {
234- dai[d] = dc[d] * invm[di[d]];
235- daj[d] = dc[d] * invm[dj[d]];
236- }
237303 const applyA = (w: Float64Array, out: Float64Array) => {
238304 for (let i = 0; i < n; i++) {
239305 const a = fa[i];
@@ -244,17 +310,28 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
244310 for (let d = 0; d < nd; d++) {
245311 const i3 = 3 * di[d];
246312 const j3 = 3 * dj[d];
247- const ai = dai[d];
248- const aj = daj[d];
249313 const r0 = w[i3] - w[j3];
250314 const r1 = w[i3 + 1] - w[j3 + 1];
251315 const r2 = w[i3 + 2] - w[j3 + 2];
252- out[i3] += ai * r0;
253- out[i3 + 1] += ai * r1;
254- out[i3 + 2] += ai * r2;
255- out[j3] -= aj * r0;
256- out[j3 + 1] -= aj * r1;
257- out[j3 + 2] -= aj * r2;
316+ // force on i: -c r - ca (d . r) d, with r = w_i - w_j
317+ let f0 = dc[d] * r0;
318+ let f1 = dc[d] * r1;
319+ let f2 = dc[d] * r2;
320+ if (dca[d] !== 0) {
321+ const u0 = dd[3 * d], u1 = dd[3 * d + 1], u2 = dd[3 * d + 2];
322+ const s = dca[d] * (u0 * r0 + u1 * r1 + u2 * r2);
323+ f0 += s * u0;
324+ f1 += s * u1;
325+ f2 += s * u2;
326+ }
327+ const mi = invm[di[d]];
328+ const mj = invm[dj[d]];
329+ out[i3] += mi * f0;
330+ out[i3 + 1] += mi * f1;
331+ out[i3 + 2] += mi * f2;
332+ out[j3] -= mj * f0;
333+ out[j3 + 1] -= mj * f1;
334+ out[j3 + 2] -= mj * f2;
258335 }
259336 for (let b = 0; b < nb; b++) {
260337 const cb = bcb[b];
@@ -278,6 +355,33 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
278355 out[k3 + 1] += ck * r1;
279356 out[k3 + 2] += ck * r2;
280357 }
358+ if (hingeDamped) {
359+ for (let h = 0; h < nhg; h++) {
360+ const ch = hch[h];
361+ if (ch === 0) continue;
362+ const o = 12 * h;
363+ const i3 = 3 * hi[h], j3 = 3 * hj[h], k3 = 3 * hk[h], l3 = 3 * hl[h];
364+ let s =
365+ hg[o] * w[i3] + hg[o + 1] * w[i3 + 1] + hg[o + 2] * w[i3 + 2] +
366+ hg[o + 3] * w[j3] + hg[o + 4] * w[j3 + 1] + hg[o + 5] * w[j3 + 2] +
367+ hg[o + 6] * w[k3] + hg[o + 7] * w[k3 + 1] + hg[o + 8] * w[k3 + 2] +
368+ hg[o + 9] * w[l3] + hg[o + 10] * w[l3 + 1] + hg[o + 11] * w[l3 + 2];
369+ s *= ch;
370+ const ai = s * invm[hi[h]], aj = s * invm[hj[h]], ak = s * invm[hk[h]], al = s * invm[hl[h]];
371+ out[i3] += ai * hg[o];
372+ out[i3 + 1] += ai * hg[o + 1];
373+ out[i3 + 2] += ai * hg[o + 2];
374+ out[j3] += aj * hg[o + 3];
375+ out[j3 + 1] += aj * hg[o + 4];
376+ out[j3 + 2] += aj * hg[o + 5];
377+ out[k3] += ak * hg[o + 6];
378+ out[k3 + 1] += ak * hg[o + 7];
379+ out[k3 + 2] += ak * hg[o + 8];
380+ out[l3] += al * hg[o + 9];
381+ out[l3 + 1] += al * hg[o + 10];
382+ out[l3 + 2] += al * hg[o + 11];
383+ }
384+ }
281385 };
282386 const kinetic = () => {
283387 let ke = 0;
@@ -299,7 +403,8 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
299403 }
300404 } else {
301405 // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
302- const nsub = Math.max(1, Math.ceil((Math.abs(tau) * normA) / 0.5));
406+ refreshHinges();
407+ const nsub = Math.max(1, Math.ceil((Math.abs(tau) * (normA + normAh)) / 0.5));
303408 const ts = tau / nsub;
304409 for (let r = 0; r < nsub; r++) {
305410 let vmax = 0;
@@ -337,17 +442,32 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
337442 // the dashpots and bending triplets that touch each listener bead.
338443 const touchD = Array.from(lb, b => Array.from(di.keys()).filter(d => di[d] === b || dj[d] === b));
339444 const touchB = Array.from(lb, b => Array.from(bi.keys()).filter(t => (bcb[t] !== 0 && (bi[t] === b || bj[t] === b || bk[t] === b))));
445+ const touchH = Array.from(lb, b => Array.from(hi.keys()).filter(h => hch[h] !== 0 && hb(h).includes(b)));
340446 const dampForce = (v: Float64Array, q: number, b: number, c: number) => {
341447 let f = -friction[b] * v[3 * b + c];
342448 for (const d of touchD[q]) {
343449 const other = di[d] === b ? dj[d] : di[d];
344450 f -= dc[d] * (v[3 * b + c] - v[3 * other + c]);
451+ if (dca[d] !== 0) {
452+ const rel = [0, 1, 2].reduce((acc, q) => acc + dd[3 * d + q] * (v[3 * b + q] - v[3 * other + q]), 0);
453+ f -= dca[d] * rel * dd[3 * d + c];
454+ }
345455 }
346456 for (const t of touchB[q]) {
347457 const r = v[3 * bi[t] + c] - 2 * v[3 * bj[t] + c] + v[3 * bk[t] + c];
348458 const w = bj[t] === b ? -2 : 1;
349459 f -= bcb[t] * w * r;
350460 }
461+ for (const h of touchH[q]) {
462+ const o = 12 * h;
463+ const beads = hb(h);
464+ let s = 0;
465+ for (let p = 0; p < 4; p++) {
466+ const b3 = 3 * beads[p];
467+ s += hg[o + 3 * p] * v[b3] + hg[o + 3 * p + 1] * v[b3 + 1] + hg[o + 3 * p + 2] * v[b3 + 2];
468+ }
469+ f -= hch[h] * s * hg[o + 3 * beads.indexOf(b) + c];
470+ }
351471 return f;
352472 };
353473 // At a sample boundary the damping for the time pendingD has not been
@@ -358,6 +478,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
358478 const listen = (pendingD: number) => {
359479 let s = 0;
360480 if (sys.listenQuantity !== "displacement") {
481+ refreshHinges(); // hinge damping gradients at the current positions
361482 if (pendingD !== 0 && hasDamping) {
362483 applyA(v, Av);
363484 for (let q = 0; q < N3; q++) vb[q] = v[q] - pendingD * Av[q];
@@ -395,20 +516,20 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
395516 let bv: Float64Array | null = null;
396517 let om2: Float64Array | null = null;
397518 if (modes) {
398- bz = new Float64Array(modes.free.length);
399- bv = new Float64Array(modes.free.length);
519+ bz = new Float64Array(modes.dof.length);
520+ bv = new Float64Array(modes.dof.length);
400521 om2 = new Float64Array(nModes);
401522 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
402523 }
403524 const recordModes = (f: number) => {
404525 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
405- const free = modes.free;
406- const nf = free.length;
526+ const dof = modes.dof;
527+ const nf = dof.length;
407528 for (let r = 0; r < nf; r++) {
408- const b = free[r];
409- const sm = Math.sqrt(mass[b]);
410- bz[r] = sm * (x[3 * b + 2] - pos0[3 * b + 2]);
411- bv[r] = sm * v[3 * b + 2];
529+ const q = dof[r];
530+ const sm = Math.sqrt(mass[(q / 3) | 0]);
531+ bz[r] = sm * (x[q] - pos0[q]);
532+ bv[r] = sm * v[q];
412533 }
413534 const psi = modes.psi;
414535 for (let k = 0; k < nModes; k++) {
@@ -585,6 +706,9 @@ export function estimateOmegaMax(sys: BeadSystem): number {
585706 }
586707 const ism = new Float64Array(n);
587708 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
709+ // hinges at rest: stiffness kh grad(theta) grad(theta)^T (theta = theta0 there)
710+ const hg0 = new Float64Array(12 * sys.hi.length);
711+ for (let h = 0; h < sys.hi.length; h++) hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], hg0.subarray(12 * h, 12 * h + 12));
588712 let w = new Float64Array(3 * n);
589713 let seed = 12345;
590714 for (let q = 0; q < 3 * n; q++) {
@@ -618,6 +742,14 @@ export function estimateOmegaMax(sys: BeadSystem): number {
618742 y[j3 + 1] -= f1;
619743 y[j3 + 2] -= f2;
620744 }
745+ for (let h = 0; h < sys.hi.length; h++) {
746+ const o = 12 * h;
747+ const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
748+ let s = 0;
749+ for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) s += hg0[o + 3 * p + c] * u[3 * beads[p] + c];
750+ s *= sys.hkh[h];
751+ for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) y[3 * beads[p] + c] += s * hg0[o + 3 * p + c];
752+ }
621753 for (let b = 0; b < sys.bi.length; b++) {
622754 const i3 = 3 * sys.bi[b];
623755 const j3 = 3 * sys.bj[b];
src/sim/types.tsmodified+17−3View file
@@ -32,10 +32,23 @@ export interface BeadSystem {
3232 bk: Int32Array;
3333 bkb: Float64Array;
3434 bcb: Float64Array;
35+ // Hinges across the shared edge j-k of triangles (i, j, k) and (l, k, j):
36+ // energy kh/2 (theta - theta0)^2, theta the angle between the triangles,
37+ // theta0 its value at the rest positions; damping force -ch (dtheta/dt) grad theta
38+ hi: Int32Array;
39+ hj: Int32Array;
40+ hk: Int32Array;
41+ hl: Int32Array;
42+ hkh: Float64Array;
43+ hch: Float64Array;
44+ hth0: Float64Array;
3545 // Dashpots on relative velocity: F_i = c (v_j - v_i)
3646 di: Int32Array;
3747 dj: Int32Array;
3848 dc: Float64Array;
49+ // axial part: an extra force ca (d . (v_j - v_i)) d along the current
50+ // direction d from bead i to bead j
51+ dca: Float64Array;
3952 friction: Float64Array; // n, F_i = -c0_i v_i (N s/m)
4053 hits: Hit[];
4154 listenBeads: Int32Array;
@@ -57,10 +70,11 @@ export interface SimProgress {
5770
5871 export interface Modes {
5972 freqs: Float64Array; // Hz, ascending
60- // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over free beads,
61- // column-major (nFree x nModes)
73+ // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over the degrees of
74+ // freedom in dof, column-major (dof.length x nModes)
6275 psi: Float64Array;
63- free: Int32Array; // bead index of each row of psi
76+ dof: Int32Array; // index into the 3n coordinates (3*bead + component) of each row of psi
77+ kind: "transverse" | "full"; // z motion of a flat net only, or all motion
6478 }
6579
6680 export interface SimResult {
src/ui/Plots.tsxmodified+3−3View file
@@ -186,12 +186,12 @@ export function ModeEnergy({ p }: { p: RunPayload }) {
186186 },
187187 [bands, tmax]
188188 );
189- if (!bands) return <div className="note">Mode energies are shown when linear modes are computed (flat systems only).</div>;
189+ if (!bands) return <div className="note">Mode energies are shown when linear modes are computed (see the notes above if they were skipped).</div>;
190190 return (
191191 <>
192192 <canvas ref={ref} className="plot" />
193193 <div className="legend">
194- Share of the linear (quadratic) energy in each band of transverse modes, 8 decades of color. Show the first{" "}
194+ Share of the linear (quadratic) energy in each band of {p.modeKind === "full" ? "3D" : "transverse"} modes, 8 decades of color. Show the first{" "}
195195 <select value={tmax} onChange={e => setTmax(+e.target.value)}>
196196 {[0.03, 0.1, 0.3, 1]
197197 .filter(t => t < p.sys.duration)
@@ -351,7 +351,7 @@ export function DecayTimes({ p }: { p: RunPayload }) {
351351 <>
352352 <canvas ref={ref} className="plot short" />
353353 <div className="legend">
354- Time for each transverse mode to fall 60 dB at small amplitude, from sys.friction and sys.dashpots. Use bp_decay to set it from two
354+ Time for each mode to fall 60 dB at small amplitude, from the friction, dashpots and bending and hinge damping. Use bp_decay to set it from two
355355 (frequency, T60) pairs.
356356 </div>
357357 </>
src/worker.tsmodified+17−8View file
@@ -10,7 +10,8 @@ import type { BeadSystem, SimResult } from "./sim/types";
1010 const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
1111 const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
1212
13-const MAX_MODE_BEADS = 1500;
13+const MAX_MODE_BEADS = 1500; // flat nets: z motion only
14+const MAX_FULL_DOF = 1200; // curved nets: all three components of every free bead
1415 const MAX_TWIN_BYTES = 400e6;
1516
1617 export interface RunPayload {
@@ -22,6 +23,7 @@ export interface RunPayload {
2223 sj: Int32Array;
2324 nDashpots: number;
2425 nBends: number;
26+ nHinges: number;
2527 listenQuantity: string;
2628 listenBeads: Int32Array;
2729 hitBeads: Int32Array;
@@ -30,6 +32,7 @@ export interface RunPayload {
3032 };
3133 result: SimResult;
3234 modeFreqs: Float64Array | null;
35+ modeKind: "transverse" | "full" | null;
3336 modeDecay: Float64Array | null; // small-amplitude amplitude decay rate of each mode (1/s)
3437 notes: string[];
3538 check: { dt: number; windows: { t: number; rel: number }[] } | null;
@@ -70,16 +73,20 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
7073 let modes = null;
7174 const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
7275 if (sys.modes) {
73- if (!isPlanar(sys)) notes.push("Linear modes skipped: the beads do not rest in a plane z = const.");
74- else if (nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS}).`);
75- else {
76- post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
77- modes = computeModes(sys, numblEig);
76+ const planar = isPlanar(sys);
77+ if (planar && nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS} for a flat net).`);
78+ else if (!planar && 3 * nFree > MAX_FULL_DOF) {
79+ notes.push(
80+ `Linear modes skipped: the net is curved, so its modes need all ${3 * nFree} coordinates of the free beads (limit ${MAX_FULL_DOF}, about ${MAX_FULL_DOF / 3} free beads).`
81+ );
82+ } else {
83+ post({ type: "stage", text: planar ? `Computing linear modes (${nFree} free beads)` : `Computing 3D linear modes (${3 * nFree} coordinates)` });
84+ modes = computeModes(sys, numblEig, planar ? "transverse" : "full");
7885 const fTop = modes.freqs[modes.freqs.length - 1];
7986 const nZero = modes.freqs.filter(f => f < 1e-6 * fTop).length;
8087 if (nZero > 0) {
8188 notes.push(
82- `${nZero} transverse mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` +
89+ `${nZero} mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` +
8390 "Pin more beads (a cluster rather than one bead), or tie the object to pinned anchor beads with soft springs."
8491 );
8592 }
@@ -154,8 +161,9 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
154161 pinned: sys.pinned,
155162 si: sys.si,
156163 sj: sys.sj,
157- nDashpots: sys.dc.filter(c => c !== 0).length,
164+ nDashpots: sys.dc.filter((c, d) => c !== 0 || sys.dca[d] !== 0).length,
158165 nBends: sys.bi.length,
166+ nHinges: sys.hi.length,
159167 listenQuantity: sys.listenQuantity,
160168 listenBeads: sys.listenBeads,
161169 hitBeads,
@@ -164,6 +172,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
164172 },
165173 result,
166174 modeFreqs: modes ? modes.freqs : null,
175+ modeKind: modes ? modes.kind : null,
167176 modeDecay: modes ? modalDamping(sys, modes) : null,
168177 notes,
169178 check: checkOut,
tests/bench.tsmodified+1−1View file
@@ -9,7 +9,7 @@ const { sys, warnings } = runScript(src, helpers, s => process.stdout.write(s));
99 console.log("script ms", (performance.now() - t).toFixed(0), warnings);
1010 sys.duration = Number(process.argv[3] ?? 0.2);
1111 t = performance.now();
12-const modes = isPlanar(sys) ? computeModes(sys, numblEig) : null;
12+const modes = computeModes(sys, numblEig, isPlanar(sys) ? "transverse" : "full");
1313 console.log("modes ms", (performance.now() - t).toFixed(0), modes?.freqs.slice(0, 6));
1414 const r = simulate(sys, { modes });
1515 console.log({ dt: r.dt, substeps: r.substeps, omegaMax: r.omegaMax, ms: r.elapsedMs.toFixed(0), drift: r.energyDrift, resid: r.residualForce, maxT: r.maxTension, warnings: r.warnings, aborted: r.aborted });
tests/cascade.tsmodified+2−2View file
@@ -3,12 +3,12 @@
33 import { readFileSync, readdirSync } from "fs";
44 import { numblEig, runScript } from "../src/sim/runScript";
55 import { simulate } from "../src/sim/solver";
6-import { computeModes } from "../src/sim/modes";
6+import { computeModes, isPlanar } from "../src/sim/modes";
77 const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
88 const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
99 sys.duration = Number(process.argv[3]);
1010 const fc = Number(process.argv[4]);
11-const modes = computeModes(sys, numblEig);
11+const modes = computeModes(sys, numblEig, isPlanar(sys) ? "transverse" : "full");
1212 const nm = modes.freqs.length;
1313 const scaleHits = (s: number) => sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * s) as [number, number, number], force: h.force.map(v => v * s) as [number, number, number] }));
1414 const share = (r: ReturnType<typeof simulate>, t: number) => {
tests/solver.test.tsmodified+108−4View file
@@ -1,4 +1,5 @@
11 import { readFileSync, readdirSync } from "fs";
2+import { executeCode } from "numbl";
23 import { describe, expect, it } from "vitest";
34 import { computeModes, modalDamping } from "../src/sim/modes";
45 import { numblEig, runScript } from "../src/sim/runScript";
@@ -206,7 +207,7 @@ for i = 1:N
206207 sys.hits(i).beads = i + 1;
207208 sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
208209 end
209-[sys.friction, sys.dashpots] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1]);
210+sys = bp_decay(sys, [100 0.5; 1000 0.1]);
210211 sys.listen.beads = 2:N+1;
211212 sys.listen.dir = [0 0 1];
212213 `).replace("sys.listen = struct('beads', 4);", "")
@@ -241,7 +242,7 @@ sys.listen.dir = [0 0 1];
241242 });
242243
243244 it("rejects a law the dashpots cannot make", () => {
244- expect(() => build(STRING("[f, d] = bp_decay(pos, m, sys.springs, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
245+ expect(() => build(STRING("sys = bp_decay(sys, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
245246 });
246247 });
247248
@@ -289,7 +290,7 @@ ${extra}`);
289290 });
290291
291292 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 { sys } = build(STIFF("sys = bp_decay(sys, [100 0.5; 1000 0.1]);"));
293294 const modes = computeModes(sys, numblEig);
294295 const md = modalDamping(sys, modes);
295296 const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
@@ -344,7 +345,7 @@ describe("listener quantity", () => {
344345 STRING(`
345346 [pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
346347 sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
347-[sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, m, sys.springs, [100 0.05; 1000 0.02], sys.bends);
348+sys = bp_decay(sys, [100 0.05; 1000 0.02]);
348349 sys.hits.velocity = [0 0 1e-4];
349350 sys.listen.quantity = '${q}';`)
350351 );
@@ -380,3 +381,106 @@ sys.listen.quantity = '${q}';`)
380381 expect(() => build(STRING("sys.listen.quantity = 'jerk';"))).toThrow(/quantity/);
381382 });
382383 });
384+
385+describe("hinges", () => {
386+ // A small dome: hexagonal lattice lifted to a paraboloid, clamped at the top,
387+ // with stretching springs (no tension) and hinges across interior edges.
388+ const DOME = (extra: string, H = 0.01) => `
389+R = 0.06; a = 0.012; H = ${H}; m = 1e-3; kh = 0.5; k = 2e4;
390+[pos, tri] = bp_hexdisk(R, a);
391+pos(:, 3) = H * (1 - sum(pos(:, 1:2).^2, 2) / R^2);
392+E = bp_edges(tri);
393+sys.pos = pos; sys.mass = m;
394+sys.pinned = bp_near(pos, [0 0 H], 1.05 * a);
395+sys.springs = bp_springs(pos, E, k, 0);
396+hh = bp_hinges(tri);
397+sys.hinges = [hh, kh * ones(size(hh, 1), 1)];
398+hit.beads = bp_near(pos, [0.04 0 0], 0); hit.velocity = [0 0 1e-3]; sys.hits = hit;
399+sys.listen = struct('beads', bp_near(pos, [-0.03 0.02 0], 0));
400+sys.friction = 0; sys.duration = 0.05;
401+${extra}`;
402+
403+ it("bp_hinges lists every interior edge once, with the two opposite beads", () => {
404+ const r = executeCode(
405+ `[pos, tri] = bp_hexdisk(0.06, 0.012); E = bp_edges(tri); H = bp_hinges(tri);
406+ nint = size(E, 1) - sum(accumarray(sort([tri(:,[1 2]); tri(:,[2 3]); tri(:,[3 1])], 2) * [1e4; 1], 1) == 1);`,
407+ {},
408+ helpers,
409+ "t.m"
410+ );
411+ const H = r.variableValues.H as unknown as { data: Float64Array; shape: number[] };
412+ const tri = r.variableValues.tri as unknown as { data: Float64Array; shape: number[] };
413+ const nh = H.shape[0];
414+ const nt = tri.shape[0];
415+ const tris = new Set<string>();
416+ for (let t = 0; t < nt; t++) tris.add([0, 1, 2].map(c => tri.data[t + c * nt]).sort((p, q) => p - q).join(","));
417+ for (let h = 0; h < nh; h++) {
418+ const [i, j, k, l] = [0, 1, 2, 3].map(c => H.data[h + c * nh]);
419+ expect(tris.has([i, j, k].sort((p, q) => p - q).join(","))).toBe(true);
420+ expect(tris.has([l, j, k].sort((p, q) => p - q).join(","))).toBe(true);
421+ }
422+ expect(nh).toBe(r.variableValues.nint);
423+ });
424+
425+ it("leaves a dome, and a rotated dome, unstressed at rest", () => {
426+ for (const rot of ["", "Q = [cos(0.7) 0 sin(0.7); 0 1 0; -sin(0.7) 0 cos(0.7)]; sys.pos = sys.pos * Q';"]) {
427+ const { sys } = build(DOME(rot));
428+ const r = simulate({ ...sys, duration: 0.001 });
429+ expect(r.residualForce).toBeLessThan(1e-9);
430+ expect(r.warnings).toEqual([]);
431+ }
432+ });
433+
434+ it("conserves energy to second order", () => {
435+ const { sys } = build(DOME("sys.hits.velocity = [1 0.5 3];"));
436+ const d1 = simulate({ ...sys, dt: 1 / 48000 / 8 }).energyDrift;
437+ const d2 = simulate({ ...sys, dt: 1 / 48000 / 16 }).energyDrift;
438+ expect(d1).toBeLessThan(1e-3);
439+ expect(d1 / d2).toBeGreaterThan(3);
440+ expect(d1 / d2).toBeLessThan(5);
441+ });
442+
443+ it("gives a flat net the same transverse modes in the fast and the full analysis", () => {
444+ const { sys } = build(DOME("", 0));
445+ const fz = computeModes(sys, numblEig, "transverse").freqs;
446+ const full = Array.from(computeModes(sys, numblEig, "full").freqs);
447+ for (const f of fz) expect(Math.min(...full.map(g => Math.abs(g / f - 1)))).toBeLessThan(1e-8);
448+ });
449+
450+ it("gives every 3D mode of a dome the bp_decay law", () => {
451+ const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
452+ const modes = computeModes(sys, numblEig, "full");
453+ const sigma = modalDamping(sys, modes);
454+ const s = [Math.log(1000) / 0.4, Math.log(1000) / 0.1];
455+ const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 3000) ** 2];
456+ const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
457+ const sigma0 = s[0] - sigma1 * w2[0];
458+ modes.freqs.forEach((f, k) => expect(sigma[k] / (sigma0 + sigma1 * (2 * Math.PI * f) ** 2)).toBeCloseTo(1, 8));
459+ });
460+
461+ it.each([6, -3])("keeps a dome's energy in one 3D mode (%i), decaying at the predicted rate", kSel => {
462+ const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
463+ const modes = computeModes(sys, numblEig, "full");
464+ const sigma = modalDamping(sys, modes);
465+ const k = kSel >= 0 ? kSel : modes.freqs.length + kSel; // -3: one of the highest, stretching modes
466+ const nd = modes.dof.length;
467+ // initial velocity in the shape of mode k: v = eps * psi_k / sqrt(m)
468+ const v = new Float64Array(3 * sys.n);
469+ for (let r = 0; r < nd; r++) v[modes.dof[r]] = (1e-4 * modes.psi[k * nd + r]) / Math.sqrt(sys.mass[(modes.dof[r] / 3) | 0]);
470+ const hits = Array.from({ length: sys.n }, (_, b) => ({
471+ t: 0,
472+ beads: Int32Array.of(b),
473+ velocity: [v[3 * b], v[3 * b + 1], v[3 * b + 2]] as [number, number, number],
474+ force: [0, 0, 0] as [number, number, number],
475+ duration: 0,
476+ })).filter(h => !sys.pinned[h.beads[0]]);
477+ const r = simulate({ ...sys, hits, duration: 0.1, frameRate: 1000, dt: 1 / 48000 / 8 }, { modes });
478+ const nm = modes.freqs.length;
479+ const at = (t: number) => r.frameTimes.findIndex(x => x >= t);
480+ const share = (f: number) => r.modeEnergy![f * nm + k] / r.modeEnergy!.slice(f * nm, f * nm + nm).reduce((p, q) => p + q, 0);
481+ expect(share(at(0.05))).toBeGreaterThan(0.999);
482+ const e = (f: number) => r.modeEnergy![f * nm + k];
483+ const rate = Math.log(e(at(0.02)) / e(at(0.08))) / (r.frameTimes[at(0.08)] - r.frameTimes[at(0.02)]) / 2;
484+ expect(Math.abs(rate / sigma[k] - 1)).toBeLessThan(0.02);
485+ });
486+});