1// Linear modes about the rest positions.
2//
3// For a flat net (all beads at one z, all springs in that plane), small
4// transverse (z) motion decouples from in-plane motion, and its stiffness is
5// a graph Laplacian weighted by spring tension over length,
6// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T; we use that n-by-n problem.
7// Otherwise we use the full 3n-by-3n stiffness of the springs. Either way we
8// return the eigenpairs of M^{-1/2} K M^{-1/2} over the free degrees of freedom.
10import { springBlocks } from "./forces";
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);
26}
28// Number of degrees of freedom computeModes would use.
29export function modeSize(sys: BeadSystem): number {
30 const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
31 return isPlanar(sys) ? nFree : 3 * nFree;
32}
34function stiffness(sys: BeadSystem): { K: Float64Array; dof: Int32Array } {
35 const { n, pinned, si, sj, mass } = sys;
36 const planar = isPlanar(sys);
37 const comps = planar ? [2] : [0, 1, 2];
38 const row = new Int32Array(3 * n).fill(-1);
39 const dof: number[] = [];
40 for (let i = 0; i < n; i++) if (!pinned[i]) for (const c of comps) row[3 * i + c] = dof.push(3 * i + c) - 1;
41 const nd = dof.length;
42 const K = new Float64Array(nd * nd);
43 const B = springBlocks(sys);
44 const idx = [
45 [0, 3, 4],
46 [3, 1, 5],
47 [4, 5, 2],
48 ];
49 const add = (r: number, c: number, v: number) => {
50 if (r >= 0 && c >= 0) K[r * nd + c] += v;
51 };
52 for (let s = 0; s < si.length; s++) {
53 const i3 = 3 * si[s];
54 const j3 = 3 * sj[s];
55 for (const a of comps) {
56 for (const b of comps) {
57 const v = B[6 * s + idx[a][b]];
58 add(row[i3 + a], row[i3 + b], v);
59 add(row[j3 + a], row[j3 + b], v);
60 add(row[i3 + a], row[j3 + b], -v);
61 add(row[j3 + a], row[i3 + b], -v);
62 }
63 }
64 }
65 for (let r = 0; r < nd; r++) {
66 for (let c = 0; c < nd; c++) K[r * nd + c] /= Math.sqrt(mass[Math.floor(dof[r] / 3)] * mass[Math.floor(dof[c] / 3)]);
67 }
68 return { K, dof: Int32Array.from(dof) };
69}
71export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
72 const { K, dof } = stiffness(sys);
73 const nd = dof.length;
74 const { values, vectors } = eig(K, nd);
75 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
76 const freqs = new Float64Array(nd);
77 const psi = new Float64Array(nd * nd);
78 order.forEach((src, k) => {
79 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
80 psi.set(vectors.subarray(src * nd, (src + 1) * nd), k * nd);
81 });
82 return { freqs, psi, dof, kind: isPlanar(sys) ? "transverse" : "full" };
83}
85// Small-amplitude decay rate of each mode (amplitude, 1/s), from the friction
86// and dashpots: sigma_k = (1/2) u^T C u with u the mode shape normalized so
87// u^T M u = 1. Exact when the damping matrix is a combination of M and K (as
88// bp_decay builds it for flat nets); otherwise it neglects the coupling of
89// modes by damping, which is small for light damping.
90export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
91 const { n, mass, friction, di, dj, dc } = sys;
92 const { dof, psi, freqs } = modes;
93 const nd = dof.length;
94 const u = new Float64Array(3 * n);
95 const out = new Float64Array(freqs.length);
96 for (let k = 0; k < freqs.length; k++) {
97 u.fill(0);
98 for (let r = 0; r < nd; r++) u[dof[r]] = psi[k * nd + r] / Math.sqrt(mass[Math.floor(dof[r] / 3)]);
99 let c = 0;
100 for (let i = 0; i < n; i++) c += friction[i] * (u[3 * i] ** 2 + u[3 * i + 1] ** 2 + u[3 * i + 2] ** 2);
101 for (let d = 0; d < di.length; d++) {
102 const i3 = 3 * di[d];
103 const j3 = 3 * dj[d];
104 c += dc[d] * ((u[i3] - u[j3]) ** 2 + (u[i3 + 1] - u[j3 + 1]) ** 2 + (u[i3 + 2] - u[j3 + 2]) ** 2);
105 }
106 out[k] = c / 2;
107 }
108 return out;
109}