1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// The discrete shell, as in cymsim/shell.py: thickness and height sampled on the mesh, the CST
2// slope of the nodal heights, element stiffnesses, lumped mass with the inflated rotational
3// inertia, the six M-orthonormal rigid-body modes, the node-to-element adjacency used by the
4// gather of element forces, and the Irons element bound on the largest eigenfrequency that sets
5// the stable time step. Everything here is float64 and computed once per configuration.
6import type { Material, Numerics } from "../config/schema.ts";
7import { dktStiffness } from "./dkt.ts";
8import type { CymbalProfile } from "./geometry.ts";
9import { symmetricEigenvalues } from "./linalg.ts";
10import type { Mesh } from "./mesh.ts";
12export const N_DOF = 5; // (u, v, w, theta_x, theta_y) per node
14export interface Shell {
15 mesh: Mesh; material: Material; numerics: Numerics; nl: number;
16 hE: Float64Array; hNode: Float64Array; z0Node: Float64Array; z0x: Float64Array; z0y: Float64Array;
13f4dbdexport .cymsim fileJeremy Magland 17 normalNode: Float64Array; // (N, 3) unit normal of the rest surface at the nodes, +z side (front)
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 18 De: Float64Array; Cm: Float64Array; // (E,) and (E, 3) = (C11, C12, C33)
19 Ke: Float64Array; // (E, 81)
20 rhoHA: Float64Array; // (E,) rho h A
21 M: Float64Array; Minv: Float64Array; // (N, 5)
22 Rm: Float64Array; // (6, N, 5), M-orthonormal
76fcbe0GPU: rigid modes evaluated from node positions; convert older states on importJeremy Magland 23 rigidCentre: [number, number]; // mean node position, the origin of rawRigidModes
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 24 adjOffsets: Uint32Array; adjEntries: Uint32Array; // CSR: node -> list of (3 e + local index)
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 25 mass: number; maxSlope: number; limitingElement: number;
26 omegaMax: number; // element (Irons) bound on sqrt(lambda_max(M^-1 K))
27 omegaGlobal: number; omegaGlobalIterations: number; omegaGlobalConverged: boolean; // by power iteration
28 omegaStep: number; dtStable: number; // the bound the time step uses, per numerics.stability_bound, and 2 / omegaStep
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 29 zeroModes: [number, number]; // min and max count of zero-energy element modes (expect 6)
30}
76fcbe0GPU: rigid modes evaluated from node positions; convert older states on importJeremy Magland 32/** The six rigid-body modes of the curved shell at a node at (X, Y) from the rigid centre with rest
33 * height z, before M-orthonormalisation: (u, v, w, theta_x, theta_y) of the two in-plane translations,
34 * the rotation about z, the normal translation and the two tilts. */
35export function rawRigidModes(X: number, Y: number, z: number): number[][] {
36 return [[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [-Y, X, 0, 0, 0], [0, 0, 1, 0, 0], [0, -z, Y, 1, 0], [z, 0, -X, 0, 1]];
37}
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 39/** sqrt(lambda_max(M^-1 K)) of the stiffness assembled from the element matrices Kt (E x 15 x 15,
40 * node-major DOFs), by power iteration on M^-1/2 K M^-1/2 with the lumped M (N x 5). The start
41 * vector x_i = sin(i + 1) and the stopping rule (the Rayleigh quotient changes by at most tol
42 * relative over 10 iterations) are those of cymsim's shell.global_omega_max, so that both packages
43 * arrive at the same value. The Rayleigh quotient approaches lambda_max from below. */
44export function globalOmegaMax(Kt: Float64Array, M: Float64Array, tri: Uint32Array, tol = 1e-7, maxIter = 20000):
45 { omega: number; iterations: number; converged: boolean } {
46 const E = tri.length / 3, L = M.length;
47 const dof = new Uint32Array(15 * E);
48 for (let e = 0; e < E; e++) for (let i = 0; i < 3; i++) for (let d = 0; d < N_DOF; d++) dof[15 * e + N_DOF * i + d] = N_DOF * tri[3 * e + i] + d;
49 const isq = M.map((m) => 1 / Math.sqrt(m));
50 const x = new Float64Array(L), y = new Float64Array(L);
51 const xe = new Float64Array(15), hist: number[] = [];
52 let nrm = 0;
53 for (let i = 0; i < L; i++) { x[i] = Math.sin(i + 1); nrm += x[i] * x[i]; }
54 nrm = Math.sqrt(nrm);
55 for (let i = 0; i < L; i++) x[i] /= nrm;
56 for (let it = 1; it <= maxIter; it++) {
57 y.fill(0);
58 for (let e = 0; e < E; e++) {
59 for (let a = 0; a < 15; a++) { const g = dof[15 * e + a]; xe[a] = x[g] * isq[g]; }
60 for (let a = 0; a < 15; a++) {
61 let s = 0;
62 for (let b = 0; b < 15; b++) s += Kt[225 * e + 15 * a + b] * xe[b];
63 y[dof[15 * e + a]] += s;
64 }
65 }
66 let lam = 0, yy = 0;
67 for (let i = 0; i < L; i++) { y[i] *= isq[i]; lam += x[i] * y[i]; yy += y[i] * y[i]; }
68 hist.push(lam);
69 const n = Math.sqrt(yy);
70 for (let i = 0; i < L; i++) x[i] = y[i] / n;
71 if (it >= 20 && it % 10 === 0 && Math.abs(lam - hist[it - 11]) <= tol * lam) return { omega: Math.sqrt(lam), iterations: it, converged: true };
72 }
73 return { omega: Math.sqrt(hist[maxIter - 1]), iterations: maxIter, converged: false };
74}
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 76export function buildShell(mesh: Mesh, profile: CymbalProfile, material: Material, numerics: Numerics): Shell {
77 const N = mesh.n, E = mesh.e, { xy, tri, area, Gx, Gy } = mesh;
78 const nl = numerics.nonlinear ? 1 : 0;
79 const hE = new Float64Array(E), hNode = new Float64Array(N), z0Node = new Float64Array(N);
80 for (let e = 0; e < E; e++) hE[e] = profile.thickness(mesh.centroid[2 * e], mesh.centroid[2 * e + 1]);
81 for (let n = 0; n < N; n++) { hNode[n] = profile.thickness(xy[2 * n], xy[2 * n + 1]); z0Node[n] = profile.z0(xy[2 * n], xy[2 * n + 1]); }
82 const z0x = new Float64Array(E), z0y = new Float64Array(E);
83 let maxSlope = 0;
84 for (let e = 0; e < E; e++) {
85 let sx = 0, sy = 0;
86 for (let i = 0; i < 3; i++) { const z = z0Node[tri[3 * e + i]]; sx += Gx[3 * e + i] * z; sy += Gy[3 * e + i] * z; }
87 z0x[e] = sx; z0y[e] = sy; maxSlope = Math.max(maxSlope, Math.hypot(sx, sy));
88 }
13f4dbdexport .cymsim fileJeremy Magland 89 // nodal normal from the area-weighted mean of the adjacent elements' slopes, as in cymsim/shell.py
90 const normalNode = new Float64Array(3 * N), gxN = new Float64Array(N), gyN = new Float64Array(N), aN = new Float64Array(N);
91 for (let e = 0; e < E; e++) for (let i = 0; i < 3; i++) { const n = tri[3 * e + i]; gxN[n] += area[e] * z0x[e]; gyN[n] += area[e] * z0y[e]; aN[n] += area[e]; }
92 for (let n = 0; n < N; n++) {
93 const gx = gxN[n] / aN[n], gy = gyN[n] / aN[n], nrm = Math.sqrt(1 + gx * gx + gy * gy);
94 normalNode[3 * n] = -gx / nrm; normalNode[3 * n + 1] = -gy / nrm; normalNode[3 * n + 2] = 1 / nrm;
95 }
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 96 const { E: Ey, nu, rho } = material;
97 const De = new Float64Array(E), Cm = new Float64Array(3 * E), Ke = new Float64Array(81 * E), rhoHA = new Float64Array(E);
98 for (let e = 0; e < E; e++) {
99 De[e] = (Ey * hE[e] ** 3) / (12 * (1 - nu * nu));
100 const c = (Ey * hE[e]) / (1 - nu * nu);
101 Cm[3 * e] = c; Cm[3 * e + 1] = c * nu; Cm[3 * e + 2] = c * 0.5 * (1 - nu);
102 rhoHA[e] = rho * hE[e] * area[e];
103 const xe = [xy[2 * tri[3 * e]], xy[2 * tri[3 * e + 1]], xy[2 * tri[3 * e + 2]]];
104 const ye = [xy[2 * tri[3 * e] + 1], xy[2 * tri[3 * e + 1] + 1], xy[2 * tri[3 * e + 2] + 1]];
105 Ke.set(dktStiffness(xe, ye, De[e], nu), 81 * e);
106 }
107 // lumped mass
108 const M = new Float64Array(N_DOF * N);
109 const mE = new Float64Array(E), IE = new Float64Array(E);
110 for (let e = 0; e < E; e++) {
111 mE[e] = (rho * hE[e] * area[e]) / 3;
112 IE[e] = mE[e] * ((hE[e] * hE[e]) / 12 + numerics.rot_inertia_factor * area[e]);
113 for (let k = 0; k < 3; k++) {
114 const n = tri[3 * e + k];
115 M[5 * n] += mE[e] * numerics.inplane_mass_scale; M[5 * n + 1] += mE[e] * numerics.inplane_mass_scale;
116 M[5 * n + 2] += mE[e]; M[5 * n + 3] += IE[e]; M[5 * n + 4] += IE[e];
117 }
118 }
119 const Minv = M.map((m) => 1 / m);
120 let mass = 0;
121 for (let n = 0; n < N; n++) mass += M[5 * n + 2];
122 // rigid-body modes of the curved shell, M-orthonormalised by Gram-Schmidt
123 let cx = 0, cy = 0;
124 for (let n = 0; n < N; n++) { cx += xy[2 * n]; cy += xy[2 * n + 1]; }
125 cx /= N; cy /= N;
126 const Rm = new Float64Array(6 * N_DOF * N);
127 const set = (k: number, n: number, d: number, val: number) => { Rm[(k * N + n) * N_DOF + d] = val; };
128 for (let n = 0; n < N; n++) {
76fcbe0GPU: rigid modes evaluated from node positions; convert older states on importJeremy Magland 129 const r = rawRigidModes(xy[2 * n] - cx, xy[2 * n + 1] - cy, z0Node[n]);
130 for (let k = 0; k < 6; k++) for (let d = 0; d < N_DOF; d++) set(k, n, d, r[k][d]);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 131 }
132 const L = N_DOF * N;
133 const mdot = (a: number, b: number) => { let s = 0; for (let i = 0; i < L; i++) s += Rm[a * L + i] * M[i] * Rm[b * L + i]; return s; };
134 for (let k = 0; k < 6; k++) {
135 for (let l = 0; l < k; l++) { const c = mdot(l, k); for (let i = 0; i < L; i++) Rm[k * L + i] -= c * Rm[l * L + i]; }
136 const nrm = Math.sqrt(mdot(k, k));
137 for (let i = 0; i < L; i++) Rm[k * L + i] /= nrm;
138 }
139 // adjacency (CSR)
140 const deg = new Uint32Array(N);
141 for (let e = 0; e < E; e++) for (let k = 0; k < 3; k++) deg[tri[3 * e + k]]++;
142 const adjOffsets = new Uint32Array(N + 1);
143 for (let n = 0; n < N; n++) adjOffsets[n + 1] = adjOffsets[n] + deg[n];
144 const adjEntries = new Uint32Array(adjOffsets[N]), fill = new Uint32Array(N);
145 for (let e = 0; e < E; e++) for (let k = 0; k < 3; k++) { const n = tri[3 * e + k]; adjEntries[adjOffsets[n] + fill[n]++] = 3 * e + k; }
146 // Irons bound: lambda_max(M^-1 K) <= max_e lambda_max(M_e^-1 K_e), K_e the element tangent at rest
147 let lamMax = 0, limitingElement = 0, zmin = 15, zmax = 0;
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 148 const S = new Float64Array(225), me = new Float64Array(15), Kt = new Float64Array(225 * E);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 149 for (let e = 0; e < E; e++) {
150 elementTangent(e, Gx, Gy, area[e], z0x[e], z0y[e], Cm, Ke, S);
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 151 Kt.set(S, 225 * e);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 152 for (let i = 0; i < 3; i++) {
153 me[5 * i] = mE[e] * numerics.inplane_mass_scale; me[5 * i + 1] = mE[e] * numerics.inplane_mass_scale;
154 me[5 * i + 2] = mE[e]; me[5 * i + 3] = IE[e]; me[5 * i + 4] = IE[e];
155 }
156 for (let i = 0; i < 15; i++) for (let j = 0; j < 15; j++) S[15 * i + j] /= Math.sqrt(me[i] * me[j]);
157 const ev = symmetricEigenvalues(S, 15);
158 const top = ev[14];
159 let nz = 0; for (const v of ev) if (v < 1e-9 * top) nz++;
160 zmin = Math.min(zmin, nz); zmax = Math.max(zmax, nz);
161 if (top > lamMax) { lamMax = top; limitingElement = e; }
162 }
163 const omegaMax = Math.sqrt(lamMax);
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 164 // the global lambda_max(M^-1 K), which the element bound overestimates (about 1.45x at 6 mm elements)
165 const glob = globalOmegaMax(Kt, M, tri);
166 const useGlobal = numerics.stability_bound === "global" && glob.converged;
167 const omegaStep = useGlobal ? Math.min(glob.omega, omegaMax) : omegaMax;
76fcbe0GPU: rigid modes evaluated from node positions; convert older states on importJeremy Magland 168 return { mesh, material, numerics, nl, hE, hNode, z0Node, z0x, z0y, normalNode, De, Cm, Ke, rhoHA, M, Minv, Rm, rigidCentre: [cx, cy], adjOffsets, adjEntries,
c24896fTime step from the global stability bound, as in cymsim; Resume falls back to the element bound for older runs; accept output.probesJeremy Magland 169 mass, maxSlope, omegaMax, omegaGlobal: glob.omega, omegaGlobalIterations: glob.iterations, omegaGlobalConverged: glob.converged,
170 omegaStep, dtStable: 2 / omegaStep, limitingElement, zeroModes: [zmin, zmax] };
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 171}
173/** 15 x 15 tangent stiffness of one element at rest: A B^T C B (membrane, linear part) plus DKT on the bending slots. */
174export function elementTangent(e: number, Gx: Float64Array, Gy: Float64Array, A: number, z0x: number, z0y: number,
175 Cm: Float64Array, Ke: Float64Array, out: Float64Array) {
176 out.fill(0);
177 // B rows: exx, eyy, gxy over the 15 DOFs (node-major u, v, w, thx, thy)
178 const B = [new Float64Array(15), new Float64Array(15), new Float64Array(15)];
179 for (let i = 0; i < 3; i++) {
180 const gx = Gx[3 * e + i], gy = Gy[3 * e + i];
181 B[0][5 * i] = gx; B[0][5 * i + 2] = z0x * gx;
182 B[1][5 * i + 1] = gy; B[1][5 * i + 2] = z0y * gy;
183 B[2][5 * i] = gy; B[2][5 * i + 1] = gx; B[2][5 * i + 2] = z0x * gy + z0y * gx;
184 }
185 const C = [[Cm[3 * e], Cm[3 * e + 1], 0], [Cm[3 * e + 1], Cm[3 * e], 0], [0, 0, Cm[3 * e + 2]]];
186 for (let i = 0; i < 15; i++) for (let j = 0; j < 15; j++) {
187 let s = 0;
188 for (let k = 0; k < 3; k++) for (let l = 0; l < 3; l++) s += B[k][i] * C[k][l] * B[l][j];
189 out[15 * i + j] = A * s;
190 }
191 const slot = [2, 3, 4, 7, 8, 9, 12, 13, 14]; // (w, thx, thy) x 3 within the 15
192 for (let i = 0; i < 9; i++) for (let j = 0; j < 9; j++) out[15 * slot[i] + slot[j]] += Ke[81 * e + 9 * i + j];
193}