ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 1// Assembly of the linearised shell stiffness and static condensation onto the transverse DOFs.
2//
3// The tangent stiffness at rest couples five DOFs per node. The in-plane DOFs (u, v) carry
4// stiff membrane waves and the DKT rotations (theta_x, theta_y) carry no physical inertia worth
5// keeping, so both are eliminated statically: with q = (w, o) and K = [[Kww, Kwo], [Kow, Koo]],
6// the condensed stiffness on w is S = Kww - Kwo Koo^-1 Kow. Koo is factorised as a banded
7// Cholesky after a reverse Cuthill-McKee ordering of the nodes; three in-plane DOFs are pinned
8// to remove the in-plane rigid motions. S keeps the transverse translation and the two tilts as
9// (numerically) zero-frequency modes, which the caller drops.
10import { elementTangent, type Shell } from "../model/shell.ts";
11import type { Mesh } from "../model/mesh.ts";
12import { bandMatrix, bandedCholesky, bandedSolve, reverseCuthillMcKee, type BandMatrix } from "./banded.ts";
14export type Progress = (phase: string, frac: number) => void;
16/** Node adjacency (CSR, sorted, without self loops) from the triangles. */
17export function nodeGraph(mesh: Mesh): { ptr: Uint32Array; idx: Uint32Array } {
18 const N = mesh.n, E = mesh.e, tri = mesh.tri;
19 const sets: Set<number>[] = Array.from({ length: N }, () => new Set<number>());
20 for (let e = 0; e < E; e++) {
21 const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
22 sets[a].add(b); sets[a].add(c); sets[b].add(a); sets[b].add(c); sets[c].add(a); sets[c].add(b);
23 }
24 const ptr = new Uint32Array(N + 1);
25 for (let n = 0; n < N; n++) ptr[n + 1] = ptr[n] + sets[n].size;
26 const idx = new Uint32Array(ptr[N]);
27 for (let n = 0; n < N; n++) idx.set(Array.from(sets[n]).sort((p, q) => p - q), ptr[n]);
28 return { ptr, idx };
29}
31function bandZeroRowCol(A: BandMatrix, p: number): void {
32 const { n, b, data } = A;
33 const w = b + 1;
34 for (let j = Math.max(0, p - b); j <= p; j++) data[p * w + (j - p + b)] = 0;
35 for (let i = p + 1; i < Math.min(n, p + b + 1); i++) data[i * w + (p - i + b)] = 0;
36 data[p * w + b] = 1;
37}
39/** The condensed transverse stiffness S (row-major N x N, symmetric). */
40export function condensedStiffness(mesh: Mesh, shell: Shell, progress: Progress): Float64Array {
41 const N = mesh.n, E = mesh.e, tri = mesh.tri;
42 progress("node ordering", 0);
43 const graph = nodeGraph(mesh);
44 const perm = reverseCuthillMcKee(N, graph.ptr, graph.idx);
45 const inv = new Uint32Array(N);
46 for (let p = 0; p < N; p++) inv[perm[p]] = p;
47 let bw = 0;
48 for (let n = 0; n < N; n++) for (let k = graph.ptr[n]; k < graph.ptr[n + 1]; k++) bw = Math.max(bw, Math.abs(inv[n] - inv[graph.idx[k]]));
49 const b = 4 * bw + 3;
50 const No = 4 * N;
51 const Koo = bandMatrix(No, b);
52 const R = new Float64Array(No * N); // K_ow, later overwritten by X = Koo^-1 K_ow
53 const S = new Float64Array(N * N); // K_ww, later S
54 const cooI: number[] = [], cooJ: number[] = [], cooV: number[] = [];
56 // DOF slot d of node n: 2 = w; the others map to o = 4 p + {0, 1, 2, 3} for (u, v, thx, thy)
57 const oSlot = [0, 1, -1, 2, 3];
58 const Ke = new Float64Array(225);
59 progress("assembly", 0);
60 for (let e = 0; e < E; e++) {
61 elementTangent(e, mesh.Gx, mesh.Gy, mesh.area[e], shell.z0x[e], shell.z0y[e], shell.Cm, shell.Ke, Ke);
62 for (let a = 0; a < 3; a++) {
63 const na = tri[3 * e + a], pa = inv[na];
64 for (let bb = 0; bb < 3; bb++) {
65 const nb = tri[3 * e + bb], pb = inv[nb];
66 for (let da = 0; da < 5; da++) {
67 for (let db = 0; db < 5; db++) {
68 const v = Ke[15 * (5 * a + da) + (5 * bb + db)];
69 if (v === 0) continue;
70 const wa = da === 2, wb = db === 2;
71 if (wa && wb) {
72 S[na * N + nb] += v;
73 } else if (!wa && !wb) {
74 const i = 4 * pa + oSlot[da], j = 4 * pb + oSlot[db];
75 if (i >= j) Koo.data[i * (b + 1) + (j - i + b)] += v;
76 } else if (!wa && wb) {
77 const i = 4 * pa + oSlot[da];
78 R[i * N + nb] += v;
79 cooI.push(i); cooJ.push(nb); cooV.push(v);
80 }
81 }
82 }
83 }
84 }
85 if ((e & 1023) === 0) progress("assembly", e / E);
86 }
87 // pin u0, v0, v1 (in-plane rigid motions)
88 const pins = [4 * inv[0] + 0, 4 * inv[0] + 1, 4 * inv[1] + 1];
89 for (const p of pins) {
90 bandZeroRowCol(Koo, p);
91 R.fill(0, p * N, (p + 1) * N);
92 }
93 const pinned = new Uint8Array(No);
94 for (const p of pins) pinned[p] = 1;
95 for (let k = 0; k < cooI.length; k++) if (pinned[cooI[k]]) cooV[k] = 0;
97 progress("factorisation", 0);
98 if (!bandedCholesky(Koo, (f) => progress("factorisation", f))) throw new Error("condensation: the in-plane/rotation stiffness is not positive definite");
99 progress("condensation solve", 0);
100 bandedSolve(Koo, R, N, (f) => progress("condensation solve", f));
101 progress("Schur complement", 0);
102 for (let k = 0; k < cooI.length; k++) {
103 const v = cooV[k];
104 if (v === 0) continue;
105 const rowX = cooI[k] * N, rowS = cooJ[k] * N;
106 for (let c = 0; c < N; c++) S[rowS + c] -= v * R[rowX + c];
107 if ((k & 65535) === 0) progress("Schur complement", k / cooI.length);
108 }
109 for (let i = 0; i < N; i++) {
110 for (let j = 0; j < i; j++) {
111 const s = 0.5 * (S[i * N + j] + S[j * N + i]);
112 S[i * N + j] = s;
113 S[j * N + i] = s;
114 }
115 }
116 return S;
117}