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}
20// Upper bound on the highest angular frequency of the linearized system
21// (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
22// 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.
25export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
26 const { n, si, sj, sk, sL0 } = sys;
27 const rowSum = new Float64Array(n);
28 const X = nodePositions(sys);
29 for (let s = 0; s < si.length; s++) {
30 const a = si[s];
31 const b = sj[s];
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]);
33 const T = sk[s] * (l - sL0[s]);
34 const kmax = Math.max(sk[s], Math.abs(T) / l);
35 if (a < n) rowSum[a] += 2 * kmax;
36 if (b < n) rowSum[b] += 2 * kmax;
37 }
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 }
52 }
53 const kc = mallet ? contactStiffness(mallet) : 0;
54 let w2 = mallet ? kc / mallet.mass : 0;
55 for (let i = 0; i < n; i++) {
56 if (sys.pinned[i]) continue;
57 const extra = mallet ? kc * mallet.weights[i] : 0;
58 w2 = Math.max(w2, (rowSum[i] + extra) / sys.mass[i]);
59 }
60 return Math.sqrt(w2);
61}
63// per hinge: sum of |M| entries of its rows over the triangles that hold it
64function 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}
75// bending energy at positions X, optionally only over triangles with keep(t)
76export 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}
92export type Progress = (fraction: number) => void;
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.
98export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, progress?: Progress): SimResult {
99 const t0 = performance.now();
100 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
101 const nn = n + anchors.length;
102 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);
111 const mallet = exc.mallet;
112 const omegaMax = omegaMaxBound(sys, mallet);
113 const dtMax = (opt.cfl * 2) / omegaMax;
114 const stepsPerSample = Math.max(1, Math.ceil(1 / (opt.fs * dtMax)));
115 const dt = 1 / (opt.fs * stepsPerSample);
116 const nSamples = Math.round(opt.duration * opt.fs);
117 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
118 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
120 const X = new Float64Array(3 * nn);
121 const V = new Float64Array(3 * nn);
122 const F = new Float64Array(3 * n); // force on beads
123 const A = new Float64Array(3 * n);
124 for (let i = 0; i < n; i++) {
125 for (let c = 0; c < 3; c++) {
126 const k = 3 * i + c;
127 X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
128 V[k] = pinned[i] ? 0 : exc.v0[k];
129 }
130 }
131 const slide: number[] = [];
132 for (let a = 0; a < anchors.length; a++) {
133 for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
134 if (anchors[a].slideBead >= 0) slide.push(a);
135 }
136 const syncAnchors = () => {
137 for (const a of slide) {
138 const b = anchors[a].slideBead;
139 const k = 3 * (n + a);
140 X[k + 1] = X[3 * b + 1];
141 X[k + 2] = X[3 * b + 2];
142 V[k + 1] = V[3 * b + 1];
143 V[k + 2] = V[3 * b + 2];
144 }
145 };
146 syncAnchors();
148 const force = exc.force;
149 const fShape = force?.shape;
150 const pulse = (t: number) =>
151 force && t >= 0 && t <= force.tau ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau)) : 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);
157 let fContact = 0;
158 let eContact = 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];
172 const contact = () => {
173 fContact = 0;
174 eContact = 0;
175 pFriction = 0;
176 qa.fill(0);
177 if (!mallet) return;
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) {
181 const k = 3 * i;
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];
188 }
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;
224 }
225 };
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.
231 let pDamp = 0;
232 let pExt = 0;
233 const computeForces = (t: number) => {
234 F.fill(0);
235 contact();
236 pDamp = pFriction;
237 for (let s = 0; s < ns; s++) {
238 const a = 3 * si[s];
239 const b = 3 * sj[s];
240 const dx = X[b] - X[a];
241 const dy = X[b + 1] - X[a + 1];
242 const dz = X[b + 2] - X[a + 2];
243 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
244 const g = (sk[s] * (l - sL0[s])) / l;
245 let fx = g * dx;
246 let fy = g * dy;
247 let fz = g * dz;
248 const c = sc[s];
249 if (c !== 0) {
250 const rx = V[b] - V[a];
251 const ry = V[b + 1] - V[a + 1];
252 const rz = V[b + 2] - V[a + 2];
253 fx += c * rx;
254 fy += c * ry;
255 fz += c * rz;
256 pDamp += c * (rx * rx + ry * ry + rz * rz);
257 }
258 if (si[s] < n) {
259 F[a] += fx;
260 F[a + 1] += fy;
261 F[a + 2] += fz;
262 }
263 if (sj[s] < n) {
264 F[b] -= fx;
265 F[b + 1] -= fy;
266 F[b + 2] -= fz;
267 }
268 }
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 }
330 const s0 = sys.sigma0;
331 const p = fShape ? pulse(t) : 0;
332 pExt = 0;
333 for (let i = 0; i < n; i++) {
334 const k = 3 * i;
335 if (pinned[i]) {
336 F[k] = F[k + 1] = F[k + 2] = 0;
337 continue;
338 }
339 if (s0 !== 0) {
340 const m2 = 2 * s0 * mass[i];
341 for (let c = 0; c < 3; c++) {
342 F[k + c] -= m2 * V[k + c];
343 pDamp += m2 * V[k + c] * V[k + c];
344 }
345 }
346 if (p !== 0) {
347 for (let c = 0; c < 3; c++) {
348 const fe = p * fShape![k + c];
349 F[k + c] += fe;
350 pExt += fe * V[k + c];
351 }
352 }
353 }
354 for (let i = 0; i < n; i++) {
355 const im = 1 / mass[i];
356 A[3 * i] = F[3 * i] * im;
357 A[3 * i + 1] = F[3 * i + 1] * im;
358 A[3 * i + 2] = F[3 * i + 2] * im;
359 }
360 };
362 const velocityForces =
363 sys.sigma0 !== 0 || sc.some((c) => c !== 0) || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
364 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)
368 const energyNow = () => {
369 let e = 0;
370 for (let i = 0; i < n; i++) {
371 const k = 3 * i;
372 e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
373 }
374 for (let s = 0; s < ns; s++) {
375 const a = 3 * si[s];
376 const b = 3 * sj[s];
377 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
378 const d = l - sL0[s];
379 e += 0.5 * sk[s] * d * d;
380 if (sj[s] < n && si[s] < n) maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
381 }
382 return e + eHinge + malletEnergy();
383 };
385 const audio = new Float32Array(nSamples);
386 const frames = new Float32Array(nFrames * 3 * n);
387 const frameTimes = new Float64Array(nFrames);
388 const energy = new Float64Array(nFrames);
389 const balance = new Float64Array(nFrames);
390 const contactForce = mallet ? new Float32Array(nSamples) : undefined;
391 const malletPos = mallet ? new Float32Array(3 * nFrames) : undefined;
392 const w = opt.pickup;
393 const delay = opt.pickupDelay;
395 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
396 let Eeq = 0;
397 for (let s = 0; s < ns; s++) {
398 const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
399 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));
400 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
401 }
402 computeForces(0);
403 const E0 = energyNow();
404 let dissipated = 0;
405 let workIn = 0;
406 let t = 0;
407 let frame = 0;
408 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
410 for (let smp = 0; smp < nSamples; smp++) {
411 // record at t = smp / fs; A holds forces at full-step velocity here
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 }
428 if (contactForce) contactForce[smp] = fContact;
429 if (smp % samplesPerFrame === 0 && frame < nFrames) {
430 if (malletPos) malletPos.set(q, 3 * frame);
431 const off = frame * 3 * n;
432 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
433 frameTimes[frame] = t;
434 const eTotal = energyNow();
435 // the plot shows the instrument's own energy; the balance includes the mallet
436 energy[frame] = eTotal - Eeq - malletEnergy();
437 balance[frame] = eTotal + dissipated - workIn - E0;
438 frame++;
439 }
440 if (smp === nSamples - 1) break;
442 for (let st = 0; st < stepsPerSample; st++) {
443 // kick (half), drift
444 for (let k = 0; k < 3 * n; k++) {
445 V[k] += 0.5 * dt * A[k];
446 X[k] += dt * V[k];
447 }
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 }
453 syncAnchors();
454 t += dt;
455 // forces at new position with half-step velocity; powers use it too
456 computeForces(t);
457 dissipated += dt * pDamp;
458 workIn += dt * pExt;
459 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
460 if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];
461 syncAnchors();
462 }
463 // at a sample time: recompute with the full-step velocity
464 if (velocityForces) computeForces(t);
465 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
466 }
468 return {
469 audio,
470 fs: opt.fs,
471 frames,
472 frameTimes,
473 energy,
474 balance,
475 dt,
476 stepsPerSample,
477 omegaMax,
478 wallSeconds: (performance.now() - t0) / 1000,
479 maxStretch,
480 contactForce,
481 malletPos,
482 };
483}