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// 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 {
24 const { n, si, sj, sk, sL0 } = sys;
25 const X = nodePositions(sys);
26 const ns = si.length;
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];
39 }
40 const b = sys.bending;
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) {
45 const g = new Float64Array(13);
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);
50 }
51 }
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;
135 }
136 return Math.sqrt(w2);
137}
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 };
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}
169export type Progress = (fraction: number) => void;
170// receives the result so far (a prefix of every output) while a run is in progress
171export type PartialResult = (r: SimResult) => void;
173// wall time before the first partial result, and between later ones (ms)
174export const PARTIAL_FIRST_MS = 1000;
175export const PARTIAL_EVERY_MS = 2000;
177// Velocity Verlet (kick-drift-kick), explicit. Spring and hinge forces are
178// exact (geometric nonlinearity kept); damping and friction use the half-step
179// velocity in the kick, and the full-step velocity whenever an audio sample is
180// recorded.
181export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, progress?: Progress, partial?: PartialResult): SimResult {
182 const t0 = performance.now();
183 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
184 const nn = n + anchors.length;
185 const ns = si.length;
186 const bend = sys.bending;
187 const nh = bend ? bend.w0.length : 0;
188 const nt = bend ? bend.tri.length / 3 : 0;
189 const hG = new Float64Array(12 * nh); // gradients of phi
190 const hPhi = new Float64Array(nh);
191 const hRate = new Float64Array(nh);
192 const hMom = new Float64Array(nh);
193 const hMomEl = new Float64Array(nh); // elastic part of the moments
195 const mallet = exc.mallet;
196 const { omegaMax, dt, stepsPerSample } = timeStep(sys, mallet, opt.fs, opt.cfl);
197 const nSamples = Math.round(opt.duration * opt.fs);
198 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
199 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
201 const X = new Float64Array(3 * nn);
202 const V = new Float64Array(3 * nn);
203 const F = new Float64Array(3 * n); // force on beads
204 const Fel = new Float64Array(3 * n); // elastic (position-dependent) part
205 const A = new Float64Array(3 * n);
206 for (let i = 0; i < n; i++) {
207 for (let c = 0; c < 3; c++) {
208 const k = 3 * i + c;
209 X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
210 V[k] = pinned[i] ? 0 : exc.v0[k];
211 }
212 }
213 const slide: number[] = [];
214 for (let a = 0; a < anchors.length; a++) {
215 for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
216 if (anchors[a].slideBead >= 0) slide.push(a);
217 }
218 const syncAnchors = () => {
219 for (const a of slide) {
220 const b = anchors[a].slideBead;
221 const k = 3 * (n + a);
222 X[k + 1] = X[3 * b + 1];
223 X[k + 2] = X[3 * b + 2];
224 V[k + 1] = V[3 * b + 1];
225 V[k + 2] = V[3 * b + 2];
226 }
227 };
228 syncAnchors();
230 const force = exc.force;
231 const fShape = force?.shape;
232 const pulse = (t: number) =>
233 force && t >= 0 && t <= force.tau ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau)) : 0;
235 // mallet state: 3D position q, velocity qv, acceleration qa
236 const q = new Float64Array(3);
237 const qv = new Float64Array(mallet ? mallet.vel0 : [0, 0, 0]);
238 const qa = new Float64Array(3);
239 let fContact = 0;
240 let eContact = 0;
241 let pFriction = 0;
242 const patch: number[] = [];
243 let mEff = 1;
244 if (mallet) {
245 let inv = 0;
246 for (let i = 0; i < n; i++)
247 if (mallet.weights[i] > 0) {
248 patch.push(i);
249 inv += (mallet.weights[i] * mallet.weights[i]) / mass[i];
250 }
251 mEff = 1 / (1 / mallet.mass + inv);
252 }
253 const [nx, ny, nz] = mallet?.normal ?? [0, 0, 0];
254 const contact = () => {
255 fContact = 0;
256 eContact = 0;
257 pFriction = 0;
258 qa.fill(0);
259 if (!mallet) return;
260 const w = mallet.weights;
261 let ux = 0, uy = 0, uz = 0, vx = 0, vy = 0, vz = 0;
262 for (const i of patch) {
263 const k = 3 * i;
264 ux += w[i] * (X[k] - sys.rest[k]);
265 uy += w[i] * (X[k + 1] - sys.rest[k + 1]);
266 uz += w[i] * (X[k + 2] - sys.rest[k + 2]);
267 vx += w[i] * V[k];
268 vy += w[i] * V[k + 1];
269 vz += w[i] * V[k + 2];
270 }
271 const d = (q[0] - ux) * nx + (q[1] - uy) * ny + (q[2] - uz) * nz;
272 if (!(d > 0)) return;
273 if (mallet.reach !== undefined) {
274 const ox = q[0] - ux - d * nx, oy = q[1] - uy - d * ny, oz = q[2] - uz - d * nz;
275 if (ox * ox + oy * oy + oz * oz > mallet.reach * mallet.reach) return;
276 }
277 const fn = mallet.K * d ** mallet.p;
278 fContact = fn;
279 eContact = (fn * d) / (mallet.p + 1);
280 // force on the mallet: -fn normal + friction against tangential slip
281 let fx = -fn * nx;
282 let fy = -fn * ny;
283 let fz = -fn * nz;
284 if (mallet.mu > 0) {
285 let rx = qv[0] - vx, ry = qv[1] - vy, rz = qv[2] - vz;
286 const rn = rx * nx + ry * ny + rz * nz;
287 rx -= rn * nx;
288 ry -= rn * ny;
289 rz -= rn * nz;
290 const slip = Math.hypot(rx, ry, rz);
291 // regularized Coulomb; the slip scale keeps the explicit step stable
292 const eps = Math.max(1e-3, (4 * mallet.mu * fn * dt) / mEff);
293 const g = (mallet.mu * fn) / Math.sqrt(slip * slip + eps * eps);
294 fx -= g * rx;
295 fy -= g * ry;
296 fz -= g * rz;
297 pFriction = g * slip * slip;
298 }
299 qa[0] = fx / mallet.mass;
300 qa[1] = fy / mallet.mass;
301 qa[2] = fz / mallet.mass;
302 for (const i of patch) {
303 F[3 * i] -= w[i] * fx;
304 F[3 * i + 1] -= w[i] * fy;
305 F[3 * i + 2] -= w[i] * fz;
306 }
307 };
309 const hOut = new Float64Array(13);
310 let eHinge = 0;
311 const anyDashpot = sc.some((c) => c !== 0);
312 // Position-dependent forces: spring and bending elasticity into Fel; hinge
313 // angles, gradients and elastic moments; rigid-body geometry.
314 const elasticForces = () => {
315 Fel.fill(0);
316 for (let s = 0; s < ns; s++) {
317 const a = 3 * si[s];
318 const b = 3 * sj[s];
319 const dx = X[b] - X[a];
320 const dy = X[b + 1] - X[a + 1];
321 const dz = X[b + 2] - X[a + 2];
322 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
323 const g = (sk[s] * (l - sL0[s])) / l;
324 if (si[s] < n) {
325 Fel[a] += g * dx;
326 Fel[a + 1] += g * dy;
327 Fel[a + 2] += g * dz;
328 }
329 if (sj[s] < n) {
330 Fel[b] -= g * dx;
331 Fel[b + 1] -= g * dy;
332 Fel[b + 2] -= g * dz;
333 }
334 }
335 eHinge = 0;
336 if (bend) {
337 for (let e = 0; e < nh; e++) {
338 const th = hingeAngleGrad(X, bend.w0[e], bend.w1[e], bend.e0[e], bend.e1[e], hOut);
339 const sg = bend.sign[e];
340 hPhi[e] = sg * (th - bend.theta0[e]);
341 const o = 12 * e;
342 for (let q = 0; q < 12; q++) hG[o + q] = sg * hOut[1 + q];
343 }
344 hMomEl.fill(0);
345 const tri = bend.tri;
346 const Mt = bend.M;
347 for (let t = 0; t < nt; t++) {
348 for (let a = 0; a < 3; a++) {
349 const ea = tri[3 * t + a];
350 if (ea < 0) continue;
351 let m = 0;
352 for (let c = 0; c < 3; c++) {
353 const ec = tri[3 * t + c];
354 if (ec >= 0) m += Mt[9 * t + 3 * a + c] * hPhi[ec];
355 }
356 eHinge += 0.5 * hPhi[ea] * m;
357 hMomEl[ea] += m;
358 }
359 }
360 }
361 if (sys.freeRigid) rigidGeometry();
362 };
364 // Fills F with total force on beads and A with accelerations, from the
365 // elastic part (elasticForces must have run at this position) and the
366 // velocity-dependent part at V; sets damping power (W, >= 0) and external
367 // power (W), both at V.
368 let pDamp = 0;
369 let pExt = 0;
370 const dynamicForces = (t: number) => {
371 F.set(Fel);
372 contact();
373 pDamp = pFriction;
374 if (anyDashpot)
375 for (let s = 0; s < ns; s++) {
376 const c = sc[s];
377 if (c === 0) continue;
378 const a = 3 * si[s];
379 const b = 3 * sj[s];
380 const rx = V[b] - V[a];
381 const ry = V[b + 1] - V[a + 1];
382 const rz = V[b + 2] - V[a + 2];
383 pDamp += c * (rx * rx + ry * ry + rz * rz);
384 if (si[s] < n) {
385 F[a] += c * rx;
386 F[a + 1] += c * ry;
387 F[a + 2] += c * rz;
388 }
389 if (sj[s] < n) {
390 F[b] -= c * rx;
391 F[b + 1] -= c * ry;
392 F[b + 2] -= c * rz;
393 }
394 }
395 if (bend) {
396 const s1 = bend.sigma1;
397 hMom.set(hMomEl);
398 if (s1 !== 0) {
399 for (let e = 0; e < nh; e++) {
400 const o = 12 * e;
401 const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
402 hRate[e] =
403 hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
404 hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
405 hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
406 hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
407 }
408 const tri = bend.tri;
409 const Mt = bend.M;
410 for (let t = 0; t < nt; t++) {
411 for (let a = 0; a < 3; a++) {
412 const ea = tri[3 * t + a];
413 if (ea < 0) continue;
414 let r = 0;
415 for (let c = 0; c < 3; c++) {
416 const ec = tri[3 * t + c];
417 if (ec >= 0) r += Mt[9 * t + 3 * a + c] * hRate[ec];
418 }
419 pDamp += 2 * s1 * hRate[ea] * r;
420 hMom[ea] += 2 * s1 * r;
421 }
422 }
423 }
424 for (let e = 0; e < nh; e++) {
425 const m = hMom[e];
426 if (m === 0) continue;
427 const o = 12 * e;
428 const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
429 F[w0] -= m * hG[o];
430 F[w0 + 1] -= m * hG[o + 1];
431 F[w0 + 2] -= m * hG[o + 2];
432 F[w1] -= m * hG[o + 3];
433 F[w1 + 1] -= m * hG[o + 4];
434 F[w1 + 2] -= m * hG[o + 5];
435 F[e0] -= m * hG[o + 6];
436 F[e0 + 1] -= m * hG[o + 7];
437 F[e0 + 2] -= m * hG[o + 8];
438 F[e1] -= m * hG[o + 9];
439 F[e1 + 1] -= m * hG[o + 10];
440 F[e1 + 2] -= m * hG[o + 11];
441 }
442 }
443 const s0 = sys.sigma0;
444 const p = fShape ? pulse(t) : 0;
445 pExt = 0;
446 for (let i = 0; i < n; i++) {
447 const k = 3 * i;
448 if (pinned[i]) {
449 F[k] = F[k + 1] = F[k + 2] = 0;
450 continue;
451 }
452 if (s0 !== 0) {
453 const m2 = 2 * s0 * mass[i];
454 for (let c = 0; c < 3; c++) {
455 F[k + c] -= m2 * V[k + c];
456 pDamp += m2 * V[k + c] * V[k + c];
457 }
458 }
459 if (p !== 0) {
460 for (let c = 0; c < 3; c++) {
461 const fe = p * fShape![k + c];
462 F[k + c] += fe;
463 pExt += fe * V[k + c];
464 }
465 }
466 }
467 for (let i = 0; i < n; i++) {
468 const im = 1 / mass[i];
469 A[3 * i] = F[3 * i] * im;
470 A[3 * i + 1] = F[3 * i + 1] * im;
471 A[3 * i + 2] = F[3 * i + 2] * im;
472 }
473 if (sys.freeRigid) removeRigid(A);
474 };
475 const computeForces = (t: number) => {
476 elasticForces();
477 dynamicForces(t);
478 };
480 // centre of mass and inverse inertia tensor at the current positions
481 let totalMass = 0;
482 for (let i = 0; i < n; i++) totalMass += mass[i];
483 const cm = new Float64Array(3);
484 const Iinv = new Float64Array(6);
485 const rigidGeometry = () => {
486 let cx = 0, cy = 0, cz = 0;
487 for (let i = 0; i < n; i++) {
488 cx += mass[i] * X[3 * i];
489 cy += mass[i] * X[3 * i + 1];
490 cz += mass[i] * X[3 * i + 2];
491 }
492 cx /= totalMass;
493 cy /= totalMass;
494 cz /= totalMass;
495 let ixx = 0, iyy = 0, izz = 0, ixy = 0, ixz = 0, iyz = 0;
496 for (let i = 0; i < n; i++) {
497 const m = mass[i];
498 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
499 ixx += m * (ry * ry + rz * rz);
500 iyy += m * (rx * rx + rz * rz);
501 izz += m * (rx * rx + ry * ry);
502 ixy -= m * rx * ry;
503 ixz -= m * rx * rz;
504 iyz -= m * ry * rz;
505 }
506 cm[0] = cx;
507 cm[1] = cy;
508 cm[2] = cz;
509 invert3(ixx, ixy, ixz, iyy, iyz, izz, Iinv);
510 };
512 // subtract the rigid-body part (mass-weighted) of a per-bead field u, about
513 // the geometry from rigidGeometry; returns the kinetic energy of that part
514 // when u is a velocity
515 const removeRigid = (u: Float64Array): number => {
516 const [cx, cy, cz] = cm;
517 let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;
518 for (let i = 0; i < n; i++) {
519 const m = mass[i];
520 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
521 const ux = u[3 * i], uy = u[3 * i + 1], uz = u[3 * i + 2];
522 px += m * ux;
523 py += m * uy;
524 pz += m * uz;
525 lx += m * (ry * uz - rz * uy);
526 ly += m * (rz * ux - rx * uz);
527 lz += m * (rx * uy - ry * ux);
528 }
529 const [a, b, c, d, e, f] = Iinv;
530 const wx = a * lx + b * ly + c * lz;
531 const wy = b * lx + d * ly + e * lz;
532 const wz = c * lx + e * ly + f * lz;
533 const vx = px / totalMass, vy = py / totalMass, vz = pz / totalMass;
534 for (let i = 0; i < n; i++) {
535 const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
536 u[3 * i] -= vx + (wy * rz - wz * ry);
537 u[3 * i + 1] -= vy + (wz * rx - wx * rz);
538 u[3 * i + 2] -= vz + (wx * ry - wy * rx);
539 }
540 return 0.5 * (px * vx + py * vy + pz * vz) + 0.5 * (wx * lx + wy * ly + wz * lz);
541 };
542 if (sys.freeRigid) {
543 rigidGeometry();
544 removeRigid(V); // an initial strike carries no rigid motion either
545 }
547 const velocityForces = sys.sigma0 !== 0 || anyDashpot || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
548 let maxStretch = 0;
549 const malletEnergy = () => (mallet ? 0.5 * mallet.mass * (qv[0] ** 2 + qv[1] ** 2 + qv[2] ** 2) + eContact : 0);
550 // total mechanical energy at the current state; computeForces must have run
551 // at this state (it sets the hinge and contact energies)
552 const energyNow = () => {
553 let e = 0;
554 for (let i = 0; i < n; i++) {
555 const k = 3 * i;
556 e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
557 }
558 for (let s = 0; s < ns; s++) {
559 const a = 3 * si[s];
560 const b = 3 * sj[s];
561 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
562 const d = l - sL0[s];
563 e += 0.5 * sk[s] * d * d;
564 if (sj[s] < n && si[s] < n) maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
565 }
566 return e + eHinge + malletEnergy();
567 };
569 const audio = new Float32Array(nSamples);
570 const frames = new Float32Array(nFrames * 3 * n);
571 const frameTimes = new Float64Array(nFrames);
572 const energy = new Float64Array(nFrames);
573 const balance = new Float64Array(nFrames);
574 const contactForce = mallet ? new Float32Array(nSamples) : undefined;
575 const malletPos = mallet ? new Float32Array(3 * nFrames) : undefined;
576 const w = opt.pickup;
577 const delay = opt.pickupDelay;
579 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
580 let Eeq = 0;
581 for (let s = 0; s < ns; s++) {
582 const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
583 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));
584 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
585 }
586 computeForces(0);
587 const E0 = energyNow();
588 let dissipated = 0;
589 let workIn = 0;
590 let t = 0;
591 let frame = 0;
592 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
593 let nextPartial = t0 + PARTIAL_FIRST_MS;
594 // The outputs up to sample count S and frame count F: copies, except for the
595 // finished run (a partial result may be transferred away). With a delayed
596 // pickup, audio[i] only gets terms from samples <= i, so the first S are complete.
597 const result = (S: number, F: number): SimResult => {
598 const all = S === nSamples;
599 return {
600 audio: all ? audio : audio.slice(0, S),
601 fs: opt.fs,
602 frames: all ? frames : frames.slice(0, F * 3 * n),
603 frameTimes: all ? frameTimes : frameTimes.slice(0, F),
604 energy: all ? energy : energy.slice(0, F),
605 balance: all ? balance : balance.slice(0, F),
606 dt,
607 stepsPerSample,
608 omegaMax,
609 wallSeconds: (performance.now() - t0) / 1000,
610 maxStretch,
611 contactForce: all ? contactForce : contactForce?.slice(0, S),
612 malletPos: all ? malletPos : malletPos?.slice(0, 3 * F),
613 };
614 };
616 for (let smp = 0; smp < nSamples; smp++) {
617 // record at t = smp / fs; A holds forces at full-step velocity here
618 if (delay) {
619 for (let i = 0; i < n; i++) {
620 const k = 3 * i;
621 const c = w[k] * A[k] + w[k + 1] * A[k + 1] + w[k + 2] * A[k + 2];
622 if (c === 0) continue;
623 const pos = smp + delay[i];
624 const i0 = Math.floor(pos);
625 const fr = pos - i0;
626 if (i0 < nSamples) audio[i0] += (1 - fr) * c;
627 if (i0 + 1 < nSamples) audio[i0 + 1] += fr * c;
628 }
629 } else {
630 let acc = 0;
631 for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
632 audio[smp] = acc;
633 }
634 if (contactForce) contactForce[smp] = fContact;
635 if (smp % samplesPerFrame === 0 && frame < nFrames) {
636 if (malletPos) malletPos.set(q, 3 * frame);
637 const off = frame * 3 * n;
638 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
639 frameTimes[frame] = t;
640 const eTotal = energyNow();
641 // the plot shows the instrument's own energy; the balance includes the mallet
642 energy[frame] = eTotal - Eeq - malletEnergy();
643 balance[frame] = eTotal + dissipated - workIn - E0;
644 frame++;
645 }
646 if (smp === nSamples - 1) break;
647 if (partial && smp % 32 === 0 && performance.now() >= nextPartial) {
648 partial(result(smp + 1, frame));
649 nextPartial = performance.now() + PARTIAL_EVERY_MS;
650 }
652 for (let st = 0; st < stepsPerSample; st++) {
653 // kick (half), drift
654 for (let k = 0; k < 3 * n; k++) {
655 V[k] += 0.5 * dt * A[k];
656 X[k] += dt * V[k];
657 }
658 if (mallet)
659 for (let c = 0; c < 3; c++) {
660 qv[c] += 0.5 * dt * qa[c];
661 q[c] += dt * qv[c];
662 }
663 syncAnchors();
664 t += dt;
665 // forces at new position with half-step velocity; powers use it too
666 computeForces(t);
667 dissipated += dt * pDamp;
668 workIn += dt * pExt;
669 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
670 if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];
671 // the rigid-free subspace turns with the body; drop what drifts out of it
672 if (sys.freeRigid) dissipated += removeRigid(V);
673 syncAnchors();
674 }
675 // at a sample time: recompute with the full-step velocity (positions unchanged)
676 if (velocityForces) dynamicForces(t);
677 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
678 }
680 return result(nSamples, nFrames);
681}
683// inverse of the symmetric 3x3 matrix [a b c; b d e; c e f], as its six
684// distinct entries (zero if singular)
685function invert3(a: number, b: number, c: number, d: number, e: number, f: number, out: Float64Array) {
686 const A00 = d * f - e * e, A01 = c * e - b * f, A02 = b * e - c * d;
687 const A11 = a * f - c * c, A12 = b * c - a * e, A22 = a * d - b * b;
688 const det = a * A00 + b * A01 + c * A02;
689 if (!(Math.abs(det) > 0)) {
690 out.fill(0);
691 return;
692 }
693 out[0] = A00 / det;
694 out[1] = A01 / det;
695 out[2] = A02 / det;
696 out[3] = A11 / det;
697 out[4] = A12 / det;
698 out[5] = A22 / det;
699}