Cymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphone
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012zsyegMfQnzJEzgMiFLzew
16 changed files+1323−154
README.mdmodified+11−0View file
@@ -34,6 +34,17 @@ String tests (`test/string.test.ts`): modes vs the discrete-string formula; smal
3434 - Presets: a fast tom with in-plane stiffness softened ~70x (about 8x slower than real time for 571 beads in Chromium) and a realistic Mylar one (~6x more steps).
3535 - Tests (`test/drum.test.ts`): every bead has six springs; a small centre press rings within 1.3% of the ideal membrane fundamental 2.405 c / (2 pi R); the mallet bounces off and energy balances.
3636
37+## Cymbal
38+
39+- Thin bronze plate (E = 100 GPa, 8600 kg/m^3, 0.8 mm by default), round, free rim, centre hole; flat by default, optional bow (dome, y = h (1 - (r/R)^2)) and bell (raised cosine). Beads on a triangular lattice lifted to the profile; the mesh is scaled so its area equals the disk's.
40+- In-plane: springs with k = sqrt(3)/2 E t / softening (an edge with one triangle gets half). Softening (default 30x in the fast preset) speeds the run; bending is unaffected.
41+- Bending: dihedral angles across interior edges (`src/sim/hinge.ts`, Bridson et al. 2003), sign-normalized so a positive change folds toward the upward normal. Each triangle builds a shape operator S = sum_e phi_e l_e / (2A) m_e m_e^T and stores the plate energy A D/2 [(tr S)^2 - 2 (1 - nu) det S] (nu = 1/3) as a 3x3 matrix in its hinge angles, condensed over free edges so they carry no moment. Plain hinge springs get the twisting term wrong on a regular mesh (effective nu = -1/3), which matters for a free rim. Damping uses the same matrices (Rayleigh, sigma1 omega^2).
42+- Centre support: beads within the felt radius get six short springs with dashpots to fixed points above, below and to the sides.
43+- Excitation: drumstick (Hertz, K delta^1.5) or felt mallet as a 3D point mass with contact along the surface normal, regularized Coulomb friction and a contact reach; hit angle tilts the stroke toward the rim. Also strike and force pulse.
44+- Sound: far-field microphone via the Rayleigh integral, p = rho / (2 pi) sum A_i a_n,i(t - R_i / c) / R_i (baffled-plate assumption: no diffraction around the rim).
45+- Tests (`test/cymbal.test.ts`, `test/hinge.test.ts`): hinge gradient vs finite differences; bending energy matches the plate for cylinder, sphere, saddle and twist to 0.1%; a free flat plate's (2,0) mode within 3.4% of Leissa's value on a 9-ring lattice; energy balance with a glancing stick, friction, felts and damping.
46+- Limits: about 45x slower than real time for ~370 beads in Chromium; bending is resolved only to ~1.3 kHz at the default resolution, so it sounds like a small gong, not a crash. A single-bead stick contact chatters on this lattice, so presets spread the contact over ~1.5 bead spacings.
47+
3748 ## Layout
3849
3950 - `src/sim/` generic 3D bead-spring system, solver, worker.
src/instruments/bead-cymbal/model.tsadded+460−0View file
@@ -0,0 +1,460 @@
1+import type { Anchor, BeadSystem, Bending, Excitation } from '../../sim/types';
2+import { hingeAngleGrad } from '../../sim/hinge';
3+import { dampingFromT60 } from '../bead-string/model';
4+
5+// A cymbal: a thin bronze plate, round, free at the rim, with a hole in the
6+// centre where it sits between felt washers. Beads on a triangular lattice in
7+// the x-z plane, lifted to the profile height y(r) (flat by default; optional
8+// bow and bell). In-plane springs on every edge; plate bending from the
9+// dihedral angles across interior edges, per triangle (see Bending in sim/types). Vibration is mostly along y.
10+
11+export type CymbalExcitation = 'stick' | 'mallet' | 'strike' | 'force';
12+
13+export interface CymbalParams {
14+ // geometry
15+ diameter: number; // m
16+ rings: number; // beads from the centre to the rim
17+ holeRadius: number; // m
18+ bowHeight: number; // m, dome rise from rim to centre (0 = flat)
19+ bellDiameter: number; // m
20+ bellHeight: number; // m (0 = no bell)
21+ // material
22+ thickness: number; // m
23+ youngs: number; // Pa
24+ density: number; // kg/m^3
25+ soften: number; // in-plane stiffness divided by this (1 = real bronze)
26+ // centre support
27+ feltRadius: number; // m, beads within this radius sit in the felts
28+ feltK: number; // N/m, total vertical stiffness
29+ feltC: number; // N s/m, total damping
30+ feltGap: number; // m, felt thickness (spring length)
31+ // damping
32+ t60Low: number;
33+ fLow: number;
34+ t60High: number;
35+ fHigh: number;
36+ // excitation
37+ excitation: CymbalExcitation;
38+ excR: number; // fraction of radius
39+ excAngle: number; // degrees, azimuth in the x-z plane
40+ hitAngle: number; // degrees from the surface normal, tilted toward the rim
41+ excWidth: number; // m, Gaussian sigma of the patch (0 = one bead)
42+ stickMass: number;
43+ stickSpeed: number;
44+ stickK: number; // N/m^1.5 (Hertz)
45+ malletMass: number;
46+ malletSpeed: number;
47+ malletK: number; // N/m^p
48+ malletP: number;
49+ friction: number;
50+ strikeVel: number;
51+ forcePeak: number;
52+ forceMs: number;
53+ // microphone
54+ micDistance: number; // m
55+ micElevation: number; // degrees above the plane
56+ micAzimuth: number; // degrees
57+ // simulation
58+ duration: number;
59+ fs: number;
60+ cfl: number;
61+ frameRate: number;
62+}
63+
64+export function profileHeight(p: CymbalParams, r: number): number {
65+ const R = p.diameter / 2;
66+ let y = p.bowHeight * (1 - (r / R) ** 2);
67+ const rb = p.bellDiameter / 2;
68+ if (p.bellHeight > 0 && rb > 0 && r < rb) y += p.bellHeight * 0.5 * (1 + Math.cos((Math.PI * r) / rb));
69+ return y;
70+}
71+
72+export interface CymbalMesh {
73+ a: number;
74+ pos: [number, number][]; // (x, z)
75+ tris: [number, number, number][];
76+ edges: { i: number; j: number; opp: number[] }[]; // opposite vertices of the 1 or 2 triangles
77+}
78+
79+export function cymbalMesh(R: number, rings: number, hole: number): CymbalMesh {
80+ const a = R / (rings + 0.5);
81+ const at = (i: number, j: number): [number, number] => [a * (i + 0.5 * j), a * (Math.sqrt(3) / 2) * j];
82+ const index = new Map<string, number>();
83+ const pos: [number, number][] = [];
84+ const M = rings + 2;
85+ for (let j = -M; j <= M; j++)
86+ for (let i = -M - Math.abs(j); i <= M + Math.abs(j); i++) {
87+ const r = Math.hypot(...at(i, j));
88+ if (r <= R - 0.5 * a + 1e-9 && !(r < hole)) {
89+ index.set(`${i},${j}`, pos.length);
90+ pos.push(at(i, j));
91+ }
92+ }
93+ const id = (i: number, j: number) => index.get(`${i},${j}`);
94+ const tris: [number, number, number][] = [];
95+ for (const k of index.keys()) {
96+ const [i, j] = k.split(',').map(Number);
97+ const t1 = [id(i, j), id(i + 1, j), id(i, j + 1)];
98+ const t2 = [id(i + 1, j), id(i + 1, j + 1), id(i, j + 1)];
99+ for (const t of [t1, t2]) if (t.every((v) => v !== undefined)) tris.push(t as [number, number, number]);
100+ }
101+ const em = new Map<string, { i: number; j: number; opp: number[] }>();
102+ for (const [x, y, z] of tris)
103+ for (const [u, v, o] of [[x, y, z], [y, z, x], [z, x, y]]) {
104+ const key = u < v ? `${u},${v}` : `${v},${u}`;
105+ let e = em.get(key);
106+ if (!e) em.set(key, (e = { i: Math.min(u, v), j: Math.max(u, v), opp: [] }));
107+ e.opp.push(o);
108+ }
109+ // scale so the triangulated area equals the disk's (minus the hole): the
110+ // outermost row of beads otherwise leaves the plate about a/2 short
111+ let meshArea = 0;
112+ for (const [u, v, w] of tris) {
113+ const [p, q, r] = [pos[u], pos[v], pos[w]];
114+ meshArea += 0.5 * Math.abs((q[0] - p[0]) * (r[1] - p[1]) - (r[0] - p[0]) * (q[1] - p[1]));
115+ }
116+ const f = Math.sqrt((Math.PI * (R * R - hole * hole)) / meshArea);
117+ for (const p of pos) {
118+ p[0] *= f;
119+ p[1] *= f;
120+ }
121+ return { a: a * f, pos, tris, edges: [...em.values()] };
122+}
123+
124+export interface CymbalDerived {
125+ mesh: CymbalMesh;
126+ n: number;
127+ rest: Float64Array; // 3n
128+ area: Float64Array; // n, lumped area per bead
129+ normal: Float64Array; // 3n, upward unit normals
130+ mass: Float64Array;
131+ k: number; // interior spring stiffness
132+ D: number; // bending stiffness
133+ support: number[]; // beads in the felts
134+ sigma0: number;
135+ sigma1: number;
136+ waveSpeedP: number; // in-plane P waves (after softening)
137+ bendingCoef: number; // sqrt(D / rho t): bending waves omega = coef * k^2
138+ warnings: string[];
139+}
140+
141+// Triangular lattice, spacing a, Poisson ratio 1/3:
142+// in-plane: spring stiffness k = sqrt(3)/2 E t (an edge with one triangle gets half)
143+// bending: D = E t^3 / (12 (1 - nu^2)), through per-triangle shape operators
144+export function deriveCymbal(p: CymbalParams): CymbalDerived {
145+ const R = p.diameter / 2;
146+ const mesh = cymbalMesh(R, Math.max(2, Math.round(p.rings)), p.holeRadius);
147+ const n = mesh.pos.length;
148+ const rest = new Float64Array(3 * n);
149+ mesh.pos.forEach(([x, z], i) => {
150+ rest[3 * i] = x;
151+ rest[3 * i + 1] = profileHeight(p, Math.hypot(x, z));
152+ rest[3 * i + 2] = z;
153+ });
154+ const area = new Float64Array(n);
155+ const normal = new Float64Array(3 * n);
156+ for (const [u, v, w] of mesh.tris) {
157+ const e1 = [0, 1, 2].map((c) => rest[3 * v + c] - rest[3 * u + c]);
158+ const e2 = [0, 1, 2].map((c) => rest[3 * w + c] - rest[3 * u + c]);
159+ let nx = e1[1] * e2[2] - e1[2] * e2[1];
160+ let ny = e1[2] * e2[0] - e1[0] * e2[2];
161+ let nz = e1[0] * e2[1] - e1[1] * e2[0];
162+ const A2 = Math.hypot(nx, ny, nz);
163+ if (ny < 0) [nx, ny, nz] = [-nx, -ny, -nz];
164+ for (const b of [u, v, w]) {
165+ area[b] += A2 / 6;
166+ normal[3 * b] += nx;
167+ normal[3 * b + 1] += ny;
168+ normal[3 * b + 2] += nz;
169+ }
170+ }
171+ for (let i = 0; i < n; i++) {
172+ const l = Math.hypot(normal[3 * i], normal[3 * i + 1], normal[3 * i + 2]) || 1;
173+ for (let c = 0; c < 3; c++) normal[3 * i + c] /= l;
174+ }
175+ const rhoT = p.density * p.thickness;
176+ const mass = area.map((A) => rhoT * A);
177+ const Et = (p.youngs * p.thickness) / Math.max(p.soften, 1e-9);
178+ const k = (Math.sqrt(3) / 2) * Et;
179+ const nu = 1 / 3;
180+ const D = (p.youngs * p.thickness ** 3) / (12 * (1 - nu * nu));
181+ const support: number[] = [];
182+ mesh.pos.forEach(([x, z], i) => {
183+ if (Math.hypot(x, z) <= p.feltRadius + 1e-9) support.push(i);
184+ });
185+ const warnings: string[] = [];
186+ if (!support.length) warnings.push('No beads sit inside the felt radius: the cymbal is unsupported and will drift. Increase the felt radius.');
187+ if (p.holeRadius >= p.feltRadius) warnings.push('The centre hole is as large as the felt radius.');
188+ if (p.soften > 1)
189+ warnings.push(
190+ `In-plane stiffness is softened ${p.soften}×. A flat plate bends the same at small amplitude; curved shapes lose curvature stiffening and big hits exaggerate the nonlinear coupling.`,
191+ );
192+ const { sigma0, sigma1, note } = dampingFromT60(p);
193+ if (note) warnings.push(note);
194+ const lambdaPlus2Mu = (3 * Math.sqrt(3) * k) / 4;
195+ return {
196+ mesh,
197+ n,
198+ rest,
199+ area,
200+ normal,
201+ mass,
202+ k,
203+ D,
204+ support,
205+ sigma0,
206+ sigma1,
207+ waveSpeedP: Math.sqrt(lambdaPlus2Mu / rhoT),
208+ bendingCoef: Math.sqrt(D / rhoT),
209+ warnings,
210+ };
211+}
212+
213+// highest frequency whose bending wavelength spans six lattice spacings (where
214+// the lattice is still a fair plate)
215+export function bendingLimitHz(d: CymbalDerived): number {
216+ const kb = (2 * Math.PI) / (6 * d.mesh.a);
217+ return (d.bendingCoef * kb * kb) / (2 * Math.PI);
218+}
219+
220+export function buildCymbal(p: CymbalParams, d: CymbalDerived): BeadSystem {
221+ const { n, rest, mesh } = d;
222+ const si: number[] = [];
223+ const sj: number[] = [];
224+ const sk: number[] = [];
225+ const sL0: number[] = [];
226+ const sc: number[] = [];
227+ const len = (i: number, j: number) =>
228+ Math.hypot(rest[3 * j] - rest[3 * i], rest[3 * j + 1] - rest[3 * i + 1], rest[3 * j + 2] - rest[3 * i + 2]);
229+ for (const e of mesh.edges) {
230+ si.push(e.i);
231+ sj.push(e.j);
232+ sk.push((d.k / 2) * e.opp.length);
233+ sL0.push(len(e.i, e.j));
234+ sc.push(0);
235+ }
236+ // felt washers: each supported bead gets six short springs to fixed points
237+ // above, below and to the sides; dashpots on all of them
238+ const anchors: Anchor[] = [];
239+ const ns = d.support.length;
240+ const kf = ns ? p.feltK / (2 * ns) : 0;
241+ const cf = ns ? p.feltC / (6 * ns) : 0;
242+ for (const i of d.support) {
243+ for (let c = 0; c < 3; c++)
244+ for (const s of [1, -1]) {
245+ const pos: [number, number, number] = [rest[3 * i], rest[3 * i + 1], rest[3 * i + 2]];
246+ pos[c] += s * p.feltGap;
247+ anchors.push({ pos, slideBead: -1 });
248+ si.push(i);
249+ sj.push(n + anchors.length - 1);
250+ sk.push(kf);
251+ sL0.push(p.feltGap);
252+ sc.push(cf);
253+ }
254+ }
255+ const bending = buildBending(d, 1 / 3);
256+ return {
257+ n,
258+ rest,
259+ mass: d.mass,
260+ pinned: new Uint8Array(n),
261+ anchors,
262+ si: Int32Array.from(si),
263+ sj: Int32Array.from(sj),
264+ sk: Float64Array.from(sk),
265+ sL0: Float64Array.from(sL0),
266+ sc: Float64Array.from(sc),
267+ sigma0: d.sigma0,
268+ bending,
269+ };
270+}
271+
272+export function hitPoint(p: CymbalParams): [number, number] {
273+ const r = (p.diameter / 2) * p.excR;
274+ const th = (p.excAngle * Math.PI) / 180;
275+ return [r * Math.cos(th), r * Math.sin(th)];
276+}
277+
278+function nearest(d: CymbalDerived, x: number, z: number): number {
279+ let best = 0;
280+ let bd = Infinity;
281+ d.mesh.pos.forEach(([px, pz], i) => {
282+ const q = (px - x) ** 2 + (pz - z) ** 2;
283+ if (q < bd) [bd, best] = [q, i];
284+ });
285+ return best;
286+}
287+
288+// patch weights (sum 1), the inward normal there and the hit direction
289+// (tilted from the normal toward the rim by hitAngle)
290+export function hitGeometry(p: CymbalParams, d: CymbalDerived) {
291+ const [cx, cz] = hitPoint(p);
292+ const w = new Float64Array(d.n);
293+ if (p.excWidth > 0)
294+ d.mesh.pos.forEach(([x, z], i) => {
295+ const r2 = ((x - cx) ** 2 + (z - cz) ** 2) / (p.excWidth * p.excWidth);
296+ w[i] = r2 < 9 ? Math.exp(-0.5 * r2) : 0;
297+ });
298+ if (!w.some((v) => v > 0)) w[nearest(d, cx, cz)] = 1;
299+ const sum = w.reduce((s, v) => s + v, 0);
300+ for (let i = 0; i < d.n; i++) w[i] /= sum;
301+ const up = [0, 0, 0];
302+ for (let i = 0; i < d.n; i++) for (let c = 0; c < 3; c++) up[c] += w[i] * d.normal[3 * i + c];
303+ const ul = Math.hypot(...up);
304+ const N = up.map((v) => v / ul);
305+ // outward tangent: radial direction minus its normal part
306+ const rl = Math.hypot(cx, cz) || 1;
307+ const er = [cx / rl, 0, cz / rl];
308+ const dot = er[0] * N[0] + er[2] * N[2];
309+ let t = [er[0] - dot * N[0], -dot * N[1], er[2] - dot * N[2]];
310+ const tl = Math.hypot(...t) || 1;
311+ t = t.map((v) => v / tl);
312+ const al = (p.hitAngle * Math.PI) / 180;
313+ const inward: [number, number, number] = [-N[0], -N[1], -N[2]];
314+ const dir: [number, number, number] = [0, 1, 2].map((c) => Math.cos(al) * inward[c] + Math.sin(al) * t[c]) as [number, number, number];
315+ return { w, inward, dir };
316+}
317+
318+export function buildCymbalExcitation(p: CymbalParams, d: CymbalDerived): Excitation {
319+ const n = d.n;
320+ const u0 = new Float64Array(3 * n);
321+ const v0 = new Float64Array(3 * n);
322+ const { w, inward, dir } = hitGeometry(p, d);
323+ const wmax = Math.max(...w);
324+ if (p.excitation === 'strike') {
325+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) v0[3 * i + c] = (p.strikeVel * w[i] * dir[c]) / wmax;
326+ return { u0, v0 };
327+ }
328+ if (p.excitation === 'force') {
329+ const shape = new Float64Array(3 * n);
330+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) shape[3 * i + c] = w[i] * dir[c];
331+ return { u0, v0, force: { shape, F0: p.forcePeak, tau: p.forceMs / 1000 } };
332+ }
333+ const stick = p.excitation === 'stick';
334+ const speed = stick ? p.stickSpeed : p.malletSpeed;
335+ return {
336+ u0,
337+ v0,
338+ mallet: {
339+ weights: w,
340+ normal: inward,
341+ vel0: [dir[0] * speed, dir[1] * speed, dir[2] * speed],
342+ mass: stick ? p.stickMass : p.malletMass,
343+ K: stick ? p.stickK : p.malletK,
344+ p: stick ? 1.5 : p.malletP,
345+ mu: p.friction,
346+ reach: Math.max(2.5 * p.excWidth, 0.75 * d.mesh.a),
347+ },
348+ };
349+}
350+
351+const RHO_AIR = 1.2;
352+const C_AIR = 343;
353+
354+export function micDirection(p: CymbalParams): [number, number, number] {
355+ const el = (p.micElevation * Math.PI) / 180;
356+ const az = (p.micAzimuth * Math.PI) / 180;
357+ return [Math.cos(el) * Math.cos(az), Math.sin(el), Math.cos(el) * Math.sin(az)];
358+}
359+
360+// Rayleigh integral: p(t) = rho / (2 pi) sum_i A_i a_n,i(t - R_i / c) / R_i,
361+// with a_n the acceleration along the normal on the microphone's side.
362+export function buildMic(p: CymbalParams, d: CymbalDerived): { weights: Float64Array; delay: Float64Array } {
363+ const [mx, my, mz] = micDirection(p).map((v) => v * p.micDistance);
364+ const side = p.micElevation >= 0 ? 1 : -1;
365+ const weights = new Float64Array(3 * d.n);
366+ const delay = new Float64Array(d.n);
367+ for (let i = 0; i < d.n; i++) {
368+ const R = Math.hypot(mx - d.rest[3 * i], my - d.rest[3 * i + 1], mz - d.rest[3 * i + 2]);
369+ delay[i] = (R / C_AIR) * p.fs;
370+ const g = (side * RHO_AIR * d.area[i]) / (2 * Math.PI * R);
371+ for (let c = 0; c < 3; c++) weights[3 * i + c] = g * d.normal[3 * i + c];
372+ }
373+ return { weights, delay };
374+}
375+
376+// Hinges on interior edges, oriented so a positive change folds the surface
377+// toward its upward normal; per-triangle plate-energy matrices, condensed over
378+// free edges (the missing angle takes the value that leaves no moment there).
379+function buildBending(d: CymbalDerived, nu: number): Bending {
380+ const { rest, mesh, normal } = d;
381+ const key = (i: number, j: number) => (i < j ? `${i},${j}` : `${j},${i}`);
382+ const hingeOf = new Map<string, number>();
383+ const w0: number[] = [], w1: number[] = [], e0: number[] = [], e1: number[] = [];
384+ for (const e of mesh.edges)
385+ if (e.opp.length === 2) {
386+ hingeOf.set(key(e.i, e.j), w0.length);
387+ w0.push(e.opp[0]);
388+ w1.push(e.opp[1]);
389+ e0.push(e.i);
390+ e1.push(e.j);
391+ }
392+ const nh = w0.length;
393+ const g = new Float64Array(13);
394+ const theta0 = new Float64Array(nh);
395+ const sign = new Float64Array(nh);
396+ const lift = 1e-4 * mesh.a;
397+ for (let e = 0; e < nh; e++) {
398+ theta0[e] = hingeAngleGrad(rest, w0[e], w1[e], e0[e], e1[e], g);
399+ const Y = rest.slice();
400+ for (const w of [w0[e], w1[e]]) for (let c = 0; c < 3; c++) Y[3 * w + c] += lift * normal[3 * w + c];
401+ sign[e] = Math.sign(hingeAngleGrad(Y, w0[e], w1[e], e0[e], e1[e], g) - theta0[e]) || 1;
402+ }
403+ const nt = mesh.tris.length;
404+ const tri = new Int32Array(3 * nt).fill(-1);
405+ const M = new Float64Array(9 * nt);
406+ const P = (i: number) => [rest[3 * i], rest[3 * i + 1], rest[3 * i + 2]];
407+ mesh.tris.forEach(([u, v, w], t) => {
408+ const ed = [[v, w], [w, u], [u, v]];
409+ const vec = ed.map(([i, j]) => P(j).map((x, c) => x - P(i)[c]));
410+ const len = vec.map((q) => Math.hypot(...q));
411+ const tu = vec.map((q, k) => q.map((x) => x / len[k]));
412+ const [x1, x2] = [vec[2], vec[1].map((x) => -x)]; // u->v and u->w
413+ const A = 0.5 * Math.hypot(x1[1] * x2[2] - x1[2] * x2[1], x1[2] * x2[0] - x1[0] * x2[2], x1[0] * x2[1] - x1[1] * x2[0]);
414+ const beta = len.map((l) => l / (2 * A));
415+ const full = new Float64Array(9);
416+ for (let a = 0; a < 3; a++)
417+ for (let b = 0; b < 3; b++) {
418+ const c2 = (tu[a][0] * tu[b][0] + tu[a][1] * tu[b][1] + tu[a][2] * tu[b][2]) ** 2;
419+ full[3 * a + b] = A * d.D * beta[a] * beta[b] * (nu + (1 - nu) * c2);
420+ }
421+ const have = ed.map(([i, j]) => hingeOf.get(key(i, j)) ?? -1);
422+ have.forEach((h, a) => (tri[3 * t + a] = h));
423+ const K = [0, 1, 2].filter((a) => have[a] >= 0);
424+ const B = [0, 1, 2].filter((a) => have[a] < 0);
425+ // condense: M' = M_KK - M_KB M_BB^-1 M_BK
426+ const inv = invertSmall(B.map((a) => B.map((b) => full[3 * a + b])));
427+ for (const a of K)
428+ for (const b of K) {
429+ let m = full[3 * a + b];
430+ for (let x = 0; x < B.length; x++)
431+ for (let y = 0; y < B.length; y++) m -= full[3 * a + B[x]] * inv[x][y] * full[3 * B[y] + b];
432+ M[9 * t + 3 * a + b] = m;
433+ }
434+ });
435+ return {
436+ w0: Int32Array.from(w0),
437+ w1: Int32Array.from(w1),
438+ e0: Int32Array.from(e0),
439+ e1: Int32Array.from(e1),
440+ sign,
441+ theta0,
442+ tri,
443+ M,
444+ sigma1: d.sigma1,
445+ };
446+}
447+
448+function invertSmall(m: number[][]): number[][] {
449+ if (m.length === 0) return [];
450+ if (m.length === 1) return [[1 / m[0][0]]];
451+ if (m.length === 2) {
452+ const det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
453+ return [
454+ [m[1][1] / det, -m[0][1] / det],
455+ [-m[1][0] / det, m[0][0] / det],
456+ ];
457+ }
458+ // all three edges free: an isolated triangle has no bending energy
459+ return m.map((r) => r.map(() => 0));
460+}
src/instruments/bead-cymbal/page.tsadded+189−0View file
@@ -0,0 +1,189 @@
1+import type { Section } from '../../ui/controls';
2+import { fmtNum } from '../../ui/plot';
3+import { MembraneView } from '../../ui/membraneview';
4+import { mountInstrumentPage, type Built } from '../../ui/instrumentPage';
5+import {
6+ bendingLimitHz,
7+ buildCymbal,
8+ buildCymbalExcitation,
9+ buildMic,
10+ deriveCymbal,
11+ hitPoint,
12+ micDirection,
13+ profileHeight,
14+ type CymbalParams,
15+} from './model';
16+import { CYMBAL_PRESETS } from './presets';
17+
18+const isExc = (...t: CymbalParams['excitation'][]) => (p: CymbalParams) => t.includes(p.excitation);
19+
20+const SECTIONS: Section<CymbalParams>[] = [
21+ {
22+ title: 'Geometry',
23+ params: [
24+ { key: 'diameter', label: 'Diameter', kind: 'number', unit: 'cm', scale: 100 },
25+ { key: 'rings', label: 'Rings of beads', kind: 'int', min: 2, help: 'Beads from the centre to the rim; bead count grows as about 3 × rings². More rings resolve higher frequencies but cost much more.' },
26+ { key: 'holeRadius', label: 'Centre hole radius', kind: 'number', unit: 'mm', scale: 1000 },
27+ { key: 'bowHeight', label: 'Bow height', kind: 'number', unit: 'mm', scale: 1000, help: 'Rise of the dome from the rim to the centre. 0 = flat.' },
28+ { key: 'bellDiameter', label: 'Bell diameter', kind: 'number', unit: 'cm', scale: 100 },
29+ { key: 'bellHeight', label: 'Bell height', kind: 'number', unit: 'mm', scale: 1000, help: '0 = no bell.' },
30+ ],
31+ },
32+ {
33+ title: 'Material',
34+ params: [
35+ { key: 'thickness', label: 'Thickness', kind: 'number', unit: 'mm', scale: 1000 },
36+ { key: 'youngs', label: "Young's modulus", kind: 'number', unit: 'GPa', scale: 1e-9, help: 'B20 bronze: about 100 GPa.' },
37+ { key: 'density', label: 'Density', kind: 'number', unit: 'kg/m³', help: 'B20 bronze: about 8600 kg/m³.' },
38+ { key: 'soften', label: 'In-plane softening', kind: 'number', unit: '×', help: 'Divide the in-plane stiffness by this for speed (bending stays real). 1 = real bronze. Flat plates are unaffected at small amplitude; curved shapes and big hits are.' },
39+ ],
40+ },
41+ {
42+ title: 'Centre support (felts)',
43+ params: [
44+ { key: 'feltRadius', label: 'Felt radius', kind: 'number', unit: 'mm', scale: 1000, help: 'Beads within this radius sit between the felt washers.' },
45+ { key: 'feltK', label: 'Felt stiffness', kind: 'number', unit: 'N/m', help: 'Total vertical stiffness of the felts (sideways the same).' },
46+ { key: 'feltC', label: 'Felt damping', kind: 'number', unit: 'N·s/m', help: 'Total damping of the felts.' },
47+ { key: 'feltGap', label: 'Felt thickness', kind: 'number', unit: 'mm', scale: 1000 },
48+ ],
49+ },
50+ {
51+ title: 'Damping (0 = off)',
52+ params: [
53+ { key: 't60Low', label: 'T60 low', kind: 'number', unit: 's', help: 'Decay time of bending vibration at the low frequency.' },
54+ { key: 'fLow', label: 'at', kind: 'number', unit: 'Hz' },
55+ { key: 't60High', label: 'T60 high', kind: 'number', unit: 's' },
56+ { key: 'fHigh', label: 'at', kind: 'number', unit: 'Hz' },
57+ ],
58+ },
59+ {
60+ title: 'Excitation',
61+ params: [
62+ {
63+ key: 'excitation',
64+ label: 'Type',
65+ kind: 'select',
66+ options: [
67+ ['stick', 'Drumstick (hard contact)'],
68+ ['mallet', 'Felt mallet'],
69+ ['strike', 'Strike (initial velocity)'],
70+ ['force', 'Force pulse'],
71+ ],
72+ },
73+ { key: 'excR', label: 'Position', kind: 'number', unit: '× radius', help: '0 = centre, 1 = rim.' },
74+ { key: 'excAngle', label: 'Azimuth', kind: 'number', unit: '°' },
75+ { key: 'hitAngle', label: 'Hit angle', kind: 'number', unit: '° from normal', help: '0 = straight into the surface; larger angles tilt the stroke toward the rim, like a crash stroke.' },
76+ { key: 'excWidth', label: 'Patch width', kind: 'number', unit: 'mm', scale: 1000, help: 'Gaussian width of the contact patch; 0 = a single bead. A real stick tip is smaller than a bead, but a single light bead makes the stick chatter on this coarse lattice, so the presets use about 1.5 bead spacings.' },
77+ { key: 'stickMass', label: 'Stick mass', kind: 'number', unit: 'g', scale: 1000, show: isExc('stick'), help: 'Effective mass of the stick at the tip.' },
78+ { key: 'stickSpeed', label: 'Stick speed', kind: 'number', unit: 'm/s', show: isExc('stick') },
79+ { key: 'stickK', label: 'Tip stiffness K', kind: 'number', unit: 'N/m^1.5', show: isExc('stick'), help: 'Hertz contact F = K δ^1.5. Wood on bronze with a 5 mm tip: about 5e8 to 1e9.' },
80+ { key: 'malletMass', label: 'Mallet mass', kind: 'number', unit: 'g', scale: 1000, show: isExc('mallet') },
81+ { key: 'malletSpeed', label: 'Mallet speed', kind: 'number', unit: 'm/s', show: isExc('mallet') },
82+ { key: 'malletK', label: 'Felt K', kind: 'number', unit: 'N/mᵖ', show: isExc('mallet') },
83+ { key: 'malletP', label: 'Felt exponent p', kind: 'number', show: isExc('mallet') },
84+ { key: 'friction', label: 'Friction', kind: 'number', show: isExc('stick', 'mallet'), help: 'Coulomb friction coefficient between tip and cymbal; matters for angled hits.' },
85+ { key: 'strikeVel', label: 'Velocity', kind: 'number', unit: 'm/s', show: isExc('strike') },
86+ { key: 'forcePeak', label: 'Peak force', kind: 'number', unit: 'N', show: isExc('force') },
87+ { key: 'forceMs', label: 'Duration', kind: 'number', unit: 'ms', show: isExc('force') },
88+ ],
89+ },
90+ {
91+ title: 'Microphone',
92+ params: [
93+ { key: 'micDistance', label: 'Distance', kind: 'number', unit: 'm' },
94+ { key: 'micElevation', label: 'Elevation', kind: 'number', unit: '°', help: 'Angle above the cymbal plane (90 = straight above, negative = below). Straight above, every mode with nodal diameters cancels.' },
95+ { key: 'micAzimuth', label: 'Azimuth', kind: 'number', unit: '°' },
96+ ],
97+ },
98+ {
99+ title: 'Simulation',
100+ params: [
101+ { key: 'duration', label: 'Duration', kind: 'number', unit: 's' },
102+ { key: 'fs', label: 'Sample rate', kind: 'select', options: [['44100', '44.1 kHz'], ['48000', '48 kHz'], ['96000', '96 kHz']] },
103+ { key: 'cfl', label: 'Step safety', kind: 'number', help: 'Fraction of the leapfrog stability limit dt = 2/ω_max.' },
104+ { key: 'frameRate', label: 'Replay frames', kind: 'number', unit: 'per s' },
105+ ],
106+ },
107+];
108+
109+export function mountBeadCymbal(root: HTMLElement) {
110+ const build = (p: CymbalParams): Built | null => {
111+ const d = deriveCymbal(p);
112+ if (d.n < 3 || !(p.thickness > 0)) return null;
113+ const sys = buildCymbal(p, d);
114+ const exc = buildCymbalExcitation(p, d);
115+ const mic = buildMic(p, d);
116+ return {
117+ sys,
118+ exc,
119+ pickup: mic.weights,
120+ pickupDelay: mic.delay,
121+ scene: {
122+ sys,
123+ exc,
124+ excKind: p.excitation,
125+ pickup: mic.weights,
126+ pickupLabel: 'microphone',
127+ felt: { radius: Math.max(p.feltRadius, d.mesh.a * 0.6), gap: p.feltGap, y: profileHeight(p, 0) },
128+ mic: {
129+ dir: micDirection(p),
130+ from: [0, profileHeight(p, 0), 0],
131+ label: `microphone direction (${fmtNum(p.micDistance)} m away, ${fmtNum(p.micElevation)}° up)`,
132+ },
133+ },
134+ };
135+ };
136+
137+ mountInstrumentPage<CymbalParams>(root, {
138+ title: 'Cymbal',
139+ wavName: 'bead-cymbal.wav',
140+ audioQuantity: { label: 'pressure at the microphone', unit: 'Pa' },
141+ presets: CYMBAL_PRESETS,
142+ sections: SECTIONS,
143+ setup(p) {
144+ const d = deriveCymbal(p);
145+ const R = p.diameter / 2;
146+ const rows: [string, string][] = [
147+ ['Beads', `${d.n} (${d.support.length} in the felts)`],
148+ ['Lattice spacing', `${fmtNum(d.mesh.a * 1000)} mm`],
149+ ['Mass', `${fmtNum(d.mass.reduce((s, m) => s + m, 0) * 1000)} g`],
150+ ['Bending stiffness D', `${fmtNum(d.D)} N·m`],
151+ ['In-plane wave speed', `${fmtNum(d.waveSpeedP)} m/s${p.soften > 1 ? ` (softened ${fmtNum(p.soften)}×)` : ''}`],
152+ ['Free flat plate (2,0) mode', `${fmtNum((5.253 / (2 * Math.PI * R * R)) * d.bendingCoef)} Hz (ideal, no felts)`],
153+ ['Bending resolved to', `~${fmtNum(bendingLimitHz(d))} Hz`],
154+ ];
155+ const warnings = [...d.warnings];
156+ const fmax = bendingLimitHz(d);
157+ if (fmax < 3000)
158+ warnings.push(
159+ `Bending waves are only resolved up to about ${fmtNum(fmax)} Hz: above that the lattice is too coarse and the sound is not cymbal-like. More rings help, at a steep cost.`,
160+ );
161+ if (Math.abs(p.micElevation) > 80) warnings.push('The microphone is almost straight above: modes with nodal diameters cancel there.');
162+ if (p.hitAngle >= 85) warnings.push('A hit angle near 90° barely touches the surface.');
163+ return { rows, warnings, built: build(p) };
164+ },
165+ summaryRows: (_p, _b, r) => [['Largest strain in the plate', `${fmtNum(r.maxStretch * 100)}%`]],
166+ extraCards: '<section class="card" data-r="out" hidden><h2>Cymbal (top view and cross-section)</h2><div data-r="membrane"></div></section>',
167+ mountExtras(q, timeline) {
168+ const view = new MembraneView(q('membrane'), timeline);
169+ return {
170+ show(p, _built, r) {
171+ const d = deriveCymbal(p);
172+ view.set(
173+ {
174+ pos: d.mesh.pos,
175+ radius: p.diameter / 2,
176+ spacing: d.mesh.a,
177+ strike: hitPoint(p),
178+ pickup: null,
179+ pickupNote:
180+ 'Sound: far-field microphone (Rayleigh integral: each bead’s normal acceleration × area, delayed by its distance; assumes a baffled plate, so sound bending around the rim is ignored).',
181+ sectionAngle: p.excR > 0.01 ? p.excAngle : 0,
182+ },
183+ r.frames,
184+ );
185+ },
186+ };
187+ },
188+ });
189+}
src/instruments/bead-cymbal/presets.tsadded+86−0View file
@@ -0,0 +1,86 @@
1+import type { CymbalParams } from './model';
2+
3+// 8-inch toy cymbal in B20 bronze (E ~ 100 GPa, 8600 kg/m^3), 0.8 mm thick.
4+const toy: CymbalParams = {
5+ diameter: 0.2,
6+ rings: 10,
7+ holeRadius: 0.006,
8+ bowHeight: 0,
9+ bellDiameter: 0.07,
10+ bellHeight: 0,
11+ thickness: 0.0008,
12+ youngs: 1.0e11,
13+ density: 8600,
14+ soften: 30,
15+ feltRadius: 0.02,
16+ feltK: 2e5,
17+ feltC: 5,
18+ feltGap: 0.006,
19+ t60Low: 4,
20+ fLow: 200,
21+ t60High: 1,
22+ fHigh: 2000,
23+ excitation: 'stick',
24+ excR: 0.85,
25+ excAngle: 0,
26+ hitAngle: 0,
27+ excWidth: 0.015, // wider than a stick tip: one light bead alone makes the stick chatter
28+ stickMass: 0.03,
29+ stickSpeed: 4,
30+ stickK: 5e8,
31+ malletMass: 0.03,
32+ malletSpeed: 2,
33+ malletK: 1e7,
34+ malletP: 2.5,
35+ friction: 0.3,
36+ strikeVel: 1,
37+ forcePeak: 50,
38+ forceMs: 0.5,
39+ micDistance: 1,
40+ micElevation: 30,
41+ micAzimuth: 30,
42+ duration: 1,
43+ fs: 48000,
44+ cfl: 0.5,
45+ frameRate: 2000,
46+};
47+
48+export interface CymbalPreset {
49+ name: string;
50+ description: string;
51+ params: CymbalParams;
52+}
53+
54+export const CYMBAL_PRESETS: CymbalPreset[] = [
55+ {
56+ name: 'Flat 8" toy, stick (fast)',
57+ description:
58+ 'Flat 0.8 mm bronze plate, 20 cm, ~330 beads, on felts. In-plane stiffness softened 30× for speed (bending is real). Stick hit near the edge, spread over a 15 mm patch (a single bead would make the stick chatter). Bending is resolved to about 1.3 kHz.',
59+ params: toy,
60+ },
61+ {
62+ name: 'Flat 8" toy, real bronze (slow)',
63+ description: 'Same plate with real in-plane stiffness: about 5× more time steps.',
64+ params: { ...toy, soften: 1 },
65+ },
66+ {
67+ name: 'Bow and bell',
68+ description: '12 mm bow (dome) and a 7 cm bell, 12 mm high. Curvature couples bending to in-plane stretching, so softening matters more here.',
69+ params: { ...toy, bowHeight: 0.012, bellHeight: 0.012 },
70+ },
71+ {
72+ name: 'Glancing crash stroke',
73+ description: 'Stick at 6 m/s, 60° from the normal toward the rim, with friction: pushes the plate sideways as well as down.',
74+ params: { ...toy, bowHeight: 0.012, bellHeight: 0.012, stickSpeed: 6, hitAngle: 60 },
75+ },
76+ {
77+ name: 'Felt mallet',
78+ description: 'Soft 30 g mallet at 2 m/s, longer contact, darker sound.',
79+ params: { ...toy, excitation: 'mallet', excR: 0.7, excWidth: 0.008 },
80+ },
81+ {
82+ name: 'Coarse quick look',
83+ description: '6 rings (~120 beads): fast, but bending is only resolved to a few hundred hertz.',
84+ params: { ...toy, rings: 6 },
85+ },
86+];
src/instruments/bead-drum/model.tsmodified+1−1View file
@@ -242,7 +242,7 @@ export function buildDrumExcitation(p: DrumParams, d: DrumDerived): Excitation {
242242 return {
243243 u0,
244244 v0,
245- mallet: { weights, dir: [0, -1, 0], mass: p.malletMass, K: p.feltK, p: p.feltP, v0: p.malletSpeed },
245+ mallet: { weights, normal: [0, -1, 0], vel0: [0, -p.malletSpeed, 0], mass: p.malletMass, K: p.feltK, p: p.feltP, mu: 0 },
246246 };
247247 }
248248
src/instruments/bead-drum/page.tsmodified+3−53View file
@@ -1,5 +1,5 @@
11 import type { Section } from '../../ui/controls';
2-import { LinePlot, css, fmtNum } from '../../ui/plot';
2+import { fmtNum } from '../../ui/plot';
33 import { MembraneView } from '../../ui/membraneview';
44 import { mountInstrumentPage, type Built } from '../../ui/instrumentPage';
55 import {
@@ -87,21 +87,6 @@ const SECTIONS: Section<DrumParams>[] = [
8787 },
8888 ];
8989
90-// separate contact intervals in the contact-force record
91-function contacts(cf: Float32Array, fs: number): { start: number; ms: number; peak: number }[] {
92- const out: { start: number; ms: number; peak: number }[] = [];
93- let i = 0;
94- while (i < cf.length) {
95- if (cf[i] > 0) {
96- const a = i;
97- let peak = 0;
98- while (i < cf.length && cf[i] > 0) peak = Math.max(peak, cf[i++]);
99- out.push({ start: a / fs, ms: ((i - a) / fs) * 1000, peak });
100- } else i++;
101- }
102- return out;
103-}
104-
10590 export function mountBeadDrum(root: HTMLElement) {
10691 const build = (p: DrumParams): Built | null => {
10792 const d = deriveDrum(p);
@@ -145,29 +130,11 @@ export function mountBeadDrum(root: HTMLElement) {
145130 if (d.n > 3000) warnings.push(`${d.n} beads: expect a slow run.`);
146131 return { rows, warnings, built: build(p) };
147132 },
148- summaryRows(p, _b, r) {
149- const rows: [string, string][] = [['Largest spring strain', `${fmtNum(r.maxStretch * 100)}% (pre-stretch ${fmtNum(deriveDrum(p).prestretch * 100)}%)`]];
150- if (r.contactForce) {
151- const c = contacts(r.contactForce, r.fs);
152- rows.push([
153- 'Mallet contact',
154- c.length
155- ? `${fmtNum(c[0].ms)} ms, peak ${fmtNum(c[0].peak)} N${c.length > 1 ? ` (+${c.length - 1} rebound${c.length > 2 ? 's' : ''})` : ''}`
156- : 'no contact',
157- ]);
158- }
159- return rows;
160- },
133+ summaryRows: (p, _b, r) => [['Largest spring strain', `${fmtNum(r.maxStretch * 100)}% (pre-stretch ${fmtNum(deriveDrum(p).prestretch * 100)}%)`]],
161134 extraCards: `
162- <section class="card" data-r="out" hidden><h2>Drumhead (top view and cross-section)</h2><div data-r="membrane"></div></section>
163- <section class="card" data-r="out-contact" hidden>
164- <h2>Mallet contact force</h2>
165- <div class="canvas-wrap short"><canvas data-r="contact"></canvas></div>
166- <p class="note" data-r="contact-note"></p>
167- </section>`,
135+ <section class="card" data-r="out" hidden><h2>Drumhead (top view and cross-section)</h2><div data-r="membrane"></div></section>`,
168136 mountExtras(q, timeline) {
169137 const view = new MembraneView(q('membrane'), timeline);
170- const contactPlot = new LinePlot(q('contact'));
171138 return {
172139 show(p, _built, r) {
173140 const d = deriveDrum(p);
@@ -182,23 +149,6 @@ export function mountBeadDrum(root: HTMLElement) {
182149 },
183150 r.frames,
184151 );
185- const card = q('out-contact');
186- card.hidden = !r.contactForce;
187- if (!r.contactForce) return;
188- const cf = r.contactForce;
189- const c = contacts(cf, r.fs);
190- const end = c.length ? c[c.length - 1].start + c[c.length - 1].ms / 1000 : 0.01;
191- const nShow = Math.min(cf.length, Math.ceil((end * 1.3 + 0.002) * r.fs));
192- const t = new Float64Array(nShow);
193- for (let i = 0; i < nShow; i++) t[i] = (i / r.fs) * 1000;
194- contactPlot.set([{ x: t, y: cf.subarray(0, nShow), color: css('--series-2'), label: 'force (N)' }], {
195- xLabel: 'time (ms)',
196- yLabel: 'contact force (N)',
197- yRange: [0, Math.max(...cf.subarray(0, nShow)) * 1.08 || 1],
198- });
199- q('contact-note').textContent = c.length
200- ? `${c.length} contact${c.length > 1 ? 's' : ''}: ` + c.map((k) => `${fmtNum(k.ms)} ms from ${fmtNum(k.start * 1000)} ms, peak ${fmtNum(k.peak)} N`).join('; ') + `. During contact the head takes energy from the mallet and gives much of it back as the mallet rebounds, so the energy plot spikes first; the mallet then leaves at constant speed (no gravity).`
201- : 'The mallet never touched the head.';
202152 },
203153 };
204154 },
src/instruments/registry.tsmodified+7−1View file
@@ -1,5 +1,6 @@
11 import { mountBeadString } from './bead-string/page';
22 import { mountBeadDrum } from './bead-drum/page';
3+import { mountBeadCymbal } from './bead-cymbal/page';
34
45 export interface Instrument {
56 id: string;
@@ -21,6 +22,11 @@ export const INSTRUMENTS: Instrument[] = [
2122 blurb: 'A round drumhead: a tensioned hex net of beads in a clamped rim. Hit it with a felt mallet, strike, push or press it.',
2223 mount: mountBeadDrum,
2324 },
25+ {
26+ id: 'bead-cymbal',
27+ name: 'Cymbal',
28+ blurb: 'A thin bronze plate on felts at its centre, flat or with bow and bell. Hit it with a stick or mallet, at any angle; listen through a microphone.',
29+ mount: mountBeadCymbal,
30+ },
2431 { id: 'bead-bar', name: 'Bar', blurb: 'A stiff chain with bending springs, free ends: xylophone-like overtones.' },
25- { id: 'bead-plate', name: 'Plate', blurb: 'A 2D net with bending, free or clamped edge.' },
2632 ];
src/sim/hinge.tsadded+62−0View file
@@ -0,0 +1,62 @@
1+// Dihedral angle across an edge shared by two triangles, and its gradient
2+// (Bridson, Marino & Fedkiw 2003). Edge (e0, e1); wings w0 and w1 on either
3+// side. theta = 0 when the two triangles are coplanar; the sign tells which way
4+// the hinge folds. Positions are read from X (3 per node).
5+
6+export const hingeOut = new Float64Array(13); // [theta, g_w0(3), g_w1(3), g_e0(3), g_e1(3)]
7+
8+export function hingeAngleGrad(X: Float64Array, w0: number, w1: number, e0: number, e1: number, out = hingeOut): number {
9+ const a = 3 * w0;
10+ const b = 3 * w1;
11+ const c = 3 * e0;
12+ const d = 3 * e1;
13+ // E = x_e1 - x_e0
14+ const Ex = X[d] - X[c];
15+ const Ey = X[d + 1] - X[c + 1];
16+ const Ez = X[d + 2] - X[c + 2];
17+ // N1 = (x_w0 - x_e0) x (x_w0 - x_e1)
18+ const p1x = X[a] - X[c], p1y = X[a + 1] - X[c + 1], p1z = X[a + 2] - X[c + 2];
19+ const q1x = X[a] - X[d], q1y = X[a + 1] - X[d + 1], q1z = X[a + 2] - X[d + 2];
20+ const N1x = p1y * q1z - p1z * q1y;
21+ const N1y = p1z * q1x - p1x * q1z;
22+ const N1z = p1x * q1y - p1y * q1x;
23+ // N2 = (x_w1 - x_e1) x (x_w1 - x_e0)
24+ const p2x = X[b] - X[d], p2y = X[b + 1] - X[d + 1], p2z = X[b + 2] - X[d + 2];
25+ const q2x = X[b] - X[c], q2y = X[b + 1] - X[c + 1], q2z = X[b + 2] - X[c + 2];
26+ const N2x = p2y * q2z - p2z * q2y;
27+ const N2y = p2z * q2x - p2x * q2z;
28+ const N2z = p2x * q2y - p2y * q2x;
29+ const n1sq = N1x * N1x + N1y * N1y + N1z * N1z;
30+ const n2sq = N2x * N2x + N2y * N2y + N2z * N2z;
31+ const eLen = Math.sqrt(Ex * Ex + Ey * Ey + Ez * Ez);
32+ const n1 = Math.sqrt(n1sq);
33+ const n2 = Math.sqrt(n2sq);
34+ // angle between unit normals, signed about the edge direction
35+ const cos = (N1x * N2x + N1y * N2y + N1z * N2z) / (n1 * n2);
36+ const cx = N1y * N2z - N1z * N2y;
37+ const cy = N1z * N2x - N1x * N2z;
38+ const cz = N1x * N2y - N1y * N2x;
39+ const sin = (cx * Ex + cy * Ey + cz * Ez) / (n1 * n2 * eLen);
40+ const theta = Math.atan2(sin, cos);
41+ // gradients
42+ const s1 = eLen / n1sq;
43+ const s2 = eLen / n2sq;
44+ out[0] = theta;
45+ out[1] = -s1 * N1x;
46+ out[2] = -s1 * N1y;
47+ out[3] = -s1 * N1z;
48+ out[4] = -s2 * N2x;
49+ out[5] = -s2 * N2y;
50+ out[6] = -s2 * N2z;
51+ const t1 = (q1x * Ex + q1y * Ey + q1z * Ez) / eLen / n1sq; // (x_w0 - x_e1).E / |E| / |N1|^2
52+ const t2 = (p2x * Ex + p2y * Ey + p2z * Ez) / eLen / n2sq; // (x_w1 - x_e1).E / |E| / |N2|^2
53+ const r1 = (p1x * Ex + p1y * Ey + p1z * Ez) / eLen / n1sq; // (x_w0 - x_e0).E
54+ const r2 = (q2x * Ex + q2y * Ey + q2z * Ez) / eLen / n2sq; // (x_w1 - x_e0).E
55+ out[7] = -(t1 * N1x + t2 * N2x);
56+ out[8] = -(t1 * N1y + t2 * N2y);
57+ out[9] = -(t1 * N1z + t2 * N2z);
58+ out[10] = r1 * N1x + r2 * N2x;
59+ out[11] = r1 * N1y + r2 * N2y;
60+ out[12] = r1 * N1z + r2 * N2z;
61+ return theta;
62+}
src/sim/solver.tsmodified+243−70View file
@@ -1,32 +1,56 @@
1-import type { BeadSystem, Excitation, Mallet, SimOptions, SimResult } from './types';
1+import type { BeadSystem, Bending, Excitation, Mallet, SimOptions, SimResult } from './types';
2+import { hingeAngleGrad } from './hinge';
3+
4+// Normal speed of a mallet into the surface, and its linearized contact
5+// stiffness at the largest compression a rigid surface would see.
6+function contactStiffness(m: Mallet): number {
7+ const vn = m.vel0[0] * m.normal[0] + m.vel0[1] * m.normal[1] + m.vel0[2] * m.normal[2];
8+ if (!(vn > 0)) return 0;
9+ const dmax = (((m.p + 1) * m.mass * vn * vn) / (2 * m.K)) ** (1 / (m.p + 1));
10+ return m.p * m.K * dmax ** (m.p - 1);
11+}
12+
13+function nodePositions(sys: BeadSystem): Float64Array {
14+ const X = new Float64Array(3 * (sys.n + sys.anchors.length));
15+ X.set(sys.rest);
16+ sys.anchors.forEach((a, k) => X.set(a.pos, 3 * (sys.n + k)));
17+ return X;
18+}
219
320 // Upper bound on the highest angular frequency of the linearized system
421 // (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
522 // its axis and T/l across it, so each spring contributes 2 max(k, |T|/l) / m.
23+// A hinge contributes k g g^T (g = gradient of its angle); a mallet a rank-one
24+// contact stiffness between it and the patch.
625 export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
726 const { n, si, sj, sk, sL0 } = sys;
827 const rowSum = new Float64Array(n);
9- const p = (node: number, c: number) =>
10- node < n ? sys.rest[3 * node + c] : sys.anchors[node - n].pos[c];
28+ const X = nodePositions(sys);
1129 for (let s = 0; s < si.length; s++) {
1230 const a = si[s];
1331 const b = sj[s];
14- let l2 = 0;
15- for (let c = 0; c < 3; c++) l2 += (p(b, c) - p(a, c)) ** 2;
16- const l = Math.sqrt(l2);
32+ const l = Math.hypot(X[3 * b] - X[3 * a], X[3 * b + 1] - X[3 * a + 1], X[3 * b + 2] - X[3 * a + 2]);
1733 const T = sk[s] * (l - sL0[s]);
1834 const kmax = Math.max(sk[s], Math.abs(T) / l);
1935 if (a < n) rowSum[a] += 2 * kmax;
2036 if (b < n) rowSum[b] += 2 * kmax;
2137 }
22- // a mallet adds a rank-one contact stiffness k w w^T between it and the beads,
23- // linearized at the largest compression a rigid surface would see
24- let kc = 0;
25- if (mallet) {
26- const { K, p, mass, v0 } = mallet;
27- const dmax = (((p + 1) * mass * v0 * v0) / (2 * K)) ** (1 / (p + 1));
28- kc = p * K * dmax ** (p - 1);
38+ const h = sys.bending;
39+ if (h) {
40+ const kh = hingeStiffness(h);
41+ const g = new Float64Array(13);
42+ for (let e = 0; e < kh.length; e++) {
43+ const nodes = [h.w0[e], h.w1[e], h.e0[e], h.e1[e]];
44+ hingeAngleGrad(X, nodes[0], nodes[1], nodes[2], nodes[3], g);
45+ let l1 = 0;
46+ for (let q = 1; q < 13; q++) l1 += Math.abs(g[q]);
47+ nodes.forEach((node, k) => {
48+ const inf = Math.max(Math.abs(g[1 + 3 * k]), Math.abs(g[2 + 3 * k]), Math.abs(g[3 + 3 * k]));
49+ if (node < n) rowSum[node] += kh[e] * inf * l1;
50+ });
51+ }
2952 }
53+ const kc = mallet ? contactStiffness(mallet) : 0;
3054 let w2 = mallet ? kc / mallet.mass : 0;
3155 for (let i = 0; i < n; i++) {
3256 if (sys.pinned[i]) continue;
@@ -36,21 +60,53 @@ export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
3660 return Math.sqrt(w2);
3761 }
3862
63+// per hinge: sum of |M| entries of its rows over the triangles that hold it
64+function hingeStiffness(b: Bending): Float64Array {
65+ const k = new Float64Array(b.w0.length);
66+ for (let t = 0; t < b.tri.length / 3; t++)
67+ for (let a = 0; a < 3; a++) {
68+ const e = b.tri[3 * t + a];
69+ if (e < 0) continue;
70+ for (let c = 0; c < 3; c++) k[e] += Math.abs(b.M[9 * t + 3 * a + c]);
71+ }
72+ return k;
73+}
74+
75+// bending energy at positions X, optionally only over triangles with keep(t)
76+export function bendingEnergy(b: Bending, X: Float64Array, keep?: (t: number) => boolean): number {
77+ const g = new Float64Array(13);
78+ const phi = Float64Array.from(b.w0, (_, e) => b.sign[e] * (hingeAngleGrad(X, b.w0[e], b.w1[e], b.e0[e], b.e1[e], g) - b.theta0[e]));
79+ let E = 0;
80+ for (let t = 0; t < b.tri.length / 3; t++) {
81+ if (keep && !keep(t)) continue;
82+ for (let a = 0; a < 3; a++)
83+ for (let c = 0; c < 3; c++) {
84+ const ea = b.tri[3 * t + a];
85+ const ec = b.tri[3 * t + c];
86+ if (ea >= 0 && ec >= 0) E += 0.5 * phi[ea] * b.M[9 * t + 3 * a + c] * phi[ec];
87+ }
88+ }
89+ return E;
90+}
91+
3992 export type Progress = (fraction: number) => void;
4093
41-// Velocity Verlet (kick-drift-kick), explicit. Spring forces are exact
42-// (geometric nonlinearity kept); damping uses the half-step velocity in the
43-// kick, and the full-step velocity whenever an audio sample is recorded.
44-export function simulate(
45- sys: BeadSystem,
46- exc: Excitation,
47- opt: SimOptions,
48- progress?: Progress,
49-): SimResult {
94+// Velocity Verlet (kick-drift-kick), explicit. Spring and hinge forces are
95+// exact (geometric nonlinearity kept); damping and friction use the half-step
96+// velocity in the kick, and the full-step velocity whenever an audio sample is
97+// recorded.
98+export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, progress?: Progress): SimResult {
5099 const t0 = performance.now();
51100 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
52101 const nn = n + anchors.length;
53102 const ns = si.length;
103+ const bend = sys.bending;
104+ const nh = bend ? bend.w0.length : 0;
105+ const nt = bend ? bend.tri.length / 3 : 0;
106+ const hG = new Float64Array(12 * nh); // gradients of phi
107+ const hPhi = new Float64Array(nh);
108+ const hRate = new Float64Array(nh);
109+ const hMom = new Float64Array(nh);
54110
55111 const mallet = exc.mallet;
56112 const omegaMax = omegaMaxBound(sys, mallet);
@@ -92,53 +148,92 @@ export function simulate(
92148 const force = exc.force;
93149 const fShape = force?.shape;
94150 const pulse = (t: number) =>
95- force && t >= 0 && t <= force.tau
96- ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau))
97- : 0;
151+ force && t >= 0 && t <= force.tau ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau)) : 0;
98152
99- // mallet state: position s and velocity along dir; contact force and its potential
100- let ms = 0;
101- let mv = mallet ? mallet.v0 : 0;
102- let ma = 0;
153+ // mallet state: 3D position q, velocity qv, acceleration qa
154+ const q = new Float64Array(3);
155+ const qv = new Float64Array(mallet ? mallet.vel0 : [0, 0, 0]);
156+ const qa = new Float64Array(3);
103157 let fContact = 0;
104158 let eContact = 0;
105- const mw = mallet?.weights;
106- const [mdx, mdy, mdz] = mallet?.dir ?? [0, 0, 0];
159+ let pFriction = 0;
160+ const patch: number[] = [];
161+ let mEff = 1;
162+ if (mallet) {
163+ let inv = 0;
164+ for (let i = 0; i < n; i++)
165+ if (mallet.weights[i] > 0) {
166+ patch.push(i);
167+ inv += (mallet.weights[i] * mallet.weights[i]) / mass[i];
168+ }
169+ mEff = 1 / (1 / mallet.mass + inv);
170+ }
171+ const [nx, ny, nz] = mallet?.normal ?? [0, 0, 0];
107172 const contact = () => {
108173 fContact = 0;
109174 eContact = 0;
175+ pFriction = 0;
176+ qa.fill(0);
110177 if (!mallet) return;
111- let h = 0;
112- for (let i = 0; i < n; i++) {
113- const wi = mw![i];
114- if (wi === 0) continue;
178+ const w = mallet.weights;
179+ let ux = 0, uy = 0, uz = 0, vx = 0, vy = 0, vz = 0;
180+ for (const i of patch) {
115181 const k = 3 * i;
116- h += wi * ((X[k] - sys.rest[k]) * mdx + (X[k + 1] - sys.rest[k + 1]) * mdy + (X[k + 2] - sys.rest[k + 2]) * mdz);
182+ ux += w[i] * (X[k] - sys.rest[k]);
183+ uy += w[i] * (X[k + 1] - sys.rest[k + 1]);
184+ uz += w[i] * (X[k + 2] - sys.rest[k + 2]);
185+ vx += w[i] * V[k];
186+ vy += w[i] * V[k + 1];
187+ vz += w[i] * V[k + 2];
117188 }
118- const d = ms - h;
119- if (d > 0) {
120- const dp = mallet.K * d ** mallet.p;
121- fContact = dp;
122- eContact = (dp * d) / (mallet.p + 1);
123- for (let i = 0; i < n; i++) {
124- const wi = mw![i];
125- if (wi === 0) continue;
126- F[3 * i] += fContact * wi * mdx;
127- F[3 * i + 1] += fContact * wi * mdy;
128- F[3 * i + 2] += fContact * wi * mdz;
129- }
189+ const d = (q[0] - ux) * nx + (q[1] - uy) * ny + (q[2] - uz) * nz;
190+ if (!(d > 0)) return;
191+ if (mallet.reach !== undefined) {
192+ const ox = q[0] - ux - d * nx, oy = q[1] - uy - d * ny, oz = q[2] - uz - d * nz;
193+ if (ox * ox + oy * oy + oz * oz > mallet.reach * mallet.reach) return;
194+ }
195+ const fn = mallet.K * d ** mallet.p;
196+ fContact = fn;
197+ eContact = (fn * d) / (mallet.p + 1);
198+ // force on the mallet: -fn normal + friction against tangential slip
199+ let fx = -fn * nx;
200+ let fy = -fn * ny;
201+ let fz = -fn * nz;
202+ if (mallet.mu > 0) {
203+ let rx = qv[0] - vx, ry = qv[1] - vy, rz = qv[2] - vz;
204+ const rn = rx * nx + ry * ny + rz * nz;
205+ rx -= rn * nx;
206+ ry -= rn * ny;
207+ rz -= rn * nz;
208+ const slip = Math.hypot(rx, ry, rz);
209+ // regularized Coulomb; the slip scale keeps the explicit step stable
210+ const eps = Math.max(1e-3, (4 * mallet.mu * fn * dt) / mEff);
211+ const g = (mallet.mu * fn) / Math.sqrt(slip * slip + eps * eps);
212+ fx -= g * rx;
213+ fy -= g * ry;
214+ fz -= g * rz;
215+ pFriction = g * slip * slip;
216+ }
217+ qa[0] = fx / mallet.mass;
218+ qa[1] = fy / mallet.mass;
219+ qa[2] = fz / mallet.mass;
220+ for (const i of patch) {
221+ F[3 * i] -= w[i] * fx;
222+ F[3 * i + 1] -= w[i] * fy;
223+ F[3 * i + 2] -= w[i] * fz;
130224 }
131- ma = -fContact / mallet.mass;
132225 };
133226
134- // Fills F with total force on beads; returns damping power (W, >= 0)
135- // and external power (W), both at velocity V.
227+ const hOut = new Float64Array(13);
228+ let eHinge = 0;
229+ // Fills F with total force on beads; sets damping power (W, >= 0) and
230+ // external power (W), both at velocity V.
136231 let pDamp = 0;
137232 let pExt = 0;
138233 const computeForces = (t: number) => {
139234 F.fill(0);
140235 contact();
141- pDamp = 0;
236+ pDamp = pFriction;
142237 for (let s = 0; s < ns; s++) {
143238 const a = 3 * si[s];
144239 const b = 3 * sj[s];
@@ -171,6 +266,67 @@ export function simulate(
171266 F[b + 2] -= fz;
172267 }
173268 }
269+ eHinge = 0;
270+ if (bend) {
271+ const s1 = bend.sigma1;
272+ for (let e = 0; e < nh; e++) {
273+ const th = hingeAngleGrad(X, bend.w0[e], bend.w1[e], bend.e0[e], bend.e1[e], hOut);
274+ const sg = bend.sign[e];
275+ hPhi[e] = sg * (th - bend.theta0[e]);
276+ const o = 12 * e;
277+ for (let q = 0; q < 12; q++) hG[o + q] = sg * hOut[1 + q];
278+ if (s1 !== 0) {
279+ const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
280+ hRate[e] =
281+ hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
282+ hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
283+ hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
284+ hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
285+ }
286+ }
287+ hMom.fill(0);
288+ const tri = bend.tri;
289+ const Mt = bend.M;
290+ for (let t = 0; t < nt; t++) {
291+ for (let a = 0; a < 3; a++) {
292+ const ea = tri[3 * t + a];
293+ if (ea < 0) continue;
294+ let m = 0;
295+ let r = 0;
296+ for (let c = 0; c < 3; c++) {
297+ const ec = tri[3 * t + c];
298+ if (ec < 0) continue;
299+ const M = Mt[9 * t + 3 * a + c];
300+ m += M * hPhi[ec];
301+ r += M * hRate[ec];
302+ }
303+ eHinge += 0.5 * hPhi[ea] * m;
304+ if (s1 !== 0) {
305+ pDamp += 2 * s1 * hRate[ea] * r;
306+ m += 2 * s1 * r;
307+ }
308+ hMom[ea] += m;
309+ }
310+ }
311+ for (let e = 0; e < nh; e++) {
312+ const m = hMom[e];
313+ if (m === 0) continue;
314+ const o = 12 * e;
315+ const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
316+ F[w0] -= m * hG[o];
317+ F[w0 + 1] -= m * hG[o + 1];
318+ F[w0 + 2] -= m * hG[o + 2];
319+ F[w1] -= m * hG[o + 3];
320+ F[w1 + 1] -= m * hG[o + 4];
321+ F[w1 + 2] -= m * hG[o + 5];
322+ F[e0] -= m * hG[o + 6];
323+ F[e0 + 1] -= m * hG[o + 7];
324+ F[e0 + 2] -= m * hG[o + 8];
325+ F[e1] -= m * hG[o + 9];
326+ F[e1 + 1] -= m * hG[o + 10];
327+ F[e1 + 2] -= m * hG[o + 11];
328+ }
329+ }
174330 const s0 = sys.sigma0;
175331 const p = fShape ? pulse(t) : 0;
176332 pExt = 0;
@@ -203,8 +359,12 @@ export function simulate(
203359 }
204360 };
205361
206- const damped = sys.sigma0 !== 0 || sc.some((c) => c !== 0);
362+ const velocityForces =
363+ sys.sigma0 !== 0 || sc.some((c) => c !== 0) || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
207364 let maxStretch = 0;
365+ const malletEnergy = () => (mallet ? 0.5 * mallet.mass * (qv[0] ** 2 + qv[1] ** 2 + qv[2] ** 2) + eContact : 0);
366+ // total mechanical energy at the current state; computeForces must have run
367+ // at this state (it sets the hinge and contact energies)
208368 const energyNow = () => {
209369 let e = 0;
210370 for (let i = 0; i < n; i++) {
@@ -217,10 +377,9 @@ export function simulate(
217377 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
218378 const d = l - sL0[s];
219379 e += 0.5 * sk[s] * d * d;
220- maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
380+ if (sj[s] < n && si[s] < n) maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
221381 }
222- if (mallet) e += 0.5 * mallet.mass * mv * mv + eContact;
223- return e;
382+ return e + eHinge + malletEnergy();
224383 };
225384
226385 const audio = new Float32Array(nSamples);
@@ -229,8 +388,9 @@ export function simulate(
229388 const energy = new Float64Array(nFrames);
230389 const balance = new Float64Array(nFrames);
231390 const contactForce = mallet ? new Float32Array(nSamples) : undefined;
232- const malletPos = mallet ? new Float32Array(nFrames) : undefined;
391+ const malletPos = mallet ? new Float32Array(3 * nFrames) : undefined;
233392 const w = opt.pickup;
393+ const delay = opt.pickupDelay;
234394
235395 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
236396 let Eeq = 0;
@@ -239,29 +399,41 @@ export function simulate(
239399 const l = Math.hypot(pa(sj[s], 0) - pa(si[s], 0), pa(sj[s], 1) - pa(si[s], 1), pa(sj[s], 2) - pa(si[s], 2));
240400 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
241401 }
242- computeForces(0); // contact energy at t = 0
402+ computeForces(0);
243403 const E0 = energyNow();
244404 let dissipated = 0;
245405 let workIn = 0;
246406 let t = 0;
247- computeForces(0);
248407 let frame = 0;
249408 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
250409
251410 for (let smp = 0; smp < nSamples; smp++) {
252411 // record at t = smp / fs; A holds forces at full-step velocity here
253- let acc = 0;
254- for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
255- audio[smp] = acc;
412+ if (delay) {
413+ for (let i = 0; i < n; i++) {
414+ const k = 3 * i;
415+ const c = w[k] * A[k] + w[k + 1] * A[k + 1] + w[k + 2] * A[k + 2];
416+ if (c === 0) continue;
417+ const pos = smp + delay[i];
418+ const i0 = Math.floor(pos);
419+ const fr = pos - i0;
420+ if (i0 < nSamples) audio[i0] += (1 - fr) * c;
421+ if (i0 + 1 < nSamples) audio[i0 + 1] += fr * c;
422+ }
423+ } else {
424+ let acc = 0;
425+ for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
426+ audio[smp] = acc;
427+ }
256428 if (contactForce) contactForce[smp] = fContact;
257429 if (smp % samplesPerFrame === 0 && frame < nFrames) {
258- if (malletPos) malletPos[frame] = ms;
430+ if (malletPos) malletPos.set(q, 3 * frame);
259431 const off = frame * 3 * n;
260432 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
261433 frameTimes[frame] = t;
262434 const eTotal = energyNow();
263435 // the plot shows the instrument's own energy; the balance includes the mallet
264- energy[frame] = eTotal - Eeq - (mallet ? 0.5 * mallet.mass * mv * mv + eContact : 0);
436+ energy[frame] = eTotal - Eeq - malletEnergy();
265437 balance[frame] = eTotal + dissipated - workIn - E0;
266438 frame++;
267439 }
@@ -273,10 +445,11 @@ export function simulate(
273445 V[k] += 0.5 * dt * A[k];
274446 X[k] += dt * V[k];
275447 }
276- if (mallet) {
277- mv += 0.5 * dt * ma;
278- ms += dt * mv;
279- }
448+ if (mallet)
449+ for (let c = 0; c < 3; c++) {
450+ qv[c] += 0.5 * dt * qa[c];
451+ q[c] += dt * qv[c];
452+ }
280453 syncAnchors();
281454 t += dt;
282455 // forces at new position with half-step velocity; powers use it too
@@ -284,11 +457,11 @@ export function simulate(
284457 dissipated += dt * pDamp;
285458 workIn += dt * pExt;
286459 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
287- if (mallet) mv += 0.5 * dt * ma;
460+ if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];
288461 syncAnchors();
289462 }
290463 // at a sample time: recompute with the full-step velocity
291- if (damped) computeForces(t);
464+ if (velocityForces) computeForces(t);
292465 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
293466 }
294467
src/sim/types.tsmodified+34−7View file
@@ -21,6 +21,27 @@ export interface BeadSystem {
2121 sL0: Float64Array; // rest length (m)
2222 sc: Float64Array; // dashpot on relative velocity, all 3 components (N s/m)
2323 sigma0: number; // mass-proportional damping: force -2 sigma0 m v
24+ bending?: Bending;
25+}
26+
27+// Plate bending (Poisson ratio built into M). Hinge e measures the dihedral
28+// angle across edge (e0, e1) shared by triangles (e0, e1, w0) and (e1, e0, w1);
29+// phi_e = sign_e (theta_e - theta0_e) is its change, positive when the surface
30+// folds toward its upward normal. Triangle t holds up to three hinges
31+// (-1 = free edge) and has energy 1/2 phi^T M phi; M is the plate energy
32+// A D/2 [(tr S)^2 - 2 (1 - nu) det S] of the shape operator
33+// S = sum_e phi_e l_e / (2A) m_e m_e^T, condensed over free edges (no moment there).
34+// Damping: the moments use phi + 2 sigma1 dphi/dt (Rayleigh, sigma1 omega^2).
35+export interface Bending {
36+ w0: Int32Array;
37+ w1: Int32Array;
38+ e0: Int32Array;
39+ e1: Int32Array;
40+ sign: Float64Array;
41+ theta0: Float64Array;
42+ tri: Int32Array; // 3 per triangle
43+ M: Float64Array; // 9 per triangle
44+ sigma1: number;
2445 }
2546
2647 export interface Excitation {
@@ -28,19 +49,23 @@ export interface Excitation {
2849 v0: Float64Array; // 3n initial velocity
2950 // External force F(t) = shape * F0 * 0.5 (1 - cos(2 pi t / tau)), 0 <= t <= tau
3051 force?: { shape: Float64Array; F0: number; tau: number };
31- // A mallet: point mass moving along dir, touching the beads through a felt
32- // spring with force K * delta^p, delta = s - sum_i w_i (x_i - rest_i) . dir.
33- // Starts at s = 0 (just touching the undeformed surface) with speed v0 along dir.
52+ // A mallet or stick: point mass q in 3D touching a patch of beads.
53+ // Normal force K delta^p with delta = (q - sum_i w_i (x_i - rest_i)) . normal;
54+ // Coulomb friction mu * normal force against tangential slip.
55+ // Starts at q = 0 (just touching the undeformed patch) with velocity vel0.
3456 mallet?: Mallet;
3557 }
3658
3759 export interface Mallet {
3860 weights: Float64Array; // n, sums to 1 (contact patch)
39- dir: [number, number, number]; // unit
61+ normal: [number, number, number]; // unit, pointing into the surface
62+ vel0: [number, number, number]; // m/s; its component along normal must be > 0
4063 mass: number; // kg
4164 K: number; // N / m^p
4265 p: number;
43- v0: number; // m/s, > 0 toward the surface
66+ mu: number; // friction coefficient (0 = frictionless)
67+ // no contact once the mallet has slid this far sideways from the patch (m); default unlimited
68+ reach?: number;
4469 }
4570
4671 export interface SimOptions {
@@ -49,6 +74,8 @@ export interface SimOptions {
4974 cfl: number; // fraction of the leapfrog stability limit 2/omega_max
5075 frameRate: number; // displacement frames per second
5176 pickup: Float64Array; // 3n weights; audio = sum w . acceleration
77+ // optional per-bead delay in samples (a microphone): bead i's term arrives delay[i] later
78+ pickupDelay?: Float64Array;
5279 }
5380
5481 export interface SimResult {
@@ -62,7 +89,7 @@ export interface SimResult {
6289 stepsPerSample: number;
6390 omegaMax: number;
6491 wallSeconds: number;
65- maxStretch: number; // largest |l - L0| / L0 seen, at frame times
92+ maxStretch: number; // largest |l - L0| / L0 of bead-to-bead springs, at frame times
6693 contactForce?: Float32Array; // mallet contact force at audio rate (N)
67- malletPos?: Float32Array; // mallet position s along dir at frame times (m)
94+ malletPos?: Float32Array; // mallet position q (3 per frame) at frame times (m)
6895 }
src/ui/instrumentPage.tsmodified+72−7View file
@@ -24,6 +24,7 @@ export interface Built {
2424 sys: BeadSystem;
2525 exc: Excitation;
2626 pickup: Float64Array;
27+ pickupDelay?: Float64Array; // per bead, in samples (a microphone)
2728 scene: SceneSpec;
2829 }
2930
@@ -40,6 +41,7 @@ export interface Extras<P> {
4041 export interface InstrumentConfig<P extends SimParams> {
4142 title: string;
4243 wavName: string;
44+ audioQuantity?: { label: string; unit: string }; // default: acceleration, m/s²
4345 presets: Preset<P>[];
4446 sections: Section<P>[];
4547 // instrument-specific summary of the setup; built = null when the parameters are invalid
@@ -77,7 +79,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
7779 <p class="note" data-r="audio-note"></p>
7880 </section>
7981 <section class="card" data-r="out" hidden>
80- <div class="card-head"><h2>Pickup acceleration</h2>
82+ <div class="card-head"><h2>${cap((cfg.audioQuantity ?? { label: 'pickup acceleration' }).label)}</h2>
8183 <div class="seg" data-r="wave-range"><button data-v="0" class="on">Full</button><button data-v="0.05">First 50 ms</button><button data-v="0.005">First 5 ms</button></div></div>
8284 <div class="canvas-wrap"><canvas data-r="wave"></canvas></div>
8385 </section>
@@ -91,6 +93,11 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
9193 <div class="canvas-wrap"><canvas data-r="spectrum"></canvas></div>
9294 </section>
9395 ${cfg.extraCards ?? ''}
96+ <section class="card" data-r="out-contact" hidden>
97+ <h2>Contact force</h2>
98+ <div class="canvas-wrap short"><canvas data-r="contact"></canvas></div>
99+ <p class="note" data-r="contact-note"></p>
100+ </section>
94101 <section class="card" data-r="out" hidden>
95102 <h2>Energy</h2>
96103 <div class="canvas-wrap short"><canvas data-r="energy"></canvas></div>
@@ -111,6 +118,8 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
111118 const spec = new Spectrogram(q('spec'));
112119 const specPlot = new LinePlot(q('spectrum'));
113120 const energyPlot = new LinePlot(q('energy'));
121+ const contactPlot = new LinePlot(q('contact'));
122+ const quantity = cfg.audioQuantity ?? { label: 'acceleration', unit: 'm/s²' };
114123 const timeline = new Timeline();
115124 const scene = new BeadScene(q('scene'), timeline);
116125 const extras = cfg.mountExtras?.(q, timeline);
@@ -136,11 +145,13 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
136145 const spp = Math.max(1, Math.ceil(1 / (p.fs * ((p.cfl * 2) / wmax))));
137146 rows.push(['Highest frequency bound', `${fmtNum(wmax / (2 * Math.PI))} Hz`]);
138147 rows.push(['Time step', `${fmtNum(1 / (p.fs * spp))} s (${spp} per sample)`]);
139- const work = spp * p.fs * p.duration * built.sys.si.length;
140- rows.push(['Spring updates', fmtNum(work)]);
148+ // a bending hinge costs about four springs
149+ const elements = built.sys.si.length + 4 * (built.sys.bending?.w0.length ?? 0);
150+ const work = spp * p.fs * p.duration * elements;
151+ rows.push(['Element updates', fmtNum(work)]);
141152 const frameMB = (Math.round(p.duration * p.frameRate) * 3 * built.sys.n * 4) / 1e6;
142153 if (frameMB > 400) warnings.push(`Replay frames would take ${fmtNum(frameMB)} MB; lower the replay frame rate or the duration.`);
143- if (work > 3e9) warnings.push('This run is large (more than 3e9 spring updates) and may take minutes.');
154+ if (work > 3e9) warnings.push(`This run is large (${fmtNum(work)} element updates): expect roughly ${fmtNum(work / 2.2e7 / 60)} min on a desktop.`);
144155 }
145156 q('derived').innerHTML = statRows(rows);
146157 q('warnings').innerHTML = warnings.map((w) => `<li>${w}</li>`).join('');
@@ -168,7 +179,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
168179 for (let i = 0; i < t.length; i++) t[i] = i / r.fs;
169180 wave.set([{ x: t, y: r.audio, color: css('--series-1'), label: 'a' }], {
170181 xLabel: 'time (s)',
171- yLabel: 'acceleration (m/s²)',
182+ yLabel: `${quantity.label} (${quantity.unit})`,
172183 xRange: waveRange ? [0, Math.min(waveRange, t[t.length - 1])] : undefined,
173184 });
174185 }
@@ -206,7 +217,8 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
206217 ['Wall time', `${fmtNum(r.wallSeconds)} s (${rt >= 1 ? fmtNum(rt) + '× faster than' : fmtNum(1 / rt) + '× slower than'} real time)`],
207218 ['Energy balance error', emax > 0 ? `${fmtNum((bal / emax) * 100)}% of peak energy` : '–'],
208219 ...(cfg.summaryRows?.(p, built, r) ?? []),
209- ['Peak pickup acceleration', `${fmtNum(peak)} m/s²`],
220+ ...contactRows(r),
221+ [`Peak ${quantity.label}`, `${fmtNum(peak)} ${quantity.unit}`],
210222 ]);
211223 if (audioUrl) URL.revokeObjectURL(audioUrl);
212224 audioUrl = URL.createObjectURL(encodeWav(r.audio, r.fs));
@@ -219,6 +231,7 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
219231 timeline.setTimes(r.frameTimes);
220232 scene.showReplay({ ...built.scene, frames: r.frames, malletPos: r.malletPos });
221233 extras?.show(p, built, r);
234+ showContact(r);
222235 energyPlot.set([{ x: r.frameTimes, y: r.energy, color: css('--series-1'), label: 'energy (J)' }], {
223236 xLabel: 'time (s)',
224237 yLabel: 'energy (J)',
@@ -226,11 +239,33 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
226239 });
227240 }
228241
242+ function showContact(r: SimResult) {
243+ const card = q('out-contact');
244+ card.hidden = !r.contactForce;
245+ if (!r.contactForce) return;
246+ const cf = r.contactForce;
247+ const c = contacts(cf, r.fs);
248+ const end = c.length ? c[c.length - 1].start + c[c.length - 1].ms / 1000 : 0.01;
249+ const nShow = Math.min(cf.length, Math.ceil((end * 1.3 + 0.002) * r.fs));
250+ const t = new Float64Array(nShow);
251+ for (let i = 0; i < nShow; i++) t[i] = (i / r.fs) * 1000;
252+ contactPlot.set([{ x: t, y: cf.subarray(0, nShow), color: css('--series-2'), label: 'force (N)' }], {
253+ xLabel: 'time (ms)',
254+ yLabel: 'contact force (N)',
255+ yRange: [0, Math.max(...cf.subarray(0, nShow)) * 1.08 || 1],
256+ });
257+ q('contact-note').textContent = c.length
258+ ? `${c.length} contact${c.length > 1 ? 's' : ''}: ` +
259+ c.map((k) => `${fmtNum(k.ms)} ms from ${fmtNum(k.start * 1000)} ms, peak ${fmtNum(k.peak)} N`).join('; ') +
260+ '. During contact the instrument takes energy from the mallet and gives part of it back as the mallet rebounds, so the energy plot can spike first; the mallet then leaves at constant speed (no gravity).'
261+ : 'The mallet never touched the instrument.';
262+ }
263+
229264 runBtn.addEventListener('click', async () => {
230265 const p = controls.get();
231266 const { built } = cfg.setup(p);
232267 if (!built) return;
233- const opt = { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: built.pickup };
268+ const opt = { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: built.pickup, pickupDelay: built.pickupDelay };
234269 runBtn.disabled = true;
235270 errorBox.hidden = true;
236271 cancelBtn.hidden = false;
@@ -253,3 +288,33 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
253288 });
254289 cancelBtn.addEventListener('click', () => running?.cancel());
255290 }
291+
292+// separate contact intervals in the contact-force record
293+function contacts(cf: Float32Array, fs: number): { start: number; ms: number; peak: number }[] {
294+ const out: { start: number; ms: number; peak: number }[] = [];
295+ let i = 0;
296+ while (i < cf.length) {
297+ if (cf[i] > 0) {
298+ const a = i;
299+ let peak = 0;
300+ while (i < cf.length && cf[i] > 0) peak = Math.max(peak, cf[i++]);
301+ out.push({ start: a / fs, ms: ((i - a) / fs) * 1000, peak });
302+ } else i++;
303+ }
304+ return out;
305+}
306+
307+function contactRows(r: SimResult): [string, string][] {
308+ if (!r.contactForce) return [];
309+ const c = contacts(r.contactForce, r.fs);
310+ return [
311+ [
312+ 'Contact',
313+ c.length
314+ ? `${fmtNum(c[0].ms)} ms, peak ${fmtNum(c[0].peak)} N${c.length > 1 ? ` (+${c.length - 1} more contact${c.length > 2 ? 's' : ''})` : ''}`
315+ : 'none',
316+ ],
317+ ];
318+}
319+
320+const cap = (s: string) => s[0].toUpperCase() + s.slice(1);
src/ui/membraneview.tsmodified+3−2View file
@@ -11,7 +11,8 @@ export interface MembraneGeometry {
1111 radius: number; // rim radius
1212 spacing: number; // lattice spacing
1313 strike: [number, number] | null; // excitation point
14- pickup: number | null; // pickup bead, null = sum
14+ pickup: number | null; // pickup bead, null = not a single bead
15+ pickupNote?: string; // shown when pickup is null (default: sum over all beads)
1516 sectionAngle: number; // degrees, direction of the cross-section diameter
1617 }
1718
@@ -144,6 +145,6 @@ export class MembraneView {
144145 });
145146 this.note.textContent =
146147 `Colour: vertical (y) displacement, blue down and red up, scaled to ±${fmtNum(this.maxY * 1000)} mm (the largest over the run). ` +
147- `${g.strike ? 'Orange ×: excitation point. ' : ''}${g.pickup !== null ? 'Green ring: pickup bead.' : 'Pickup: sum over all beads.'}`;
148+ `${g.strike ? 'Orange ×: excitation point. ' : ''}${g.pickup !== null ? 'Green ring: pickup bead.' : (g.pickupNote ?? 'Pickup: sum over all beads.')}`;
148149 }
149150 }
src/ui/scene3d.tsmodified+36−12View file
@@ -11,15 +11,17 @@ import { timelineBar, type Timeline } from './timeline';
1111 export interface SceneSpec {
1212 sys: BeadSystem;
1313 exc: Excitation;
14- excKind: 'pluck' | 'strike' | 'hammer' | 'mallet' | 'force' | 'press';
14+ excKind: 'pluck' | 'strike' | 'hammer' | 'mallet' | 'stick' | 'force' | 'press';
1515 pickup: Float64Array; // 3n weights
1616 pickupLabel: string;
1717 hoop?: { radius: number }; // round rim in the x-z plane, centred on the origin (instead of walls)
18+ felt?: { radius: number; gap: number; y: number }; // centre support: felt washers above and below (0, y, 0)
19+ mic?: { dir: [number, number, number]; label: string; from: [number, number, number] }; // direction to a far-field microphone
1820 }
1921
2022 export interface ReplaySpec extends SceneSpec {
2123 frames: Float32Array; // nFrames x 3n, matches the timeline
22- malletPos?: Float32Array; // mallet position along its direction, per frame
24+ malletPos?: Float32Array; // mallet position (3 per frame)
2325 }
2426
2527 type Mode = 'setup' | 'replay';
@@ -58,7 +60,7 @@ export class BeadScene {
5860 private slideRings: { mesh: THREE.Mesh; bead: number }[] = [];
5961 private hammerArrows: THREE.Group | null = null;
6062 private pickupArrow: { arrow: THREE.ArrowHelper; bead: number } | null = null;
61- private mallet: { group: THREE.Group; contact: THREE.Vector3; dir: THREE.Vector3; restS: number } | null = null;
63+ private mallet: { group: THREE.Group; contact: THREE.Vector3; dir: THREE.Vector3 } | null = null;
6264 private spec: SceneSpec | null = null;
6365 private frames: Float32Array | null = null;
6466
@@ -289,7 +291,26 @@ export class BeadScene {
289291 hoop.position.set(-this.center.x * this.scale, -this.center.y * this.scale, -this.center.z * this.scale);
290292 this.model.add(hoop);
291293 }
292- if (!spec.hoop) sys.anchors.forEach((a) => {
294+ if (spec.felt) {
295+ const o = new THREE.Vector3(-this.center.x, spec.felt.y - this.center.y, -this.center.z).multiplyScalar(this.scale);
296+ const fr = spec.felt.radius * this.scale;
297+ const gap = Math.max(spec.felt.gap * this.scale, 0.006);
298+ const feltMat = new THREE.MeshStandardMaterial({ color: cInk, roughness: 0.95, transparent: true, opacity: 0.55 });
299+ for (const side of [1, -1]) {
300+ const washer = new THREE.Mesh(new THREE.CylinderGeometry(fr, fr, 0.012, 32), feltMat);
301+ washer.position.copy(o).add(new THREE.Vector3(0, side * (gap + 0.006), 0));
302+ this.model.add(washer);
303+ }
304+ const rodLen = o.y + 0.22;
305+ const rod = new THREE.Mesh(new THREE.CylinderGeometry(0.006, 0.006, rodLen + 0.06, 10), new THREE.MeshStandardMaterial({ color: cMuted }));
306+ rod.position.copy(o).add(new THREE.Vector3(0, 0.03 - rodLen / 2, 0));
307+ this.model.add(rod);
308+ }
309+ if (spec.mic) {
310+ const o = new THREE.Vector3(...spec.mic.from).sub(this.center).multiplyScalar(this.scale);
311+ this.model.add(new THREE.ArrowHelper(new THREE.Vector3(...spec.mic.dir).normalize(), o, 0.42, new THREE.Color(css('--series-3')), 0.05, 0.025));
312+ }
313+ if (!spec.hoop && !spec.felt) sys.anchors.forEach((a) => {
293314 const at = this.nodePos(sys, sys.n + sys.anchors.indexOf(a), new THREE.Vector3());
294315 const side = at.x < 0 ? -1 : 1;
295316 if (a.slideBead < 0) {
@@ -325,6 +346,7 @@ export class BeadScene {
325346 legend.push(`<span>${sw(cBead)}bead</span>`);
326347 if (pinned.length) legend.push(`<span>${sw(cInk, 'square')}pinned bead</span>`);
327348 if (spec.hoop) legend.push(`<span>${sw(cMuted, 'ring')}clamped rim</span>`);
349+ else if (spec.felt) legend.push(`<span>${sw(cInk, 'square')}felt washers (centre support)</span>`);
328350 else if (sys.anchors.some((a) => a.slideBead < 0)) legend.push(`<span>${sw(cMuted, 'square')}clamped end (wall)</span>`);
329351 if (sys.anchors.some((a) => a.slideBead >= 0)) legend.push(`<span>${sw(cInk, 'ring')}sliding end (ring on a frictionless plane)</span>`);
330352 const size0 = modelSize(sys);
@@ -380,7 +402,7 @@ export class BeadScene {
380402 }
381403 if (exc.mallet) {
382404 const ml = exc.mallet;
383- const dir = new THREE.Vector3(...ml.dir).normalize();
405+ const dir = new THREE.Vector3(...ml.vel0).normalize(); // the stick travels along its own axis
384406 const contact = new THREE.Vector3();
385407 for (let i = 0; i < n; i++) if (ml.weights[i] > 0) contact.addScaledVector(this.nodePos(sys, i, new THREE.Vector3()), ml.weights[i]);
386408 const head = Math.max(1.6 * r, 0.025);
@@ -395,11 +417,11 @@ export class BeadScene {
395417 stick.quaternion.setFromUnitVectors(new THREE.Vector3(0, 1, 0), dir);
396418 group.add(felt, stick);
397419 this.model.add(group);
398- // in setup the mallet hovers above the head
399- const restS = -0.6 * GLYPH;
400- this.mallet = { group, contact, dir, restS };
420+ this.mallet = { group, contact, dir };
421+ const speed = Math.hypot(...ml.vel0);
422+ const what = spec.excKind === 'stick' ? 'stick' : 'mallet';
401423 legend.push(
402- `<span>${sw(cExc)}mallet: ${fmtNum(ml.mass * 1000)} g at ${fmtNum(ml.v0)} m/s, felt F = K δ^${fmtNum(ml.p)}${this.mode === 'replay' ? ' (position exaggerated with the head)' : ''}</span>`,
424+ `<span>${sw(cExc)}${what}: ${fmtNum(ml.mass * 1000)} g at ${fmtNum(speed)} m/s, contact F = K δ^${fmtNum(ml.p)}${this.mode === 'replay' ? ' (motion exaggerated like the beads)' : ''}</span>`,
403425 );
404426 }
405427 if ((spec.excKind === 'hammer' || spec.excKind === 'force') && exc.force) {
@@ -424,7 +446,7 @@ export class BeadScene {
424446 this.model.add(arrow);
425447 this.pickupArrow = { arrow, bead: i };
426448 }
427- legend.push(`<span>${sw(cPick)}pickup: ${spec.pickupLabel}</span>`);
449+ legend.push(spec.mic ? `<span>${sw(cPick, 'arrow')}${spec.mic.label}</span>` : `<span>${sw(cPick)}pickup: ${spec.pickupLabel}</span>`);
428450 this.legend.innerHTML = legend.join('');
429451 this.update();
430452 }
@@ -500,8 +522,10 @@ export class BeadScene {
500522 if (this.pickupArrow) this.pickupArrow.arrow.position.copy(P[this.pickupArrow.bead]);
501523 if (this.mallet) {
502524 const mp = this.mode === 'replay' ? this.replaySpec?.malletPos : undefined;
503- const sNorm = mp ? mp[this.tl.frame] * E : this.mallet.restS;
504- this.mallet.group.position.copy(this.mallet.contact).addScaledVector(this.mallet.dir, sNorm);
525+ const g = this.mallet.group.position.copy(this.mallet.contact);
526+ // in setup it hovers back along its path; in replay it follows the run
527+ if (mp) g.add(new THREE.Vector3(mp[3 * this.tl.frame], mp[3 * this.tl.frame + 1], mp[3 * this.tl.frame + 2]).multiplyScalar(E));
528+ else g.addScaledVector(this.mallet.dir, -0.6 * GLYPH);
505529 }
506530 if (this.hammerArrows && spec.exc.force) {
507531 const f = spec.exc.force;
test/cymbal.test.tsadded+77−0View file
@@ -0,0 +1,77 @@
1+import { describe, expect, it } from 'vitest';
2+import { buildCymbal, buildCymbalExcitation, buildMic, deriveCymbal, type CymbalParams } from '../src/instruments/bead-cymbal/model';
3+import { CYMBAL_PRESETS } from '../src/instruments/bead-cymbal/presets';
4+import { bendingEnergy, simulate } from '../src/sim/solver';
5+import { spectrum } from '../src/ui/fft';
6+
7+const base = (over: Partial<CymbalParams> = {}): CymbalParams => ({ ...CYMBAL_PRESETS[0].params, ...over });
8+
9+describe('cymbal lattice', () => {
10+ it('bending energy matches the plate (nu = 1/3) for bending, doming and twisting', () => {
11+ const p = base({ holeRadius: 0, rings: 12 });
12+ const d = deriveCymbal(p);
13+ const sys = buildCymbal(p, d);
14+ const k = 0.5; // 1/m, small: k * a << 1
15+ const R = p.diameter / 2;
16+ const nu = 1 / 3;
17+ const inner = (t: number) => d.mesh.tris[t].every((i) => Math.hypot(sys.rest[3 * i], sys.rest[3 * i + 2]) < 0.6 * R);
18+ let area = 0;
19+ d.mesh.tris.forEach((tr, t) => {
20+ if (!inner(t)) return;
21+ const [u, v, w] = tr.map((i) => [sys.rest[3 * i], sys.rest[3 * i + 2]]);
22+ area += 0.5 * Math.abs((v[0] - u[0]) * (w[1] - u[1]) - (w[0] - u[0]) * (v[1] - u[1]));
23+ });
24+ const shapes: [string, (x: number, z: number) => number, number][] = [
25+ ['cylinder', (x) => 0.5 * k * x * x, 0.5 * k * k],
26+ ['sphere', (x, z) => 0.5 * k * (x * x + z * z), 0.5 * (4 - 2 * (1 - nu)) * k * k],
27+ ['saddle', (x, z) => 0.5 * k * (x * x - z * z), (1 - nu) * k * k],
28+ ['twist', (x, z) => k * x * z, (1 - nu) * k * k],
29+ ];
30+ for (const [name, f, density] of shapes) {
31+ const X = new Float64Array(3 * (sys.n + sys.anchors.length));
32+ X.set(sys.rest);
33+ for (let i = 0; i < sys.n; i++) X[3 * i + 1] = f(sys.rest[3 * i], sys.rest[3 * i + 2]);
34+ const ratio = bendingEnergy(sys.bending!, X, inner) / (d.D * density * area);
35+ console.log(name, 'bending energy / plate', ratio.toFixed(4));
36+ expect(ratio).toBeCloseTo(1, 2);
37+ }
38+ });
39+
40+ it('a free flat plate rings near the (2,0) mode of a free circular plate', () => {
41+ // no felts, no hole, no damping; initial shape r^2 cos(2 phi), small
42+ const p = base({ holeRadius: 0, rings: 9, feltRadius: -1, t60Low: 0, t60High: 0, duration: 0.3, frameRate: 50 });
43+ const d = deriveCymbal(p);
44+ const sys = buildCymbal(p, d);
45+ const exc = buildCymbalExcitation({ ...p, excitation: 'strike', strikeVel: 0 }, d);
46+ for (let i = 0; i < sys.n; i++) {
47+ const x = sys.rest[3 * i], z = sys.rest[3 * i + 2];
48+ exc.u0[3 * i + 1] = 1e-7 * (x * x - z * z) / 0.01;
49+ }
50+ const pick = new Float64Array(3 * sys.n);
51+ let edge = 0;
52+ for (let i = 0; i < sys.n; i++) if (sys.rest[3 * i] > sys.rest[3 * edge]) edge = i;
53+ pick[3 * edge + 1] = 1;
54+ const r = simulate(sys, exc, { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: pick });
55+ const s = spectrum(r.audio, r.fs);
56+ let best = 1;
57+ for (let k = 1; k < s.freq.length && s.freq[k] < 300; k++) if (s.db[k] > s.db[best]) best = k;
58+ const R = p.diameter / 2;
59+ const ideal = (5.253 / (2 * Math.PI * R * R)) * d.bendingCoef; // Leissa, nu ~ 1/3
60+ console.log('free (2,0) measured', s.freq[best].toFixed(1), 'Hz, ideal', ideal.toFixed(1), 'Hz');
61+ expect(Math.abs(s.freq[best] / ideal - 1)).toBeLessThan(0.05);
62+ }, 60000);
63+
64+ it('energy balances with a glancing stick hit, friction, felts and damping', () => {
65+ const p = base({ rings: 7, hitAngle: 60, stickSpeed: 6, bowHeight: 0.01, bellHeight: 0.01, duration: 0.03, frameRate: 2000 });
66+ const d = deriveCymbal(p);
67+ const mic = buildMic(p, d);
68+ const r = simulate(buildCymbal(p, d), buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });
69+ const emax = Math.max(...r.energy);
70+ const bal = Math.max(...r.balance.map(Math.abs));
71+ const cf = r.contactForce!;
72+ console.log('contact samples', cf.reduce((c, v) => c + (v > 0 ? 1 : 0), 0), 'peak N', Math.max(...cf).toFixed(0), 'balance', (bal / emax).toExponential(2));
73+ expect(bal / emax).toBeLessThan(1e-2);
74+ expect(cf[cf.length - 1]).toBe(0);
75+ expect(r.audio.some((v) => v !== 0)).toBe(true);
76+ });
77+});
test/drum.test.tsmodified+1−1View file
@@ -42,7 +42,7 @@ describe('drum', () => {
4242 const nz = cf.reduce((c, v) => c + (v > 0 ? 1 : 0), 0);
4343 console.log('contact ms', ((nz / r.fs) * 1000).toFixed(2), 'peak N', Math.max(...cf).toFixed(1), 'balance', Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy));
4444 expect(cf[cf.length - 1]).toBe(0);
45- expect(r.malletPos![r.malletPos!.length - 1]).toBeLessThan(0);
45+ expect(r.malletPos![r.malletPos!.length - 2]).toBeGreaterThan(0); // y of the last frame: back above the head
4646 expect(Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy)).toBeLessThan(1e-2);
4747 });
4848 });
test/hinge.test.tsadded+38−0View file
@@ -0,0 +1,38 @@
1+import { describe, expect, it } from 'vitest';
2+import { hingeAngleGrad } from '../src/sim/hinge';
3+
4+// random bent configuration around an edge along x
5+function config(seed: number): Float64Array {
6+ let s = seed;
7+ const r = () => ((s = (s * 16807) % 2147483647) / 2147483647 - 0.5) * 0.3;
8+ // w0, w1, e0, e1
9+ return Float64Array.from([0.5 + r(), 0.8 + r(), r(), 0.4 + r(), -0.7 + r(), 0.3 + r(), r(), r(), r(), 1 + r(), r(), r()]);
10+}
11+
12+describe('hinge', () => {
13+ it('angle is zero when flat and matches a fold of known angle', () => {
14+ const flat = Float64Array.from([0.5, 0, 1, 0.5, 0, -1, 0, 0, 0, 1, 0, 0]);
15+ expect(hingeAngleGrad(flat, 0, 1, 2, 3)).toBeCloseTo(0, 12);
16+ const a = 0.3; // fold w1 down by angle a about the x axis
17+ const fold = Float64Array.from([0.5, 0, 1, 0.5, -Math.sin(a), -Math.cos(a), 0, 0, 0, 1, 0, 0]);
18+ expect(Math.abs(hingeAngleGrad(fold, 0, 1, 2, 3))).toBeCloseTo(a, 12);
19+ });
20+ it('gradient matches finite differences', () => {
21+ for (let seed = 1; seed <= 5; seed++) {
22+ const X = config(seed * 7919);
23+ const out = new Float64Array(13);
24+ hingeAngleGrad(X, 0, 1, 2, 3, out);
25+ const order = [0, 1, 2, 3]; // node k's gradient is at out[1 + 3k]
26+ const h = 1e-6;
27+ for (const node of order)
28+ for (let c = 0; c < 3; c++) {
29+ const Xp = X.slice();
30+ const Xm = X.slice();
31+ Xp[3 * node + c] += h;
32+ Xm[3 * node + c] -= h;
33+ const fd = (hingeAngleGrad(Xp, 0, 1, 2, 3) - hingeAngleGrad(Xm, 0, 1, 2, 3)) / (2 * h);
34+ expect(out[1 + 3 * node + c]).toBeCloseTo(fd, 6);
35+ }
36+ }
37+ });
38+});