/ proof-of-concept / toy-instruments
Sign in
proof-of-concept / toy-instruments
2cd36fa
Branches
main
toy-instruments / src / sim / solver.ts
310 lines · 9.8 KBCodeBlameHistory
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 1import type { BeadSystem, Excitation, Mallet, SimOptions, SimResult } from './types';
3// Upper bound on the highest angular frequency of the linearized system
4// (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
5// its axis and T/l across it, so each spring contributes 2 max(k, |T|/l) / m.
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 6export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
8 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];
11 for (let s = 0; s < si.length; s++) {
12 const a = si[s];
13 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);
17 const T = sk[s] * (l - sL0[s]);
18 const kmax = Math.max(sk[s], Math.abs(T) / l);
19 if (a < n) rowSum[a] += 2 * kmax;
20 if (b < n) rowSum[b] += 2 * kmax;
21 }
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 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);
29 }
30 let w2 = mallet ? kc / mallet.mass : 0;
31 for (let i = 0; i < n; i++) {
32 if (sys.pinned[i]) continue;
33 const extra = mallet ? kc * mallet.weights[i] : 0;
34 w2 = Math.max(w2, (rowSum[i] + extra) / sys.mass[i]);
35 }
39export type Progress = (fraction: number) => void;
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.
44export function simulate(
45 sys: BeadSystem,
46 exc: Excitation,
47 opt: SimOptions,
48 progress?: Progress,
49): SimResult {
50 const t0 = performance.now();
51 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
52 const nn = n + anchors.length;
53 const ns = si.length;
56 const omegaMax = omegaMaxBound(sys, mallet);
58 const stepsPerSample = Math.max(1, Math.ceil(1 / (opt.fs * dtMax)));
59 const dt = 1 / (opt.fs * stepsPerSample);
60 const nSamples = Math.round(opt.duration * opt.fs);
61 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
62 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
64 const X = new Float64Array(3 * nn);
65 const V = new Float64Array(3 * nn);
66 const F = new Float64Array(3 * n); // force on beads
67 const A = new Float64Array(3 * n);
68 for (let i = 0; i < n; i++) {
69 for (let c = 0; c < 3; c++) {
70 const k = 3 * i + c;
71 X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
72 V[k] = pinned[i] ? 0 : exc.v0[k];
73 }
74 }
75 const slide: number[] = [];
76 for (let a = 0; a < anchors.length; a++) {
77 for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
78 if (anchors[a].slideBead >= 0) slide.push(a);
79 }
80 const syncAnchors = () => {
81 for (const a of slide) {
82 const b = anchors[a].slideBead;
83 const k = 3 * (n + a);
84 X[k + 1] = X[3 * b + 1];
85 X[k + 2] = X[3 * b + 2];
86 V[k + 1] = V[3 * b + 1];
87 V[k + 2] = V[3 * b + 2];
88 }
89 };
90 syncAnchors();
92 const force = exc.force;
93 const fShape = force?.shape;
94 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;
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 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;
103 let fContact = 0;
104 let eContact = 0;
105 const mw = mallet?.weights;
106 const [mdx, mdy, mdz] = mallet?.dir ?? [0, 0, 0];
107 const contact = () => {
108 fContact = 0;
109 eContact = 0;
110 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;
115 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);
117 }
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 }
130 }
131 ma = -fContact / mallet.mass;
132 };
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 134 // Fills F with total force on beads; returns damping power (W, >= 0)
135 // and external power (W), both at velocity V.
136 let pDamp = 0;
137 let pExt = 0;
138 const computeForces = (t: number) => {
139 F.fill(0);
142 for (let s = 0; s < ns; s++) {
143 const a = 3 * si[s];
144 const b = 3 * sj[s];
145 const dx = X[b] - X[a];
146 const dy = X[b + 1] - X[a + 1];
147 const dz = X[b + 2] - X[a + 2];
148 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
149 const g = (sk[s] * (l - sL0[s])) / l;
150 let fx = g * dx;
151 let fy = g * dy;
152 let fz = g * dz;
153 const c = sc[s];
154 if (c !== 0) {
155 const rx = V[b] - V[a];
156 const ry = V[b + 1] - V[a + 1];
157 const rz = V[b + 2] - V[a + 2];
158 fx += c * rx;
159 fy += c * ry;
160 fz += c * rz;
161 pDamp += c * (rx * rx + ry * ry + rz * rz);
162 }
163 if (si[s] < n) {
164 F[a] += fx;
165 F[a + 1] += fy;
166 F[a + 2] += fz;
167 }
168 if (sj[s] < n) {
169 F[b] -= fx;
170 F[b + 1] -= fy;
171 F[b + 2] -= fz;
172 }
173 }
174 const s0 = sys.sigma0;
175 const p = fShape ? pulse(t) : 0;
176 pExt = 0;
177 for (let i = 0; i < n; i++) {
178 const k = 3 * i;
179 if (pinned[i]) {
180 F[k] = F[k + 1] = F[k + 2] = 0;
181 continue;
182 }
183 if (s0 !== 0) {
184 const m2 = 2 * s0 * mass[i];
185 for (let c = 0; c < 3; c++) {
186 F[k + c] -= m2 * V[k + c];
187 pDamp += m2 * V[k + c] * V[k + c];
188 }
189 }
190 if (p !== 0) {
191 for (let c = 0; c < 3; c++) {
192 const fe = p * fShape![k + c];
193 F[k + c] += fe;
194 pExt += fe * V[k + c];
195 }
196 }
197 }
198 for (let i = 0; i < n; i++) {
199 const im = 1 / mass[i];
200 A[3 * i] = F[3 * i] * im;
201 A[3 * i + 1] = F[3 * i + 1] * im;
202 A[3 * i + 2] = F[3 * i + 2] * im;
203 }
204 };
206 const damped = sys.sigma0 !== 0 || sc.some((c) => c !== 0);
207 let maxStretch = 0;
208 const energyNow = () => {
209 let e = 0;
210 for (let i = 0; i < n; i++) {
211 const k = 3 * i;
212 e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
213 }
214 for (let s = 0; s < ns; s++) {
215 const a = 3 * si[s];
216 const b = 3 * sj[s];
217 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
218 const d = l - sL0[s];
219 e += 0.5 * sk[s] * d * d;
220 maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
221 }
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 222 if (mallet) e += 0.5 * mallet.mass * mv * mv + eContact;
224 };
226 const audio = new Float32Array(nSamples);
227 const frames = new Float32Array(nFrames * 3 * n);
228 const frameTimes = new Float64Array(nFrames);
229 const energy = new Float64Array(nFrames);
230 const balance = new Float64Array(nFrames);
2cd36faBead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument pageJeremy Magland 231 const contactForce = mallet ? new Float32Array(nSamples) : undefined;
232 const malletPos = mallet ? new Float32Array(nFrames) : undefined;
235 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
236 let Eeq = 0;
237 for (let s = 0; s < ns; s++) {
238 const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
239 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));
240 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
241 }
244 let dissipated = 0;
245 let workIn = 0;
246 let t = 0;
247 computeForces(0);
248 let frame = 0;
249 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
251 for (let smp = 0; smp < nSamples; smp++) {
252 // 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;
4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 257 if (smp % samplesPerFrame === 0 && frame < nFrames) {
260 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
261 frameTimes[frame] = t;
263 // 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);
265 balance[frame] = eTotal + dissipated - workIn - E0;
267 }
268 if (smp === nSamples - 1) break;
270 for (let st = 0; st < stepsPerSample; st++) {
271 // kick (half), drift
272 for (let k = 0; k < 3 * n; k++) {
273 V[k] += 0.5 * dt * A[k];
274 X[k] += dt * V[k];
275 }
277 mv += 0.5 * dt * ma;
278 ms += dt * mv;
279 }
281 t += dt;
282 // forces at new position with half-step velocity; powers use it too
283 computeForces(t);
284 dissipated += dt * pDamp;
285 workIn += dt * pExt;
286 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
289 }
290 // at a sample time: recompute with the full-step velocity
291 if (damped) computeForces(t);
292 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
293 }
295 return {
296 audio,
297 fs: opt.fs,
298 frames,
299 frameTimes,
300 energy,
301 balance,
302 dt,
303 stepsPerSample,
304 omegaMax,
305 wallSeconds: (performance.now() - t0) / 1000,
306 maxStretch,
308 malletPos,
↑↓ move↵openescclose