c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 1import type { BeadSystem, Bending, Excitation, Mallet, SimOptions, SimResult } from './types';
2import { hingeAngleGrad } from './hinge';
4// Normal speed of a mallet into the surface, and its linearized contact
5// stiffness at the largest compression a rigid surface would see.
6function 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}
13function 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}
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 19
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 20// Linearized stiffness at rest, y = K x (x, y: 3n, beads only; anchors held).
21// A spring under tension T has stiffness k along its axis and T/l across it;
22// bending contributes G^T M G per triangle (G = hinge-angle gradients).
23function restStiffness(sys: BeadSystem): (x: Float64Array, y: Float64Array) => void {
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 24 const { n, si, sj, sk, sL0 } = sys;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 25 const X = nodePositions(sys);
27 const sd = new Float64Array(3 * ns); // unit axis
28 const ka = new Float64Array(ns); // axial stiffness minus transverse
29 const kt = new Float64Array(ns); // transverse T / l
30 for (let s = 0; s < ns; s++) {
31 const a = 3 * si[s], b = 3 * sj[s];
32 const dx = X[b] - X[a], dy = X[b + 1] - X[a + 1], dz = X[b + 2] - X[a + 2];
33 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
34 sd[3 * s] = dx / l;
35 sd[3 * s + 1] = dy / l;
36 sd[3 * s + 2] = dz / l;
37 kt[s] = (sk[s] * (l - sL0[s])) / l;
38 ka[s] = sk[s] - kt[s];
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 39 }
41 const nh = b ? b.w0.length : 0;
42 const G = new Float64Array(12 * nh);
43 const nodes = new Int32Array(4 * nh);
44 if (b) {
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 45 const g = new Float64Array(13);
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 46 for (let e = 0; e < nh; e++) {
47 hingeAngleGrad(X, b.w0[e], b.w1[e], b.e0[e], b.e1[e], g);
48 for (let q = 0; q < 12; q++) G[12 * e + q] = b.sign[e] * g[1 + q];
49 nodes.set([b.w0[e], b.w1[e], b.e0[e], b.e1[e]], 4 * e);
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 50 }
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 51 }
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 52 const phi = new Float64Array(nh);
53 const mom = new Float64Array(nh);
54 const get = (x: Float64Array, node: number, c: number) => (node < n ? x[3 * node + c] : 0);
55 return (x, y) => {
56 y.fill(0);
57 for (let s = 0; s < ns; s++) {
58 const a = si[s], c = sj[s];
59 const rx = get(x, c, 0) - get(x, a, 0), ry = get(x, c, 1) - get(x, a, 1), rz = get(x, c, 2) - get(x, a, 2);
60 const ex = sd[3 * s], ey = sd[3 * s + 1], ez = sd[3 * s + 2];
61 const pa = ka[s] * (ex * rx + ey * ry + ez * rz);
62 const fx = pa * ex + kt[s] * rx, fy = pa * ey + kt[s] * ry, fz = pa * ez + kt[s] * rz;
63 if (a < n) {
64 y[3 * a] -= fx;
65 y[3 * a + 1] -= fy;
66 y[3 * a + 2] -= fz;
67 }
68 if (c < n) {
69 y[3 * c] += fx;
70 y[3 * c + 1] += fy;
71 y[3 * c + 2] += fz;
72 }
73 }
74 if (!b) return;
75 for (let e = 0; e < nh; e++) {
76 let v = 0;
77 for (let k = 0; k < 4; k++) {
78 const nd = nodes[4 * e + k];
79 if (nd < n) v += G[12 * e + 3 * k] * x[3 * nd] + G[12 * e + 3 * k + 1] * x[3 * nd + 1] + G[12 * e + 3 * k + 2] * x[3 * nd + 2];
80 }
81 phi[e] = v;
82 }
83 mom.fill(0);
84 for (let t = 0; t < b.tri.length / 3; t++)
85 for (let a = 0; a < 3; a++) {
86 const ea = b.tri[3 * t + a];
87 if (ea < 0) continue;
88 for (let c = 0; c < 3; c++) {
89 const ec = b.tri[3 * t + c];
90 if (ec >= 0) mom[ea] += b.M[9 * t + 3 * a + c] * phi[ec];
91 }
92 }
93 for (let e = 0; e < nh; e++)
94 for (let k = 0; k < 4; k++) {
95 const nd = nodes[4 * e + k];
96 if (nd < n) for (let c = 0; c < 3; c++) y[3 * nd + c] += mom[e] * G[12 * e + 3 * k + c];
97 }
98 };
99}
101// Highest angular frequency of the linearized system at rest: power iteration
102// on M^-1/2 K M^-1/2 (seeded, so the same system always gives the same value).
103// The Rayleigh quotient approaches the top eigenvalue from below; the contact
104// stiffness of a mallet is added on top (Weyl: lambda(A + B) <= lambda(A) + lambda(B)).
105export function omegaMaxEstimate(sys: BeadSystem, mallet?: Mallet, iterations = 400): number {
106 const { n, pinned } = sys;
107 const K = restStiffness(sys);
108 const rs = Float64Array.from(sys.mass, (m, i) => (pinned[i] ? 0 : 1 / Math.sqrt(m)));
109 let seed = 12345;
110 const rand = () => ((seed = (seed * 1103515245 + 12345) % 2147483648) / 2147483648) - 0.5;
111 const v = new Float64Array(3 * n);
112 const z = new Float64Array(3 * n);
113 const y = new Float64Array(3 * n);
114 for (let k = 0; k < 3 * n; k++) v[k] = rs[Math.floor(k / 3)] > 0 ? rand() : 0;
115 let lam = 0;
116 for (let it = 0; it < iterations; it++) {
117 let nv = 0;
118 for (let k = 0; k < 3 * n; k++) nv += v[k] * v[k];
119 nv = Math.sqrt(nv);
120 if (!(nv > 0)) break;
121 for (let k = 0; k < 3 * n; k++) z[k] = (v[k] / nv) * rs[Math.floor(k / 3)];
122 K(z, y);
123 lam = 0;
124 for (let k = 0; k < 3 * n; k++) {
125 y[k] *= rs[Math.floor(k / 3)];
126 lam += (y[k] * v[k]) / nv;
127 }
128 v.set(y);
129 }
130 let w2 = Math.max(lam, 0);
131 if (mallet) {
132 let inv = 1 / mallet.mass;
133 for (let i = 0; i < n; i++) if (!pinned[i] && mallet.weights[i] > 0) inv += (mallet.weights[i] * mallet.weights[i]) / sys.mass[i];
134 w2 += contactStiffness(mallet) * inv;
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 135 }
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 136 return Math.sqrt(w2);
137}
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 139export interface TimeStep {
140 omegaMax: number;
141 dt: number;
142 stepsPerSample: number;
143}
145// dt = cfl * 2 / omega_max, rounded down to a whole number of steps per audio sample
146export function timeStep(sys: BeadSystem, mallet: Mallet | undefined, fs: number, cfl: number): TimeStep {
147 const omegaMax = omegaMaxEstimate(sys, mallet);
148 const stepsPerSample = Math.max(1, Math.ceil(omegaMax / (fs * cfl * 2)));
149 return { omegaMax, dt: 1 / (fs * stepsPerSample), stepsPerSample };
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 150}
152// bending energy at positions X, optionally only over triangles with keep(t)
153export function bendingEnergy(b: Bending, X: Float64Array, keep?: (t: number) => boolean): number {
154 const g = new Float64Array(13);
155 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]));
156 let E = 0;
157 for (let t = 0; t < b.tri.length / 3; t++) {
158 if (keep && !keep(t)) continue;
159 for (let a = 0; a < 3; a++)
160 for (let c = 0; c < 3; c++) {
161 const ea = b.tri[3 * t + a];
162 const ec = b.tri[3 * t + c];
163 if (ea >= 0 && ec >= 0) E += 0.5 * phi[ea] * b.M[9 * t + 3 * a + c] * phi[ec];
164 }
165 }
166 return E;
167}
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 169export type Progress = (fraction: number) => void;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 171// Velocity Verlet (kick-drift-kick), explicit. Spring and hinge forces are
172// exact (geometric nonlinearity kept); damping and friction use the half-step
173// velocity in the kick, and the full-step velocity whenever an audio sample is
174// recorded.
175export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, progress?: Progress): SimResult {
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 176 const t0 = performance.now();
177 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
178 const nn = n + anchors.length;
179 const ns = si.length;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 180 const bend = sys.bending;
181 const nh = bend ? bend.w0.length : 0;
182 const nt = bend ? bend.tri.length / 3 : 0;
183 const hG = new Float64Array(12 * nh); // gradients of phi
184 const hPhi = new Float64Array(nh);
185 const hRate = new Float64Array(nh);
186 const hMom = new Float64Array(nh);
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 187 const hMomEl = new Float64Array(nh); // elastic part of the moments
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 188
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 189 const mallet = exc.mallet;
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 190 const { omegaMax, dt, stepsPerSample } = timeStep(sys, mallet, opt.fs, opt.cfl);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 191 const nSamples = Math.round(opt.duration * opt.fs);
192 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
193 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
195 const X = new Float64Array(3 * nn);
196 const V = new Float64Array(3 * nn);
197 const F = new Float64Array(3 * n); // force on beads
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 198 const Fel = new Float64Array(3 * n); // elastic (position-dependent) part
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 199 const A = new Float64Array(3 * n);
200 for (let i = 0; i < n; i++) {
201 for (let c = 0; c < 3; c++) {
202 const k = 3 * i + c;
203 X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
204 V[k] = pinned[i] ? 0 : exc.v0[k];
205 }
206 }
207 const slide: number[] = [];
208 for (let a = 0; a < anchors.length; a++) {
209 for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
210 if (anchors[a].slideBead >= 0) slide.push(a);
211 }
212 const syncAnchors = () => {
213 for (const a of slide) {
214 const b = anchors[a].slideBead;
215 const k = 3 * (n + a);
216 X[k + 1] = X[3 * b + 1];
217 X[k + 2] = X[3 * b + 2];
218 V[k + 1] = V[3 * b + 1];
219 V[k + 2] = V[3 * b + 2];
220 }
221 };
222 syncAnchors();
224 const force = exc.force;
225 const fShape = force?.shape;
226 const pulse = (t: number) =>
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 227 force && t >= 0 && t <= force.tau ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau)) : 0;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 228
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 229 // mallet state: 3D position q, velocity qv, acceleration qa
230 const q = new Float64Array(3);
231 const qv = new Float64Array(mallet ? mallet.vel0 : [0, 0, 0]);
232 const qa = new Float64Array(3);
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 233 let fContact = 0;
234 let eContact = 0;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 235 let pFriction = 0;
236 const patch: number[] = [];
237 let mEff = 1;
238 if (mallet) {
239 let inv = 0;
240 for (let i = 0; i < n; i++)
241 if (mallet.weights[i] > 0) {
242 patch.push(i);
243 inv += (mallet.weights[i] * mallet.weights[i]) / mass[i];
244 }
245 mEff = 1 / (1 / mallet.mass + inv);
246 }
247 const [nx, ny, nz] = mallet?.normal ?? [0, 0, 0];
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 248 const contact = () => {
249 fContact = 0;
250 eContact = 0;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 251 pFriction = 0;
252 qa.fill(0);
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 253 if (!mallet) return;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 254 const w = mallet.weights;
255 let ux = 0, uy = 0, uz = 0, vx = 0, vy = 0, vz = 0;
256 for (const i of patch) {
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 257 const k = 3 * i;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 258 ux += w[i] * (X[k] - sys.rest[k]);
259 uy += w[i] * (X[k + 1] - sys.rest[k + 1]);
260 uz += w[i] * (X[k + 2] - sys.rest[k + 2]);
261 vx += w[i] * V[k];
262 vy += w[i] * V[k + 1];
263 vz += w[i] * V[k + 2];
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 264 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 265 const d = (q[0] - ux) * nx + (q[1] - uy) * ny + (q[2] - uz) * nz;
266 if (!(d > 0)) return;
267 if (mallet.reach !== undefined) {
268 const ox = q[0] - ux - d * nx, oy = q[1] - uy - d * ny, oz = q[2] - uz - d * nz;
269 if (ox * ox + oy * oy + oz * oz > mallet.reach * mallet.reach) return;
270 }
271 const fn = mallet.K * d ** mallet.p;
272 fContact = fn;
273 eContact = (fn * d) / (mallet.p + 1);
274 // force on the mallet: -fn normal + friction against tangential slip
275 let fx = -fn * nx;
276 let fy = -fn * ny;
277 let fz = -fn * nz;
278 if (mallet.mu > 0) {
279 let rx = qv[0] - vx, ry = qv[1] - vy, rz = qv[2] - vz;
280 const rn = rx * nx + ry * ny + rz * nz;
281 rx -= rn * nx;
282 ry -= rn * ny;
283 rz -= rn * nz;
284 const slip = Math.hypot(rx, ry, rz);
285 // regularized Coulomb; the slip scale keeps the explicit step stable
286 const eps = Math.max(1e-3, (4 * mallet.mu * fn * dt) / mEff);
287 const g = (mallet.mu * fn) / Math.sqrt(slip * slip + eps * eps);
288 fx -= g * rx;
289 fy -= g * ry;
290 fz -= g * rz;
291 pFriction = g * slip * slip;
292 }
293 qa[0] = fx / mallet.mass;
294 qa[1] = fy / mallet.mass;
295 qa[2] = fz / mallet.mass;
296 for (const i of patch) {
297 F[3 * i] -= w[i] * fx;
298 F[3 * i + 1] -= w[i] * fy;
299 F[3 * i + 2] -= w[i] * fz;
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 300 }
301 };
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 303 const hOut = new Float64Array(13);
304 let eHinge = 0;
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 305 const anyDashpot = sc.some((c) => c !== 0);
306 // Position-dependent forces: spring and bending elasticity into Fel; hinge
307 // angles, gradients and elastic moments; rigid-body geometry.
308 const elasticForces = () => {
309 Fel.fill(0);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 310 for (let s = 0; s < ns; s++) {
311 const a = 3 * si[s];
312 const b = 3 * sj[s];
313 const dx = X[b] - X[a];
314 const dy = X[b + 1] - X[a + 1];
315 const dz = X[b + 2] - X[a + 2];
316 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
317 const g = (sk[s] * (l - sL0[s])) / l;
318 if (si[s] < n) {
320 Fel[a + 1] += g * dy;
321 Fel[a + 2] += g * dz;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 322 }
323 if (sj[s] < n) {
325 Fel[b + 1] -= g * dy;
326 Fel[b + 2] -= g * dz;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 327 }
328 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 329 eHinge = 0;
330 if (bend) {
331 for (let e = 0; e < nh; e++) {
332 const th = hingeAngleGrad(X, bend.w0[e], bend.w1[e], bend.e0[e], bend.e1[e], hOut);
333 const sg = bend.sign[e];
334 hPhi[e] = sg * (th - bend.theta0[e]);
335 const o = 12 * e;
336 for (let q = 0; q < 12; q++) hG[o + q] = sg * hOut[1 + q];
337 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 339 const tri = bend.tri;
340 const Mt = bend.M;
341 for (let t = 0; t < nt; t++) {
342 for (let a = 0; a < 3; a++) {
343 const ea = tri[3 * t + a];
344 if (ea < 0) continue;
345 let m = 0;
346 for (let c = 0; c < 3; c++) {
347 const ec = tri[3 * t + c];
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 348 if (ec >= 0) m += Mt[9 * t + 3 * a + c] * hPhi[ec];
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 349 }
350 eHinge += 0.5 * hPhi[ea] * m;
352 }
353 }
354 }
355 if (sys.freeRigid) rigidGeometry();
356 };
358 // Fills F with total force on beads and A with accelerations, from the
359 // elastic part (elasticForces must have run at this position) and the
360 // velocity-dependent part at V; sets damping power (W, >= 0) and external
361 // power (W), both at V.
362 let pDamp = 0;
363 let pExt = 0;
364 const dynamicForces = (t: number) => {
365 F.set(Fel);
366 contact();
367 pDamp = pFriction;
368 if (anyDashpot)
369 for (let s = 0; s < ns; s++) {
370 const c = sc[s];
371 if (c === 0) continue;
372 const a = 3 * si[s];
373 const b = 3 * sj[s];
374 const rx = V[b] - V[a];
375 const ry = V[b + 1] - V[a + 1];
376 const rz = V[b + 2] - V[a + 2];
377 pDamp += c * (rx * rx + ry * ry + rz * rz);
378 if (si[s] < n) {
379 F[a] += c * rx;
380 F[a + 1] += c * ry;
381 F[a + 2] += c * rz;
382 }
383 if (sj[s] < n) {
384 F[b] -= c * rx;
385 F[b + 1] -= c * ry;
386 F[b + 2] -= c * rz;
387 }
388 }
389 if (bend) {
390 const s1 = bend.sigma1;
391 hMom.set(hMomEl);
392 if (s1 !== 0) {
393 for (let e = 0; e < nh; e++) {
394 const o = 12 * e;
395 const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
396 hRate[e] =
397 hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
398 hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
399 hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
400 hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
401 }
402 const tri = bend.tri;
403 const Mt = bend.M;
404 for (let t = 0; t < nt; t++) {
405 for (let a = 0; a < 3; a++) {
406 const ea = tri[3 * t + a];
407 if (ea < 0) continue;
408 let r = 0;
409 for (let c = 0; c < 3; c++) {
410 const ec = tri[3 * t + c];
411 if (ec >= 0) r += Mt[9 * t + 3 * a + c] * hRate[ec];
412 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 413 pDamp += 2 * s1 * hRate[ea] * r;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 415 }
416 }
417 }
418 for (let e = 0; e < nh; e++) {
419 const m = hMom[e];
420 if (m === 0) continue;
421 const o = 12 * e;
422 const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
423 F[w0] -= m * hG[o];
424 F[w0 + 1] -= m * hG[o + 1];
425 F[w0 + 2] -= m * hG[o + 2];
426 F[w1] -= m * hG[o + 3];
427 F[w1 + 1] -= m * hG[o + 4];
428 F[w1 + 2] -= m * hG[o + 5];
429 F[e0] -= m * hG[o + 6];
430 F[e0 + 1] -= m * hG[o + 7];
431 F[e0 + 2] -= m * hG[o + 8];
432 F[e1] -= m * hG[o + 9];
433 F[e1 + 1] -= m * hG[o + 10];
434 F[e1 + 2] -= m * hG[o + 11];
435 }
436 }
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 437 const s0 = sys.sigma0;
438 const p = fShape ? pulse(t) : 0;
439 pExt = 0;
440 for (let i = 0; i < n; i++) {
441 const k = 3 * i;
442 if (pinned[i]) {
443 F[k] = F[k + 1] = F[k + 2] = 0;
444 continue;
445 }
446 if (s0 !== 0) {
447 const m2 = 2 * s0 * mass[i];
448 for (let c = 0; c < 3; c++) {
449 F[k + c] -= m2 * V[k + c];
450 pDamp += m2 * V[k + c] * V[k + c];
451 }
452 }
453 if (p !== 0) {
454 for (let c = 0; c < 3; c++) {
455 const fe = p * fShape![k + c];
456 F[k + c] += fe;
457 pExt += fe * V[k + c];
458 }
459 }
460 }
461 for (let i = 0; i < n; i++) {
462 const im = 1 / mass[i];
463 A[3 * i] = F[3 * i] * im;
464 A[3 * i + 1] = F[3 * i + 1] * im;
465 A[3 * i + 2] = F[3 * i + 2] * im;
466 }
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 467 if (sys.freeRigid) removeRigid(A);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 468 };
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 469 const computeForces = (t: number) => {
470 elasticForces();
471 dynamicForces(t);
472 };
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 473
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 474 // centre of mass and inverse inertia tensor at the current positions
475 let totalMass = 0;
476 for (let i = 0; i < n; i++) totalMass += mass[i];
477 const cm = new Float64Array(3);
478 const Iinv = new Float64Array(6);
479 const rigidGeometry = () => {
480 let cx = 0, cy = 0, cz = 0;
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 481 for (let i = 0; i < n; i++) {
482 cx += mass[i] * X[3 * i];
483 cy += mass[i] * X[3 * i + 1];
484 cz += mass[i] * X[3 * i + 2];
485 }
487 cy /= totalMass;
488 cz /= totalMass;
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 489 let ixx = 0, iyy = 0, izz = 0, ixy = 0, ixz = 0, iyz = 0;
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 490 for (let i = 0; i < n; i++) {
491 const m = mass[i];
492 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
493 ixx += m * (ry * ry + rz * rz);
494 iyy += m * (rx * rx + rz * rz);
495 izz += m * (rx * rx + ry * ry);
496 ixy -= m * rx * ry;
497 ixz -= m * rx * rz;
498 iyz -= m * ry * rz;
499 }
500 cm[0] = cx;
501 cm[1] = cy;
502 cm[2] = cz;
503 invert3(ixx, ixy, ixz, iyy, iyz, izz, Iinv);
504 };
506 // subtract the rigid-body part (mass-weighted) of a per-bead field u, about
507 // the geometry from rigidGeometry; returns the kinetic energy of that part
508 // when u is a velocity
509 const removeRigid = (u: Float64Array): number => {
510 const [cx, cy, cz] = cm;
511 let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 512 for (let i = 0; i < n; i++) {
513 const m = mass[i];
514 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
515 const ux = u[3 * i], uy = u[3 * i + 1], uz = u[3 * i + 2];
516 px += m * ux;
517 py += m * uy;
518 pz += m * uz;
519 lx += m * (ry * uz - rz * uy);
520 ly += m * (rz * ux - rx * uz);
521 lz += m * (rx * uy - ry * ux);
522 }
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 523 const [a, b, c, d, e, f] = Iinv;
524 const wx = a * lx + b * ly + c * lz;
525 const wy = b * lx + d * ly + e * lz;
526 const wz = c * lx + e * ly + f * lz;
527 const vx = px / totalMass, vy = py / totalMass, vz = pz / totalMass;
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 528 for (let i = 0; i < n; i++) {
529 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
530 u[3 * i] -= vx + (wy * rz - wz * ry);
531 u[3 * i + 1] -= vy + (wz * rx - wx * rz);
532 u[3 * i + 2] -= vz + (wx * ry - wy * rx);
533 }
534 return 0.5 * (px * vx + py * vy + pz * vz) + 0.5 * (wx * lx + wy * ly + wz * lz);
535 };
537 rigidGeometry();
538 removeRigid(V); // an initial strike carries no rigid motion either
539 }
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 540
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 541 const velocityForces = sys.sigma0 !== 0 || anyDashpot || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 542 let maxStretch = 0;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 543 const malletEnergy = () => (mallet ? 0.5 * mallet.mass * (qv[0] ** 2 + qv[1] ** 2 + qv[2] ** 2) + eContact : 0);
544 // total mechanical energy at the current state; computeForces must have run
545 // at this state (it sets the hinge and contact energies)
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 546 const energyNow = () => {
547 let e = 0;
548 for (let i = 0; i < n; i++) {
549 const k = 3 * i;
550 e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
551 }
552 for (let s = 0; s < ns; s++) {
553 const a = 3 * si[s];
554 const b = 3 * sj[s];
555 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
556 const d = l - sL0[s];
557 e += 0.5 * sk[s] * d * d;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 558 if (sj[s] < n && si[s] < n) maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 559 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 560 return e + eHinge + malletEnergy();
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 561 };
563 const audio = new Float32Array(nSamples);
564 const frames = new Float32Array(nFrames * 3 * n);
565 const frameTimes = new Float64Array(nFrames);
566 const energy = new Float64Array(nFrames);
567 const balance = new Float64Array(nFrames);
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 568 const contactForce = mallet ? new Float32Array(nSamples) : undefined;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 569 const malletPos = mallet ? new Float32Array(3 * nFrames) : undefined;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 570 const w = opt.pickup;
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 571 const delay = opt.pickupDelay;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 572
573 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
574 let Eeq = 0;
575 for (let s = 0; s < ns; s++) {
576 const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
577 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));
578 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
579 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 580 computeForces(0);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 581 const E0 = energyNow();
582 let dissipated = 0;
583 let workIn = 0;
584 let t = 0;
585 let frame = 0;
586 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
588 for (let smp = 0; smp < nSamples; smp++) {
589 // record at t = smp / fs; A holds forces at full-step velocity here
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 590 if (delay) {
591 for (let i = 0; i < n; i++) {
592 const k = 3 * i;
593 const c = w[k] * A[k] + w[k + 1] * A[k + 1] + w[k + 2] * A[k + 2];
594 if (c === 0) continue;
595 const pos = smp + delay[i];
596 const i0 = Math.floor(pos);
597 const fr = pos - i0;
598 if (i0 < nSamples) audio[i0] += (1 - fr) * c;
599 if (i0 + 1 < nSamples) audio[i0 + 1] += fr * c;
600 }
601 } else {
602 let acc = 0;
603 for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
604 audio[smp] = acc;
605 }
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 606 if (contactForce) contactForce[smp] = fContact;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 607 if (smp % samplesPerFrame === 0 && frame < nFrames) {
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 608 if (malletPos) malletPos.set(q, 3 * frame);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 609 const off = frame * 3 * n;
610 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
611 frameTimes[frame] = t;
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 612 const eTotal = energyNow();
613 // the plot shows the instrument's own energy; the balance includes the mallet
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 614 energy[frame] = eTotal - Eeq - malletEnergy();
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 615 balance[frame] = eTotal + dissipated - workIn - E0;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 616 frame++;
617 }
618 if (smp === nSamples - 1) break;
620 for (let st = 0; st < stepsPerSample; st++) {
621 // kick (half), drift
622 for (let k = 0; k < 3 * n; k++) {
623 V[k] += 0.5 * dt * A[k];
624 X[k] += dt * V[k];
625 }
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 626 if (mallet)
627 for (let c = 0; c < 3; c++) {
628 qv[c] += 0.5 * dt * qa[c];
629 q[c] += dt * qv[c];
630 }
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 631 syncAnchors();
632 t += dt;
633 // forces at new position with half-step velocity; powers use it too
634 computeForces(t);
635 dissipated += dt * pDamp;
636 workIn += dt * pExt;
637 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 638 if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 639 // the rigid-free subspace turns with the body; drop what drifts out of it
640 if (sys.freeRigid) dissipated += removeRigid(V);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 641 syncAnchors();
642 }
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 643 // at a sample time: recompute with the full-step velocity (positions unchanged)
644 if (velocityForces) dynamicForces(t);
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 645 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
646 }
648 return {
649 audio,
650 fs: opt.fs,
651 frames,
652 frameTimes,
653 energy,
654 balance,
655 dt,
656 stepsPerSample,
657 omegaMax,
658 wallSeconds: (performance.now() - t0) / 1000,
659 maxStretch,
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 660 contactForce,
661 malletPos,
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 662 };
663}
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 664
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 665// inverse of the symmetric 3x3 matrix [a b c; b d e; c e f], as its six
666// distinct entries (zero if singular)
667function invert3(a: number, b: number, c: number, d: number, e: number, f: number, out: Float64Array) {
591834dCymbal: free plate with rigid-body motion removed (no support); lighter default dampingJeremy Magland 668 const A00 = d * f - e * e, A01 = c * e - b * f, A02 = b * e - c * d;
669 const A11 = a * f - c * c, A12 = b * c - a * e, A22 = a * d - b * b;
670 const det = a * A00 + b * A01 + c * A02;
a78afe0Add WebGPU solver with a toggle; faster CPU solverJeremy Magland 671 if (!(Math.abs(det) > 0)) {
672 out.fill(0);
673 return;
674 }
675 out[0] = A00 / det;
676 out[1] = A01 / det;
677 out[2] = A02 / det;
678 out[3] = A11 / det;
679 out[4] = A12 / det;
680 out[5] = A22 / det;