GPU: rigid modes evaluated from node positions; convert older states on import - proof-of-concept/cymsimweb

proof-of-concept / cymsimweb

GPU: rigid modes evaluated from node positions; convert older states on import

Jeremy Magland <jmagland@flatironinstitute.org> committed 5 days ago commit 76fcbe07caa9 parent 2b313b8 Browse files

5 changed files+89−21

THEORY.mdmodified+2−2View file

@@ -34,7 +34,7 @@ A sub-step needs three global synchronizations in the Python arrangement: after

3434

3535 **Element pass** (one thread per element, plus one thread per hit): reads $q$ and $v$, writes the fifteen forces of the element; accumulates the element's radiated acceleration (section 4); the hit threads evaluate the contact forces of their patch points from $q$ and the exciter's face height, applying the release at the hit's step.

3636

37-**Node pass** (one thread per node, workgroups of 256): sums the element forces of the node and the contact forces of the patch points it belongs to; removes from $v$ and from $q$ the rigid content that the previous step's update introduced, using the twelve coefficients $c_k = r_k^{\mathsf T} M v$ and $d_k = r_k^{\mathsf T} M q$ that the previous node pass summed; performs the leapfrog update $v' = (v + \Delta t M^{-1} f)/(1 + 2 \sigma_0 \Delta t)$, $q' = q + \Delta t v'$; and contributes $r_k^{\mathsf T} M v'$ and $r_k^{\mathsf T} M q'$ to a per-workgroup partial sum (a shared-memory tree reduction) for the next step. The first thread also advances the exciters, $v_s' = v_s + \Delta t F / m$, $z_s' = z_s + \Delta t v_s'$.

37+**Node pass** (one thread per node, workgroups of 256): sums the element forces of the node and the contact forces of the patch points it belongs to; removes from $v$ and from $q$ the rigid content that the previous step's update introduced, $v \leftarrow v - R G^{-1} R^{\mathsf T} M v$ and the same for $q$, using the twelve sums $R^{\mathsf T} M v$ and $R^{\mathsf T} M q$ that the previous node pass formed; performs the leapfrog update $v' = (v + \Delta t M^{-1} f)/(1 + 2 \sigma_0 \Delta t)$, $q' = q + \Delta t v'$; and contributes its share of $R^{\mathsf T} M v'$ and $R^{\mathsf T} M q'$ to a per-workgroup partial sum (a shared-memory tree reduction) for the next step. Here $R$ holds the six rigid modes before orthonormalisation (`rawRigidModes` in `src/model/shell.ts`: the two in-plane translations, the rotation about $z$, the normal translation and the two tilts of the curved shell), which the thread evaluates from the node's position and rest height, and $G = R^{\mathsf T} M R$ is their $6 \times 6$ Gram matrix, whose inverse is a constant of the shader. $R G^{-1} R^{\mathsf T} M$ is the same projection as $\sum_k r_k r_k^{\mathsf T} M$ with the M-orthonormal modes $r_k$ that the CPU loop uses, but it reads three floats per node instead of thirty. States exported before this arrangement hold the sums for the $r_k$; they are converted on import with the triangular matrix $R^{\mathsf T} M r$, which continues such a state exactly. The first thread also advances the exciters, $v_s' = v_s + \Delta t F / m$, $z_s' = z_s + \Delta t v_s'$.

3838

3939 Because the projection is applied one step late, the state carries at any time the rigid content of exactly one step's increment: a rigid velocity of order $F \Delta t / M_{shell}$ and a rigid displacement of order $\Delta t$ times that, both removed at the next step. Projecting $q$ as well as $v$ matters: without it the rigid displacement accumulates over a contact (about $10^{-5}$ m over a mallet stroke at coarse resolution, which changed the peak contact force by $7 \times 10^{-4}$); with it the linear reference run agrees with Python to $10^{-5}$ in peak force and rebound velocity and $10^{-4}$ in the energies.

4040

@@ -52,7 +52,7 @@ The acceleration that radiates is the difference of consecutive *projected* velo

5252

5353 ## 5. Buffers and precision

5454

55-The GPU keeps the set-up and the state in storage buffers, at most eight in any one pipeline (the default WebGPU limit): per-element data (CST operators, area, slopes, membrane stiffness, damping factor and the 45 entries of the upper triangle of $K_e$, 58 floats per element), element node indices, the state (q, v and the previous projected $(u, v, w)$ velocities, 13 floats per node), the element forces, per-node data (inverse mass, mass, the six rigid modes and the nodal normal, 43 floats per node, with the decimation taps at the tail), the adjacency in CSR form, the hits (static patch data and the exciter state, one float and one integer buffer), the synchronization buffer (step counters and partial sums), the radiation buffer (the $W N$ accumulator slots, the node history, the current sample's $M N$ sub-step accelerations, and the chunk's acceleration output rows when recording), the element history, the radiation coefficients (delays and interpolator taps, 36 floats per element) and the chunk's pressure output. The per-element and per-node data are stored field-major (field $f$ of element $e$ at $f E + e$), so that the threads of a workgroup read consecutive words; $K_e$, symmetric to round-off, is symmetrised in float64 before packing, and the product $K_e U_e$ is written out term by term in the shader. Together these halved the time of the element pass on an integrated GPU. A uniform holds the chunk's sample and frame bases. Everything on the GPU is float32. The set-up is float64 and is packed once. Diagnostics (energies, exciter state) are read back between chunks; the plate view reads $w$ from the state buffer directly in its vertex shader without a readback. The whole integrator state (state, sync, hits, radiation and element-history buffers, plus the step count) can be read back and written again, which is how a partial `.cymsim` file is continued.

55+The GPU keeps the set-up and the state in storage buffers, at most eight in any one pipeline (the default WebGPU limit): per-element data (CST operators, area, slopes, membrane stiffness, damping factor and the 45 entries of the upper triangle of $K_e$, 58 floats per element), element node indices, the state (q, v and the previous projected $(u, v, w)$ velocities, 13 floats per node), the element forces, per-node data (inverse mass, mass, the position relative to the rigid centre, the rest height and the nodal normal, 16 floats per node, with the decimation taps at the tail), the adjacency in CSR form, the hits (static patch data and the exciter state, one float and one integer buffer), the synchronization buffer (step counters and partial sums), the radiation buffer (the $W N$ accumulator slots, the node history, the current sample's $M N$ sub-step accelerations, and the chunk's acceleration output rows when recording), the element history, the radiation coefficients (delays and interpolator taps, 36 floats per element) and the chunk's pressure output. The per-element and per-node data are stored field-major (field $f$ of element $e$ at $f E + e$), so that the threads of a workgroup read consecutive words; $K_e$, symmetric to round-off, is symmetrised in float64 before packing, and the product $K_e U_e$ is written out term by term in the shader. Together these halved the time of the element pass on an integrated GPU. A uniform holds the chunk's sample and frame bases. Everything on the GPU is float32. The set-up is float64 and is packed once. Diagnostics (energies, exciter state) are read back between chunks; the plate view reads $w$ from the state buffer directly in its vertex shader without a readback. The whole integrator state (state, sync, hits, radiation and element-history buffers, plus the step count) can be read back and written again, which is how a partial `.cymsim` file is continued.

5656

5757 Float32 is adequate for the leapfrog at these step sizes: the displacement increment per step is about $10^{-3}$ of the displacement, so the update loses a few bits per step in $q$ without a visible effect on the energy balance, and the run is chaotic anyway once max|w|/h exceeds a few units, so no float32 trajectory tracks a float64 one for long. What is checked is the deterministic part (section 6).

5858

src/gpu/shaders.tsmodified+30−7View file

@@ -10,15 +10,19 @@

1010 export interface ShaderConstants {

1111 N: number; E: number; K: number; P: number; M: number; W: number; D: number; G: number; T: number;

1212 accRows: number; dt: number; sigma0: number; nl: number;

13+ /** Inverse of the M-Gram matrix of the six raw rigid modes (rawRigidModes in shell.ts), row-major. */

14+ rigidGramInv: ArrayLike<number>;

1315 }

1416

1517 // Per-element and per-node data are stored field-major (field f of element e at f * E + e), so that

1618 // neighbouring threads read neighbouring words.

1719 export const ELEM_STRIDE = 58; // fields: Gx[3] Gy[3] A z0x z0y C11 C12 C33 cd Ke[45] (upper triangle, row by row)

18-export const NODE_STRIDE = 43; // fields: Minv[5] M[5] Rm[6][5] normal[3]; the decimation taps follow the nodes

20+export const NODE_STRIDE = 16; // fields: Minv[5] M[5] X Y z0 normal[3] (X, Y from the rigid centre); the decimation taps follow the nodes

1921 export const HIT_HEADER = 12; // kPoint alpha mass gap velocity zs vs F Fmax pad[3]

2022 export const HITU_HEADER = 4; // releaseStep ncont pad[2]

21-export const SYNC_HEADER = 4; // stepA stepB pad[2], then partials[2][G][12]

23+export const SYNC_HEADER = 4; // stepA stepB rawBasis pad, then partials[2][G][12]

24+/** sync[2] when the partials are raw-mode sums; states exported before step 3 hold orthonormal-mode sums there and 0 here. */

25+export const RAW_BASIS = 1;

2226 export const RADCOEF_STRIDE = 36; // coefF coefB i0F i0B tapsF[16] tapsB[16]

2327 export const WORKGROUP_NODES = 256;

2428 export const GATHER_ELEMS = 256; // elements per gather workgroup

@@ -54,6 +58,7 @@ const DT: f32 = ${f(c.dt)}; const DAMP0: f32 = ${f(1 / (1 + 2 * c.sigma0 * c.dt)

5458 const HF: u32 = ${HF}u; const HU: u32 = ${HU}u; const RC: u32 = ${RADCOEF_STRIDE}u;

5559 const OFF_V: u32 = 5u * N; const OFF_VWP: u32 = 10u * N;

5660 const OFF_NHIST: u32 = W * N; const OFF_AHIST: u32 = (W + D) * N; const OFF_ACCOUT: u32 = (W + D + M) * N;

61+const GINV = array<f32, 36>(${Array.from(c.rigidGramInv, f).join(", ")});

5762 const GE: u32 = ${Math.ceil(c.E / GATHER_ELEMS)}u;

5863 const OFF_TAPS: u32 = ${NODE_STRIDE}u * N; const OFF_PART: u32 = ${SYNC_HEADER}u;

5964

@@ -141,6 +146,7 @@ ${bendingProduct()}

141146

142147 var<workgroup> red: array<f32, ${12 * WORKGROUP_NODES}>;

143148 var<workgroup> csum: array<f32, 12>;

149+var<workgroup> coef: array<f32, 12>;

144150

145151 fn modW(x: i32) -> u32 { return u32(((x % i32(W)) + i32(W)) % i32(W)); }

146152

@@ -149,15 +155,24 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca

149155 let nd = gid.x; let t = lid.x;

150156 let step = bitcast<u32>(sync[1]);

151157 let parity = step & 1u; let prevp = (step + 1u) & 1u;

152- // twelve threads sum the previous step's partials of all workgroups; the others read the sums from shared memory

158+ // twelve threads sum the previous step's partials of all workgroups, p = R^T M v and R^T M q with R

159+ // the raw rigid modes, and form the coefficients G^-1 p of the M-orthogonal projection; the others

160+ // read the coefficients from shared memory

153161 if (t < 12u) {

154162 var s = 0.0;

155163 for (var g = 0u; g < G; g++) { s += sync[OFF_PART + (prevp * G + g) * 12u + t]; }

156164 csum[t] = s;

157165 }

158166 workgroupBarrier();

167+ if (t < 12u) {

168+ let r = t % 6u; let o = t - r;

169+ var s = 0.0;

170+ for (var m = 0u; m < 6u; m++) { s += GINV[6u * r + m] * csum[o + m]; }

171+ coef[t] = s;

172+ }

173+ workgroupBarrier();

159174 var c: array<f32, 12>;

160- for (var k = 0u; k < 12u; k++) { c[k] = csum[k]; }

175+ for (var k = 0u; k < 12u; k++) { c[k] = coef[k]; }

161176 var part: array<f32, 12>;

162177 for (var k = 0u; k < 12u; k++) { part[k] = 0.0; }

163178 if (nd < N) {

@@ -178,18 +193,26 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca

178193 }

179194 }

180195 var acc = 0.0;

196+ // the rigid content to remove, sum_k c_k R_k at this node (the modes of rawRigidModes)

197+ let X = nodeF[10u * N + nd]; let Y = nodeF[11u * N + nd]; let Z = nodeF[12u * N + nd];

198+ let cv = array<f32, 5>(c[0] - c[2] * Y + c[5] * Z, c[1] + c[2] * X - c[4] * Z, c[3] + c[4] * Y - c[5] * X, c[4], c[5]);

199+ let cq = array<f32, 5>(c[6] - c[8] * Y + c[11] * Z, c[7] + c[8] * X - c[10] * Z, c[9] + c[10] * Y - c[11] * X, c[10], c[11]);

200+ var mv: array<f32, 5>; var mq: array<f32, 5>;

181201 for (var d = 0u; d < 5u; d++) {

182202 let idx = 5u * nd + d;

183203 var vd = state[OFF_V + idx]; var qd = state[idx];

184- for (var k = 0u; k < 6u; k++) { let r = nodeF[(10u + 5u * k + d) * N + nd]; vd -= c[k] * r; qd -= c[6u + k] * r; }

204+ vd -= cv[d]; qd -= cq[d];

185205 let v1 = (vd + DT * fo[d] * nodeF[d * N + nd]) * DAMP0;

186206 // normal acceleration of the previous step from consecutive projected velocities

187- if (d < 3u) { acc += nodeF[(40u + d) * N + nd] * (vd - state[OFF_VWP + 3u * nd + d]) / DT; state[OFF_VWP + 3u * nd + d] = vd; }

207+ if (d < 3u) { acc += nodeF[(13u + d) * N + nd] * (vd - state[OFF_VWP + 3u * nd + d]) / DT; state[OFF_VWP + 3u * nd + d] = vd; }

188208 let q1 = qd + DT * v1;

189209 state[idx] = q1; state[OFF_V + idx] = v1;

190210 let mm = nodeF[(5u + d) * N + nd];

191- for (var k = 0u; k < 6u; k++) { let r = nodeF[(10u + 5u * k + d) * N + nd] * mm; part[k] += r * v1; part[6u + k] += r * q1; }

211+ mv[d] = mm * v1; mq[d] = mm * q1;

192212 }

213+ // this node's share of R^T M v' and R^T M q'

214+ part[0] = mv[0]; part[1] = mv[1]; part[2] = -Y * mv[0] + X * mv[1]; part[3] = mv[2]; part[4] = -Z * mv[1] + Y * mv[2] + mv[3]; part[5] = Z * mv[0] - X * mv[2] + mv[4];

215+ part[6] = mq[0]; part[7] = mq[1]; part[8] = -Y * mq[0] + X * mq[1]; part[9] = mq[2]; part[10] = -Z * mq[1] + Y * mq[2] + mq[3]; part[11] = Z * mq[0] - X * mq[2] + mq[4];

193216 // decimation of the acceleration of sub-step s = step - 1 into the W accumulator slots of this node;

194217 // slot w belongs to audio sample j0 + w - W/2 + 1, kept at ring position (sample mod W)

195218 if (step >= 1u) {

src/gpu/simulator.tsmodified+32−5View file

@@ -3,10 +3,12 @@

33 // Node with the `webgpu` (Dawn) package, which the tests use.

44 import { HOLD_HEIGHT } from "../model/exciter.ts";

55 import { strainEnergies } from "../model/forces.ts";

6+import { invertMatrix } from "../model/linalg.ts";

67 import { DECIMATION_SPAN, DELAY_TAPS } from "../model/radiation.ts";

8+import { rawRigidModes } from "../model/shell.ts";

79 import type { Setup } from "../model/setup.ts";

810 import type { BackendOptions, BackendState, ChunkResult, Diagnostics, SimulatorBackend } from "../model/simulator.ts";

9-import { CHUNK_WORDS, ELEM_STRIDE, GATHER_ELEMS, HIT_HEADER, HITU_HEADER, NODE_STRIDE, RADCOEF_STRIDE, STATE_FLOATS, SYNC_HEADER, WORKGROUP_NODES, shaderSource } from "./shaders.ts";

11+import { CHUNK_WORDS, ELEM_STRIDE, GATHER_ELEMS, HIT_HEADER, HITU_HEADER, NODE_STRIDE, RADCOEF_STRIDE, RAW_BASIS, STATE_FLOATS, SYNC_HEADER, WORKGROUP_NODES, shaderSource } from "./shaders.ts";

1012

1113 // evaluated lazily: under Node the WebGPU globals are installed after this module is imported

1214 const storageUsage = () => GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC;

@@ -27,6 +29,7 @@ export class GpuSimulator implements SimulatorBackend {

2729 private readonly framesCapacity: number; private readonly dispCapacity: number; private readonly accRows: number;

2830 private readonly stagingPressure: GPUBuffer; private readonly stagingState: GPUBuffer; private readonly stagingHits: GPUBuffer; private readonly stagingHitsU: GPUBuffer;

2931 private readonly stagingMovie: GPUBuffer | null; private readonly stagingDisp: GPUBuffer | null; private readonly stagingAcc: GPUBuffer | null;

32+ private readonly rawFromOrtho: Float64Array; // p = L c: raw-mode partial sums from orthonormal-mode ones (6 x 6)

3033 private readonly radStateFloats: number; // accumulator slots, node history and the current sample's sub-step accelerations: the part of `rad` that is state

3134

3235 constructor(device: GPUDevice, setup: Setup, chunkCapacity = 8192, opts: BackendOptions | number = {}) {

@@ -63,12 +66,24 @@ export class GpuSimulator implements SimulatorBackend {

6366 }

6467 const elemU = new Uint32Array(3 * E); elemU.set(mesh.tri);

6568 const nodeF = new Float32Array(NODE_STRIDE * N + W * M);

69+ // the shader evaluates the raw rigid modes from (X, Y, z0) and projects with the inverse of their

70+ // M-Gram matrix, which is the projection onto the M-orthonormal modes Rm without storing them

71+ const [cx, cy] = sh.rigidCentre;

72+ const gram = new Float64Array(36), rawFromOrtho = new Float64Array(36);

6673 for (let n = 0; n < N; n++) {

6774 const set = (f: number, x: number) => { nodeF[f * N + n] = x; };

75+ const X = mesh.xy[2 * n] - cx, Y = mesh.xy[2 * n + 1] - cy, z = sh.z0Node[n];

6876 for (let d = 0; d < 5; d++) { set(d, sh.Minv[5 * n + d]); set(5 + d, sh.M[5 * n + d]); }

69- for (let k = 0; k < 6; k++) for (let d = 0; d < 5; d++) set(10 + 5 * k + d, sh.Rm[(k * N + n) * 5 + d]);

70- for (let d = 0; d < 3; d++) set(40 + d, sh.normalNode[3 * n + d]);

77+ set(10, X); set(11, Y); set(12, z);

78+ for (let d = 0; d < 3; d++) set(13 + d, sh.normalNode[3 * n + d]);

79+ const r = rawRigidModes(X, Y, z);

80+ for (let a = 0; a < 6; a++) for (let d = 0; d < 5; d++) {

81+ const rm = r[a][d] * sh.M[5 * n + d];

82+ if (rm === 0) continue;

83+ for (let b = 0; b < 6; b++) { gram[6 * a + b] += rm * r[b][d]; rawFromOrtho[6 * a + b] += rm * sh.Rm[(b * N + n) * 5 + d]; }

84+ }

7185 }

86+ this.rawFromOrtho = rawFromOrtho;

7287 nodeF.set(rad.taps, NODE_STRIDE * N);

7388 const adj = new Uint32Array(N + 1 + sh.adjEntries.length); adj.set(sh.adjOffsets); adj.set(sh.adjEntries, N + 1);

7489 const hitsF = new Float32Array(HF * K), hitsU = new Uint32Array(HU * K);

@@ -102,7 +117,7 @@ export class GpuSimulator implements SimulatorBackend {

102117 mk("elemF", elemF); mk("elemU", elemU); mk("nodeF", nodeF); mk("adj", adj); mk("radCoef", radCoef);

103118 mk("hitsF", hitsF); mk("hitsU", hitsU);

104119 mkEmpty("state", STATE_FLOATS * N); mkEmpty("elemForce", 15 * E);

105- mkEmpty("sync", SYNC_HEADER + 2 * G * 12); mkEmpty("rad", this.radStateFloats + this.accRows * N); mkEmpty("elemHist", D * E);

120+ mk("sync", Float32Array.from({ length: SYNC_HEADER + 2 * G * 12 }, (_, i) => (i === 2 ? RAW_BASIS : 0))); mkEmpty("rad", this.radStateFloats + this.accRows * N); mkEmpty("elemHist", D * E);

106121 mkEmpty("pressure", 2 * chunkCapacity * this.GE); // one partial sum per gather workgroup

107122 const uniform = (name: string) => { const b = device.createBuffer({ label: name, size: 4 * CHUNK_WORDS, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST }); this.buffers[name] = b; return b; };

108123 uniform("chunk"); uniform("chunkD");

@@ -114,7 +129,7 @@ export class GpuSimulator implements SimulatorBackend {

114129 this.stagingAcc = this.accRows > 0 ? staging(this.accRows * N * 4) : null;

115130

116131 // ---- pipelines

117- const module = device.createShaderModule({ code: shaderSource({ N, E, K, P, M, W, D, G, T, accRows: this.accRows, dt: setup.dt, sigma0: setup.config.material.sigma0, nl: sh.nl }) });

132+ const module = device.createShaderModule({ code: shaderSource({ N, E, K, P, M, W, D, G, T, rigidGramInv: invertMatrix(gram, 6), accRows: this.accRows, dt: setup.dt, sigma0: setup.config.material.sigma0, nl: sh.nl }) });

118133 type Kind = "read-only-storage" | "storage" | "uniform";

119134 const bgl = (entries: [number, Kind][]) => device.createBindGroupLayout({

120135 entries: entries.map(([binding, type]) => ({ binding, visibility: GPUShaderStage.COMPUTE, buffer: { type } })),

@@ -276,12 +291,24 @@ export class GpuSimulator implements SimulatorBackend {

276291 for (const [name, len] of Object.entries(expect)) {

277292 let a = st.arrays[name];

278293 if (name === "rad" && a?.length === legacyRad) { const full = new Float32Array(len); full.set(a); a = full; }

294+ if (name === "sync" && a?.length === len && a[2] !== RAW_BASIS) a = this.syncToRawBasis(a as Float32Array);

279295 if (!a || a.length !== len) throw new Error(`state array ${name} is missing or has the wrong size`);

280296 this.device.queue.writeBuffer(this.buffers[name], 0, a.buffer, a.byteOffset, a.byteLength);

281297 }

282298 this.step = st.step; this.samplesDone = st.samplesDone;

283299 }

284300

301+ /** A sync array of a state exported before step 3, whose partials are sums for the orthonormal modes, in the raw basis. */

302+ private syncToRawBasis(old: Float32Array): Float32Array {

303+ const out = Float32Array.from(old), L = this.rawFromOrtho;

304+ out[2] = RAW_BASIS;

305+ for (let b = 0; b < 2 * this.G; b++) for (const o of [0, 6]) {

306+ const base = SYNC_HEADER + 12 * b + o;

307+ for (let a = 0; a < 6; a++) { let s = 0; for (let k = 0; k < 6; k++) s += L[6 * a + k] * old[base + k]; out[base + a] = s; }

308+ }

309+ return out;

310+ }

311+

285312 /** The state buffer (q at offset 0, 5 floats per node), for a renderer that reads w directly. */

286313 get stateBuffer(): GPUBuffer { return this.buffers.state; }

287314

src/model/linalg.tsmodified+14−0View file

@@ -33,6 +33,20 @@ export function symmetricEigenvalues(A: Float64Array | number[], n: number): Flo

3333 return ev.sort();

3434 }

3535

36+/** Inverse of an n x n matrix (row-major) by Gauss-Jordan elimination with partial pivoting. */

37+export function invertMatrix(A: ArrayLike<number>, n: number): Float64Array {

38+ const a = Array.from({ length: n }, (_, i) => [...Array.from({ length: n }, (_, j) => A[n * i + j]), ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]);

39+ for (let c = 0; c < n; c++) {

40+ let piv = c;

41+ for (let r = c + 1; r < n; r++) if (Math.abs(a[r][c]) > Math.abs(a[piv][c])) piv = r;

42+ [a[c], a[piv]] = [a[piv], a[c]];

43+ const d = a[c][c];

44+ for (let j = 0; j < 2 * n; j++) a[c][j] /= d;

45+ for (let r = 0; r < n; r++) if (r !== c) { const f = a[r][c]; for (let j = 0; j < 2 * n; j++) a[r][j] -= f * a[c][j]; }

46+ }

47+ return Float64Array.from(a.flatMap((row) => row.slice(n)));

48+}

49+

3650 /** Solve the 3 x 3 system A x = b (row-major A) by Gaussian elimination with pivoting. */

3751 export function solve3(A: number[], b: number[]): number[] {

3852 const M = [A.slice(0, 3).concat(b[0]), A.slice(3, 6).concat(b[1]), A.slice(6, 9).concat(b[2])];

src/model/shell.tsmodified+11−7View file

@@ -20,11 +20,19 @@ export interface Shell {

2020 rhoHA: Float64Array; // (E,) rho h A

2121 M: Float64Array; Minv: Float64Array; // (N, 5)

2222 Rm: Float64Array; // (6, N, 5), M-orthonormal

23+ rigidCentre: [number, number]; // mean node position, the origin of rawRigidModes

2324 adjOffsets: Uint32Array; adjEntries: Uint32Array; // CSR: node -> list of (3 e + local index)

2425 mass: number; maxSlope: number; omegaMax: number; dtStable: number; limitingElement: number;

2526 zeroModes: [number, number]; // min and max count of zero-energy element modes (expect 6)

2627 }

2728

29+/** The six rigid-body modes of the curved shell at a node at (X, Y) from the rigid centre with rest

30+ * height z, before M-orthonormalisation: (u, v, w, theta_x, theta_y) of the two in-plane translations,

31+ * the rotation about z, the normal translation and the two tilts. */

32+export function rawRigidModes(X: number, Y: number, z: number): number[][] {

33+ 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]];

34+}

35+

2836 export function buildShell(mesh: Mesh, profile: CymbalProfile, material: Material, numerics: Numerics): Shell {

2937 const N = mesh.n, E = mesh.e, { xy, tri, area, Gx, Gy } = mesh;

3038 const nl = numerics.nonlinear ? 1 : 0;

@@ -78,12 +86,8 @@ export function buildShell(mesh: Mesh, profile: CymbalProfile, material: Materia

7886 const Rm = new Float64Array(6 * N_DOF * N);

7987 const set = (k: number, n: number, d: number, val: number) => { Rm[(k * N + n) * N_DOF + d] = val; };

8088 for (let n = 0; n < N; n++) {

81- const X = xy[2 * n] - cx, Y = xy[2 * n + 1] - cy, z = z0Node[n];

82- set(0, n, 0, 1); set(1, n, 1, 1);

83- set(2, n, 0, -Y); set(2, n, 1, X);

84- set(3, n, 2, 1);

85- set(4, n, 2, Y); set(4, n, 3, 1); set(4, n, 1, -z);

86- set(5, n, 2, -X); set(5, n, 4, 1); set(5, n, 0, z);

89+ const r = rawRigidModes(xy[2 * n] - cx, xy[2 * n + 1] - cy, z0Node[n]);

90+ for (let k = 0; k < 6; k++) for (let d = 0; d < N_DOF; d++) set(k, n, d, r[k][d]);

8791 }

8892 const L = N_DOF * N;

8993 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; };

@@ -116,7 +120,7 @@ export function buildShell(mesh: Mesh, profile: CymbalProfile, material: Materia

116120 if (top > lamMax) { lamMax = top; limitingElement = e; }

117121 }

118122 const omegaMax = Math.sqrt(lamMax);

119- return { mesh, material, numerics, nl, hE, hNode, z0Node, z0x, z0y, normalNode, De, Cm, Ke, rhoHA, M, Minv, Rm, adjOffsets, adjEntries,

123+ return { mesh, material, numerics, nl, hE, hNode, z0Node, z0x, z0y, normalNode, De, Cm, Ke, rhoHA, M, Minv, Rm, rigidCentre: [cx, cy], adjOffsets, adjEntries,

120124 mass, maxSlope, omegaMax, dtStable: 2 / omegaMax, limitingElement, zeroModes: [zmin, zmax] };

121125 }

122126