1// Linear modes of a bead system about its rest positions.
2//
3// The stiffness at rest is assembled from the elements: a spring contributes
4// the 3x3 block keff d d^T + (f/L)(I - d d^T) (d its direction, f its tension,
5// keff = k + 3 k3 e^2); a bending triplet kb (e_i - 2 e_j + e_k)(...)^T in each
6// component; a hinge kh grad(theta) grad(theta)^T (exact at rest, where
7// theta = theta0). We return eigenpairs of M^{-1/2} K M^{-1/2} over a set of
8// degrees of freedom (dof: indices 3*bead + component).
9//
10// When every bead rests in a plane z = const, small z motion decouples from
11// in-plane motion (springs in the plane have no z-to-x/y coupling, triplets
12// separate by component, and flat hinges respond only to z), so the
13// "transverse" modes use z alone. Otherwise all three components are needed
14// ("full" modes), which is about 27 times the work of a transverse problem
15// with the same number of beads.
17import { hingeAngle } from "./hinge";
18import type { BeadSystem, Modes } from "./types";
20export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
22export function isPlanar(sys: BeadSystem): boolean {
23 const { n, pos0 } = sys;
24 let zmin = Infinity;
25 let zmax = -Infinity;
26 let span = 0;
27 for (let i = 0; i < n; i++) {
28 zmin = Math.min(zmin, pos0[3 * i + 2]);
29 zmax = Math.max(zmax, pos0[3 * i + 2]);
30 span = Math.max(span, Math.abs(pos0[3 * i]), Math.abs(pos0[3 * i + 1]));
31 }
32 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
33}
35export function freeBeads(sys: BeadSystem): number[] {
36 const out: number[] = [];
37 for (let i = 0; i < sys.n; i++) if (!sys.pinned[i]) out.push(i);
38 return out;
39}
41export function modeDofs(sys: BeadSystem, kind: Modes["kind"]): Int32Array {
42 const free = freeBeads(sys);
43 return kind === "transverse"
44 ? Int32Array.from(free, b => 3 * b + 2)
45 : Int32Array.from(free.flatMap(b => [3 * b, 3 * b + 1, 3 * b + 2]));
46}
48// M^{-1/2} K M^{-1/2} restricted to dof, dense, row-major.
49export function stiffnessMatrix(sys: BeadSystem, dof: Int32Array): Float64Array {
50 const { n, pos0, mass } = sys;
51 const nd = dof.length;
52 const index = new Int32Array(3 * n).fill(-1);
53 dof.forEach((q, r) => (index[q] = r));
54 const K = new Float64Array(nd * nd);
55 const add = (qa: number, qb: number, v: number) => {
56 const a = index[qa];
57 const b = index[qb];
58 if (a >= 0 && b >= 0) K[a * nd + b] += v;
59 };
60 for (let s = 0; s < sys.si.length; s++) {
61 const i = sys.si[s];
62 const j = sys.sj[s];
63 const dx = [0, 1, 2].map(c => pos0[3 * j + c] - pos0[3 * i + c]);
64 const L = Math.hypot(dx[0], dx[1], dx[2]);
65 const e = L - sys.sL0[s];
66 const keff = sys.sk[s] + 3 * sys.sk3[s] * e * e;
67 const g = (sys.sk[s] * e + sys.sk3[s] * e * e * e) / L;
68 const u = dx.map(d => d / L);
69 for (let a = 0; a < 3; a++) {
70 for (let b = 0; b < 3; b++) {
71 const B = (keff - g) * u[a] * u[b] + (a === b ? g : 0);
72 if (B === 0) continue;
73 add(3 * i + a, 3 * i + b, B);
74 add(3 * j + a, 3 * j + b, B);
75 add(3 * i + a, 3 * j + b, -B);
76 add(3 * j + a, 3 * i + b, -B);
77 }
78 }
79 }
80 const cb = [1, -2, 1];
81 for (let t = 0; t < sys.bi.length; t++) {
82 const beads = [sys.bi[t], sys.bj[t], sys.bk[t]];
83 for (let p = 0; p < 3; p++) {
84 for (let q = 0; q < 3; q++) {
85 for (let c = 0; c < 3; c++) add(3 * beads[p] + c, 3 * beads[q] + c, sys.bkb[t] * cb[p] * cb[q]);
86 }
87 }
88 }
89 const g = new Float64Array(12);
90 for (let h = 0; h < sys.hi.length; h++) {
91 const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
92 hingeAngle(pos0, beads[0], beads[1], beads[2], beads[3], g);
93 for (let p = 0; p < 12; p++) {
94 if (g[p] === 0) continue;
95 for (let q = 0; q < 12; q++) {
96 add(3 * beads[(p / 3) | 0] + (p % 3), 3 * beads[(q / 3) | 0] + (q % 3), sys.hkh[h] * g[p] * g[q]);
97 }
98 }
99 }
100 for (let a = 0; a < nd; a++) {
101 for (let b = 0; b < nd; b++) K[a * nd + b] /= Math.sqrt(mass[(dof[a] / 3) | 0] * mass[(dof[b] / 3) | 0]);
102 }
103 return K;
104}
106export function computeModes(sys: BeadSystem, eig: SymEig, kind: Modes["kind"] = isPlanar(sys) ? "transverse" : "full"): Modes {
107 const dof = modeDofs(sys, kind);
108 const nd = dof.length;
109 const { values, vectors } = eig(stiffnessMatrix(sys, dof), nd);
110 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
111 const freqs = new Float64Array(nd);
112 const psi = new Float64Array(nd * nd);
113 order.forEach((src, k) => {
114 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
115 psi.set(vectors.subarray(src * nd, (src + 1) * nd), k * nd);
116 });
117 return { freqs, psi, dof, kind };
118}
120// Small-amplitude decay rate of each mode (amplitude, 1/s): sigma_k =
121// (1/2) u^T C u with u the mode shape normalized so u^T M u = 1 and C the
122// damping matrix at rest (friction, dashpots, bending and hinge damping).
123// Exact when C is a combination of M and K (as bp_decay builds it for flat
124// nets); otherwise it neglects the coupling of modes by damping, which is
125// small for light damping.
126export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
127 const { n, mass, friction, di, dj, dc, pos0 } = sys;
128 const { dof, psi, freqs } = modes;
129 const nd = dof.length;
130 const u = new Float64Array(3 * n);
131 const out = new Float64Array(freqs.length);
132 const hgs = new Float64Array(12 * sys.hi.length);
133 for (let h = 0; h < sys.hi.length; h++) {
134 hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], hgs.subarray(12 * h, 12 * h + 12));
135 }
136 for (let k = 0; k < freqs.length; k++) {
137 u.fill(0);
138 for (let r = 0; r < nd; r++) u[dof[r]] = psi[k * nd + r] / Math.sqrt(mass[(dof[r] / 3) | 0]);
139 let c = 0;
140 for (let i = 0; i < n; i++) c += friction[i] * (u[3 * i] ** 2 + u[3 * i + 1] ** 2 + u[3 * i + 2] ** 2);
141 for (let d = 0; d < di.length; d++) {
142 let ax = 0;
143 let len = 0;
144 for (let q = 0; q < 3; q++) {
145 const r = u[3 * di[d] + q] - u[3 * dj[d] + q];
146 c += dc[d] * r * r;
147 const e = pos0[3 * dj[d] + q] - pos0[3 * di[d] + q];
148 ax += e * r;
149 len += e * e;
150 }
151 if (sys.dca[d] !== 0) c += (sys.dca[d] * ax * ax) / len;
152 }
153 for (let t = 0; t < sys.bi.length; t++) {
154 for (let q = 0; q < 3; q++) c += sys.bcb[t] * (u[3 * sys.bi[t] + q] - 2 * u[3 * sys.bj[t] + q] + u[3 * sys.bk[t] + q]) ** 2;
155 }
156 for (let h = 0; h < sys.hi.length; h++) {
157 if (sys.hch[h] === 0) continue;
158 const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
159 let s = 0;
160 for (let p = 0; p < 4; p++) for (let q = 0; q < 3; q++) s += hgs[12 * h + 3 * p + q] * u[3 * beads[p] + q];
161 c += sys.hch[h] * s * s;
162 }
163 out[k] = c / 2;
164 }
165 return out;
166}