1// Pieces shared by the CPU and GPU solvers: the static forces and energy in
2// double precision, the checks of the rest state, the choice of time step,
3// and the projection of the motion onto linear modes.
5import { angleDiff, hingeAngle } from "./hinge";
6import { estimateOmegaMax, stabilityLimit } from "./solver";
7import type { BeadSystem, Modes } from "./types";
9// Spring, bending and hinge forces at positions x (added into F, which is
10// cleared first); returns the potential energy.
11export function staticForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
12 F.fill(0);
13 let pe = 0;
14 for (let s = 0; s < sys.si.length; s++) {
15 const i3 = 3 * sys.si[s];
16 const j3 = 3 * sys.sj[s];
17 const dx = x[j3] - x[i3];
18 const dy = x[j3 + 1] - x[i3 + 1];
19 const dz = x[j3 + 2] - x[i3 + 2];
20 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
21 const e = L - sys.sL0[s];
22 const e2 = e * e;
23 pe += 0.5 * sys.sk[s] * e2 + 0.25 * sys.sk3[s] * e2 * e2;
24 const g = (sys.sk[s] * e + sys.sk3[s] * e2 * e) / L;
25 F[i3] += g * dx;
26 F[i3 + 1] += g * dy;
27 F[i3 + 2] += g * dz;
28 F[j3] -= g * dx;
29 F[j3 + 1] -= g * dy;
30 F[j3 + 2] -= g * dz;
31 }
32 for (let b = 0; b < sys.bi.length; b++) {
33 const i3 = 3 * sys.bi[b];
34 const j3 = 3 * sys.bj[b];
35 const k3 = 3 * sys.bk[b];
36 for (let c = 0; c < 3; c++) {
37 const d = x[i3 + c] - 2 * x[j3 + c] + x[k3 + c];
38 pe += 0.5 * sys.bkb[b] * d * d;
39 F[i3 + c] -= sys.bkb[b] * d;
40 F[j3 + c] += 2 * sys.bkb[b] * d;
41 F[k3 + c] -= sys.bkb[b] * d;
42 }
43 }
44 const g = new Float64Array(12);
45 for (let h = 0; h < sys.hi.length; h++) {
46 const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
47 const d = angleDiff(hingeAngle(x, beads[0], beads[1], beads[2], beads[3], g), sys.hth0[h]);
48 pe += 0.5 * sys.hkh[h] * d * d;
49 for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) F[3 * beads[p] + c] -= sys.hkh[h] * d * g[3 * p + c];
50 }
51 return pe;
52}
54export interface Prepared {
55 peRest: number;
56 residualForce: number;
57 maxTension: number;
58 omegaMax: number;
59 substeps: number;
60 h: number;
61 nSamples: number;
62 warnings: string[];
63}
65export function prepare(sys: BeadSystem): Prepared {
66 const { n, pinned, pos0 } = sys;
67 const warnings: string[] = [];
68 const F = new Float64Array(3 * n);
69 const peRest = staticForces(sys, pos0, F);
70 let residualForce = 0;
71 for (let i = 0; i < n; i++) {
72 if (!pinned[i]) residualForce = Math.max(residualForce, Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]));
73 }
74 let maxTension = 0;
75 let forceScale = 0; // the larger of the rest tension and the force of a 1e-3 strain on any element
76 for (let s = 0; s < sys.si.length; s++) {
77 const i3 = 3 * sys.si[s];
78 const j3 = 3 * sys.sj[s];
79 const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
80 const e = L - sys.sL0[s];
81 maxTension = Math.max(maxTension, Math.abs(sys.sk[s] * e + sys.sk3[s] * e * e * e));
82 forceScale = Math.max(forceScale, 1e-3 * sys.sk[s] * sys.sL0[s]);
83 }
84 forceScale = Math.max(forceScale, maxTension);
85 for (let b = 0; b < sys.bi.length; b++) {
86 const i3 = 3 * sys.bi[b];
87 const j3 = 3 * sys.bj[b];
88 const len = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
89 forceScale = Math.max(forceScale, 1e-3 * sys.bkb[b] * len);
90 }
91 const g = new Float64Array(12);
92 for (let h = 0; h < sys.hi.length; h++) {
93 hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], g);
94 forceScale = Math.max(forceScale, 1e-3 * sys.hkh[h] * Math.max(...g.map(Math.abs)));
95 }
96 if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
97 warnings.push(
98 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
99 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
100 );
101 }
103 const omegaMax = estimateOmegaMax(sys);
104 const fs = sys.sampleRate;
105 const hLimit = stabilityLimit(sys.method) / omegaMax;
106 let substeps: number;
107 if (sys.dt > 0) {
108 substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
109 if (1 / (fs * substeps) > hLimit) {
110 warnings.push(`dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`);
111 }
112 } else {
113 substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
114 }
115 return {
116 peRest,
117 residualForce,
118 maxTension,
119 omegaMax,
120 substeps,
121 h: 1 / (fs * substeps),
122 nSamples: Math.max(1, Math.round(sys.duration * fs)),
123 warnings,
124 };
125}
127// Energy of each linear mode, E_k = (qdot_k^2 + omega_k^2 q_k^2)/2, from the
128// displacements and velocities.
129export class ModeTracker {
130 private bz: Float64Array;
131 private bv: Float64Array;
132 private om2: Float64Array;
133 constructor(private modes: Modes, private mass: Float64Array, private pos0: Float64Array) {
134 this.bz = new Float64Array(modes.dof.length);
135 this.bv = new Float64Array(modes.dof.length);
136 this.om2 = Float64Array.from(modes.freqs, f => (2 * Math.PI * f) ** 2);
137 }
138 // x: positions (or displacements if relative is true); writes nModes energies into out at offset
139 record(x: ArrayLike<number>, v: ArrayLike<number>, out: Float32Array, offset: number, relative = false) {
140 const { dof, psi } = this.modes;
141 const nf = dof.length;
142 for (let r = 0; r < nf; r++) {
143 const q = dof[r];
144 const sm = Math.sqrt(this.mass[(q / 3) | 0]);
145 this.bz[r] = sm * (relative ? x[q] : x[q] - this.pos0[q]);
146 this.bv[r] = sm * v[q];
147 }
148 const nModes = this.om2.length;
149 for (let k = 0; k < nModes; k++) {
150 let q = 0;
151 let qd = 0;
152 const off = k * nf;
153 for (let r = 0; r < nf; r++) {
154 q += psi[off + r] * this.bz[r];
155 qd += psi[off + r] * this.bv[r];
156 }
157 out[offset + k] = 0.5 * (qd * qd + this.om2[k] * q * q);
158 }
159 }
160}