a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 1// Spring forces and potential energy, shared by the time stepper and the
2// static relaxation.
4import type { BeadSystem } from "./types";
6export function springForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
7 const { si, sj, sk, sL0, sk3 } = sys;
8 const ns = si.length;
9 F.fill(0);
10 let pe = 0;
11 for (let s = 0; s < ns; s++) {
12 const i3 = 3 * si[s];
13 const j3 = 3 * sj[s];
14 const dx = x[j3] - x[i3];
15 const dy = x[j3 + 1] - x[i3 + 1];
16 const dz = x[j3 + 2] - x[i3 + 2];
17 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
18 const e = L - sL0[s];
19 const k = sk[s];
20 const k3 = sk3[s];
21 const e2 = e * e;
22 pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
23 const g = (k * e + k3 * e2 * e) / L;
24 F[i3] += g * dx;
25 F[i3 + 1] += g * dy;
26 F[i3 + 2] += g * dz;
27 F[j3] -= g * dx;
28 F[j3 + 1] -= g * dy;
29 F[j3 + 2] -= g * dz;
30 }
31 return pe;
32}
34// Constant load: adds it to F and returns its potential energy, -load . (x - pos0).
35export function loadForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
36 const load = sys.load;
37 if (!load) return 0;
38 let pe = 0;
39 for (let q = 0; q < load.length; q++) {
40 F[q] += load[q];
41 pe -= load[q] * (x[q] - sys.pos0[q]);
42 }
43 return pe;
44}
46// Per-spring 3x3 stiffness blocks at the rest positions,
47// B = keff d d^T + (f/L)(I - d d^T), stored as [xx, yy, zz, xy, xz, yz].
48// The stiffness matrix is the sum over springs of B on the (i,i) and (j,j)
49// blocks and -B on the (i,j) and (j,i) blocks.
50export function springBlocks(sys: BeadSystem): Float64Array {
51 const { si, sj, sk, sL0, sk3, pos0 } = sys;
52 const ns = si.length;
53 const B = new Float64Array(ns * 6);
54 for (let s = 0; s < ns; s++) {
55 const i3 = 3 * si[s];
56 const j3 = 3 * sj[s];
57 const dx = pos0[j3] - pos0[i3];
58 const dy = pos0[j3 + 1] - pos0[i3 + 1];
59 const dz = pos0[j3 + 2] - pos0[i3 + 2];
60 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
61 const e = L - sL0[s];
62 const keff = sk[s] + 3 * sk3[s] * e * e;
63 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
64 const u = [dx / L, dy / L, dz / L];
65 const a = keff - g;
66 B[6 * s] = a * u[0] * u[0] + g;
67 B[6 * s + 1] = a * u[1] * u[1] + g;
68 B[6 * s + 2] = a * u[2] * u[2] + g;
69 B[6 * s + 3] = a * u[0] * u[1];
70 B[6 * s + 4] = a * u[0] * u[2];
71 B[6 * s + 5] = a * u[1] * u[2];
72 }
73 return B;
74}