f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1// CPU reference implementation of the spatial discretization (float64).
2// It mirrors step.wgsl exactly and is used for the stable time-step estimate
3// and for the tests. It is not used in the time loop of the browser simulation.
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 5import { type Grid, type Material, insideCell, insideNode, nodeWeight } from "../params.ts";
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 6
7export interface Fields {
8 u: Float64Array;
9 v: Float64Array;
10 w: Float64Array;
11}
13const G3 = 1 / Math.sqrt(3);
15/**
16 * Accelerations (au, av, aw) of every node for displacement fields (u, v, w),
17 * without damping. Bending forces come from the discrete Kirchhoff energy
18 * (natural free-edge conditions); membrane forces come from bilinear (Q4)
19 * elements with 2x2 Gauss integration and the von Karman strains.
20 */
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 21export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?: Float64Array): void {
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 22 const { nx, ny, hx, hy } = g;
23 const { u, v, w } = f;
24 const hx2 = hx * hx;
25 const hy2 = hy * hy;
26 const nu = m.nu;
27 const id = (i: number, j: number) => j * nx + i;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 29 // A second difference exists where the node and both of its neighbours belong to the sheet.
30 const ins = (i: number, j: number) => insideNode(g, i, j);
31 const hasA = (i: number, j: number) => ins(i - 1, j) && ins(i, j) && ins(i + 1, j);
32 const hasB = (i: number, j: number) => ins(i, j - 1) && ins(i, j) && ins(i, j + 1);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 33 // Bending moment-like quantities at nodes (zero where the second difference is undefined).
34 const Ma = (i: number, j: number): number => {
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 36 const a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
37 let b = 0;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 38 if (hasB(i, j)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
39 return nodeWeight(g, i, j) * (a + nu * b);
41 const Mb = (i: number, j: number): number => {
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 43 const b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
44 let a = 0;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 45 if (hasA(i, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
46 return nodeWeight(g, i, j) * (b + nu * a);
48 // Mixed second difference at a cell center (zero for cells outside the plate).
49 const Cc = (ci: number, cj: number): number => {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 50 if (!insideCell(g, ci, cj)) return 0;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 51 return (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
52 };
54 const q = (hx * hy) / 4; // area per Gauss point
55 const U = new Float64Array(4);
56 const V = new Float64Array(4);
57 const Wn = new Float64Array(4);
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 58 const W0 = new Float64Array(4);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 59 const dNdx = new Float64Array(4);
60 const dNdy = new Float64Array(4);
62 // Membrane force on local node `a` of cell (ci, cj); returns [Fu, Fv, Fw].
63 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 64 if (!insideCell(g, ci, cj)) return;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 65 const n0 = id(ci, cj), n1 = id(ci + 1, cj), n2 = id(ci + 1, cj + 1), n3 = id(ci, cj + 1);
66 U[0] = u[n0]; U[1] = u[n1]; U[2] = u[n2]; U[3] = u[n3];
67 V[0] = v[n0]; V[1] = v[n1]; V[2] = v[n2]; V[3] = v[n3];
68 Wn[0] = w[n0]; Wn[1] = w[n1]; Wn[2] = w[n2]; Wn[3] = w[n3];
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 69 if (w0) { W0[0] = w0[n0]; W0[1] = w0[n1]; W0[2] = w0[n2]; W0[3] = w0[n3]; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 70 for (let gp = 0; gp < 4; gp++) {
71 const xi = gp & 1 ? G3 : -G3;
72 const eta = gp & 2 ? G3 : -G3;
73 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
74 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
75 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
76 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
77 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
78 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
79 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
80 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 81 // Marguerre strains: total slope (w0 + w) minus the stress-free slope of w0.
82 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 83 for (let b = 0; b < 4; b++) {
84 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
85 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 86 wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
87 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 89 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
90 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
91 const gxy = uy + vx + wx * wy - w0x * w0y;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 92 const Nxx = m.C * (exx + nu * eyy);
93 const Nyy = m.C * (eyy + nu * exx);
94 const Nxy = m.G * gxy;
95 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
96 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
97 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
98 }
99 };
101 const F = new Float64Array(3);
102 for (let j = 0; j < ny; j++) {
103 for (let i = 0; i < nx; i++) {
104 const k = id(i, j);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 105 const weight = nodeWeight(g, i, j);
106 if (weight === 0) {
107 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
108 continue;
109 }
111 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
112 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
113 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
114 F[0] = 0; F[1] = 0; F[2] = -m.D * hx * hy * bend;
115 cellForce(i, j, 0, F);
116 cellForce(i - 1, j, 1, F);
117 cellForce(i - 1, j - 1, 2, F);
118 cellForce(i, j - 1, 3, F);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 119 const mass = m.rhoh * hx * hy * weight;
121 out.v[k] = F[1] / mass;
122 out.w[k] = F[2] / mass;
123 }
124 }
125}
127/**
128 * Largest angular frequency of the linearized discrete system, by power
129 * iteration on a small grid with the same spacing (the extreme mode is a
130 * grid-scale mode, so the plate size does not matter).
131 */
132export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 300): number {
133 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
134 const n = gs.n;
135 const x: Fields = { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
136 const y: Fields = { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
137 let seed = 12345;
138 const rnd = () => {
139 seed = (seed * 1664525 + 1013904223) >>> 0;
140 return seed / 4294967296 - 0.5;
141 };
142 for (let k = 0; k < n; k++) {
143 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
144 }
145 const eps = 1e-9; // keeps the nonlinear terms negligible
146 let lambda = 0;
147 for (let it = 0; it < iters; it++) {
148 let norm = 0;
149 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
150 norm = Math.sqrt(norm);
151 for (let k = 0; k < n; k++) {
152 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
153 }
154 accelerations(gs, m, x, y);
155 let dot = 0;
156 for (let k = 0; k < n; k++) dot += x.u[k] * y.u[k] + x.v[k] * y.v[k] + x.w[k] * y.w[k];
157 lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1}K
158 for (let k = 0; k < n; k++) {
159 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
160 }
161 }
162 return Math.sqrt(Math.max(lambda, 0));
163}
165/** Discrete potential energy (bending + membrane) matching `accelerations`. */
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 166export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Array): number {
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 167 const { nx, ny, hx, hy } = g;
168 const { u, v, w } = f;
169 const nu = m.nu;
170 const id = (i: number, j: number) => j * nx + i;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 171 const ins = (i: number, j: number) => insideNode(g, i, j);
173 for (let j = 0; j < ny; j++) {
174 for (let i = 0; i < nx; i++) {
175 let a = 0, b = 0;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 176 if (ins(i - 1, j) && ins(i, j) && ins(i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
177 if (ins(i, j - 1) && ins(i, j) && ins(i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
178 Ub += nodeWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
180 }
181 for (let cj = 0; cj < ny - 1; cj++) {
182 for (let ci = 0; ci < nx - 1; ci++) {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 183 if (!insideCell(g, ci, cj)) continue;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 184 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
185 Ub += 2 * (1 - nu) * c * c;
186 }
187 }
188 Ub *= 0.5 * m.D * hx * hy;
190 let Um = 0;
191 const q = (hx * hy) / 4;
192 for (let cj = 0; cj < ny - 1; cj++) {
193 for (let ci = 0; ci < nx - 1; ci++) {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 194 if (!insideCell(g, ci, cj)) continue;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 195 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
196 for (let gp = 0; gp < 4; gp++) {
197 const xi = gp & 1 ? G3 : -G3;
198 const eta = gp & 2 ? G3 : -G3;
199 const dNdx = [-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)].map((s) => (s / 4) * (2 / hx));
200 const dNdy = [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi].map((s) => (s / 4) * (2 / hy));
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 201 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 202 for (let b = 0; b < 4; b++) {
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 203 const z0 = w0 ? w0[nn[b]] : 0;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 204 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
205 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 206 wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
207 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 209 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
210 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
211 const gxy = uy + vx + wx * wy - w0x * w0y;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 212 Um += q * 0.5 * (m.C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + m.G * gxy * gxy);
213 }
214 }
215 }
216 return Ub + Um;
217}