/ proof-of-concept / bead-percussion
Sign in
proof-of-concept / bead-percussion
74b1dcd
Branches
main
bead-percussion / src / sim / modes.ts
99 lines · 4.1 KBCodeBlameHistory
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 1// Linear transverse modes of a flat bead system.
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 3// When every bead rests in a plane z = const and every spring lies in that
4// plane, small transverse (z) motion decouples from in-plane motion, and its
5// stiffness matrix is a graph Laplacian weighted by spring tension over length:
7// triplet, kb (e_i - 2 e_j + e_k)(e_i - 2 e_j + e_k)^T (the triplet energy
8// separates by component, so its z part decouples too). We return the
9// eigenpairs of M^{-1/2} K M^{-1/2} over the free beads.
11import type { BeadSystem, Modes } from "./types";
13export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
15export function isPlanar(sys: BeadSystem): boolean {
16 const { n, pos0 } = sys;
17 let zmin = Infinity;
18 let zmax = -Infinity;
19 let span = 0;
20 for (let i = 0; i < n; i++) {
21 zmin = Math.min(zmin, pos0[3 * i + 2]);
22 zmax = Math.max(zmax, pos0[3 * i + 2]);
23 span = Math.max(span, Math.abs(pos0[3 * i]), Math.abs(pos0[3 * i + 1]));
24 }
25 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 28export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
29 const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
30 const freeIdx = new Int32Array(n).fill(-1);
31 const free: number[] = [];
32 for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
33 const nf = free.length;
34 const K = new Float64Array(nf * nf);
36 const i3 = 3 * si[s];
37 const j3 = 3 * sj[s];
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 38 const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
39 const e = L - sL0[s];
40 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
41 const a = freeIdx[si[s]];
42 const b = freeIdx[sj[s]];
43 if (a >= 0) K[a * nf + a] += g;
44 if (b >= 0) K[b * nf + b] += g;
45 if (a >= 0 && b >= 0) {
46 K[a * nf + b] -= g;
47 K[b * nf + a] -= g;
49 }
51 const ids = [freeIdx[sys.bi[t]], freeIdx[sys.bj[t]], freeIdx[sys.bk[t]]];
52 const c = [1, -2, 1];
53 for (let p = 0; p < 3; p++) {
54 if (ids[p] < 0) continue;
55 for (let q = 0; q < 3; q++) if (ids[q] >= 0) K[ids[p] * nf + ids[q]] += sys.bkb[t] * c[p] * c[q];
56 }
57 }
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 58 for (let a = 0; a < nf; a++) {
59 for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 61 return { K, free: Int32Array.from(free) };
64export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 65 const { K, free } = transverseStiffness(sys);
66 const nf = free.length;
67 const { values, vectors } = eig(K, nf);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 68 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 69 const freqs = new Float64Array(nf);
70 const psi = new Float64Array(nf * nf);
72 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 73 psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 75 return { freqs, psi, free };
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 78// Small-amplitude decay rate of each transverse mode (amplitude, 1/s), from
79// the friction and dashpots: sigma_k = (1/2) u^T C u with u the mode shape
80// normalized so u^T M u = 1. Exact when the damping matrix is a combination
81// of M and K (as bp_decay builds it); otherwise it neglects the coupling of
0330c41Damping by decay times: bp_decay helper, per-mode decay-time plot, examples convertedJeremy Magland 82// modes by damping, which is small for light damping.
83export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
84 const { n, mass, friction, di, dj, dc } = sys;
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 85 const { free, psi, freqs } = modes;
86 const nf = free.length;
87 const u = new Float64Array(n);
0330c41Damping by decay times: bp_decay helper, per-mode decay-time plot, examples convertedJeremy Magland 88 const out = new Float64Array(freqs.length);
89 for (let k = 0; k < freqs.length; k++) {
90 u.fill(0);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 91 for (let r = 0; r < nf; r++) u[free[r]] = psi[k * nf + r] / Math.sqrt(mass[free[r]]);
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 93 for (let r = 0; r < nf; r++) c += friction[free[r]] * u[free[r]] ** 2;
94 for (let d = 0; d < di.length; d++) c += dc[d] * (u[di[d]] - u[dj[d]]) ** 2;
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 95 for (let t = 0; t < sys.bi.length; t++) c += sys.bcb[t] * (u[sys.bi[t]] - 2 * u[sys.bj[t]] + u[sys.bk[t]]) ** 2;
97 }
98 return out;
↑↓ move↵openescclose