1// Static equilibrium of the springs under the constant load, by FIRE
2// (Bitzek et al. 2006): damped dynamics that keeps the velocity aligned with
3// the force and lengthens its step while the energy keeps falling.
5import { loadForces, springForces } from "./forces";
6import { estimateOmegaMax } from "./solver";
7import type { BeadSystem } from "./types";
9export function relax(sys: BeadSystem, maxIter = 200000): { pos: Float64Array; iterations: number; residual: number; maxTension: number } {
10 const { n, pinned, mass } = sys;
11 const N3 = 3 * n;
12 const x = new Float64Array(sys.pos0);
13 const v = new Float64Array(N3);
14 const F = new Float64Array(N3);
15 const invm = new Float64Array(n);
16 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
17 const force = () => {
18 springForces(sys, x, F);
19 loadForces(sys, x, F);
20 for (let i = 0; i < n; i++) if (pinned[i]) F.fill(0, 3 * i, 3 * i + 3);
21 };
22 const residual = () => {
23 let r = 0;
24 for (let i = 0; i < n; i++) r = Math.max(r, Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]));
25 return r;
26 };
27 let scale = 0; // force scale for the tolerance: largest spring or load force
28 for (let s = 0; s < sys.si.length; s++) scale = Math.max(scale, sys.sk[s] * Math.abs(sys.sL0[s]) * 1e-3);
29 if (sys.load) for (const l of sys.load) scale = Math.max(scale, Math.abs(l));
30 const tol = 1e-10 * Math.max(scale, 1e-30);
32 const dt0 = 0.2 / estimateOmegaMax(sys);
33 let dt = dt0;
34 const dtMax = 10 * dt0;
35 let alpha = 0.1;
36 let positive = 0;
37 force();
38 let it = 0;
39 for (; it < maxIter; it++) {
40 if (residual() < tol) break;
41 let P = 0, vn = 0, fn = 0;
42 for (let q = 0; q < N3; q++) {
43 P += F[q] * v[q];
44 vn += v[q] * v[q];
45 fn += F[q] * F[q];
46 }
47 vn = Math.sqrt(vn);
48 fn = Math.sqrt(fn);
49 if (P > 0) {
50 for (let q = 0; q < N3; q++) v[q] = (1 - alpha) * v[q] + (alpha * vn * F[q]) / (fn || 1);
51 if (++positive > 5) {
52 dt = Math.min(dt * 1.1, dtMax);
53 alpha *= 0.99;
54 }
55 } else {
56 v.fill(0);
57 dt *= 0.5;
58 alpha = 0.1;
59 positive = 0;
60 }
61 for (let i = 0; i < n; i++) {
62 const h = 0.5 * dt * invm[i];
63 v[3 * i] += h * F[3 * i];
64 v[3 * i + 1] += h * F[3 * i + 1];
65 v[3 * i + 2] += h * F[3 * i + 2];
66 }
67 for (let q = 0; q < N3; q++) x[q] += dt * v[q];
68 force();
69 for (let i = 0; i < n; i++) {
70 const h = 0.5 * dt * invm[i];
71 v[3 * i] += h * F[3 * i];
72 v[3 * i + 1] += h * F[3 * i + 1];
73 v[3 * i + 2] += h * F[3 * i + 2];
74 }
75 }
76 let maxTension = 0;
77 for (let s = 0; s < sys.si.length; s++) {
78 const i3 = 3 * sys.si[s], j3 = 3 * sys.sj[s];
79 const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
80 maxTension = Math.max(maxTension, Math.abs(sys.sk[s] * (L - sys.sL0[s])));
81 }
82 return { pos: x, iterations: it, residual: residual(), maxTension };
83}