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}
119/** The condensation as an operator: given a transverse field w, the statically consistent
120 * in-plane displacements and rotations o = -Koo^-1 Kow w, node-major (u, v, thx, thy) per node.
121 * This is what makes a mode shape from the condensed eigenproblem a complete initial condition
122 * for the five-DOF simulator: with o omitted, releasing a pure-w shape also launches spurious
123 * fast membrane and rotation transients. The assembly repeats what condensedStiffness does but
124 * keeps only Koo (factorised) and the sparse Kow, so building the operator costs a factorisation
125 * and each application one banded solve. */
126export interface CondensationOperator { complement(w: Float64Array): Float64Array }
128export function buildCondensationOperator(mesh: Mesh, shell: Shell, progress: Progress): CondensationOperator {
129 const N = mesh.n, E = mesh.e, tri = mesh.tri;
130 progress("node ordering", 0);
131 const graph = nodeGraph(mesh);
132 const perm = reverseCuthillMcKee(N, graph.ptr, graph.idx);
133 const inv = new Uint32Array(N);
134 for (let p = 0; p < N; p++) inv[perm[p]] = p;
135 let bw = 0;
136 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]]));
137 const b = 4 * bw + 3;
138 const No = 4 * N;
139 const Koo = bandMatrix(No, b);
140 const cooI: number[] = [], cooJ: number[] = [], cooV: number[] = [];
141 const oSlot = [0, 1, -1, 2, 3];
142 const Ke = new Float64Array(225);
143 progress("assembly", 0);
144 for (let e = 0; e < E; e++) {
145 elementTangent(e, mesh.Gx, mesh.Gy, mesh.area[e], shell.z0x[e], shell.z0y[e], shell.Cm, shell.Ke, Ke);
146 for (let a = 0; a < 3; a++) {
147 const pa = inv[tri[3 * e + a]];
148 for (let bb = 0; bb < 3; bb++) {
149 const nb = tri[3 * e + bb], pb = inv[nb];
150 for (let da = 0; da < 5; da++) {
151 if (da === 2) continue;
152 const i = 4 * pa + oSlot[da];
153 for (let db = 0; db < 5; db++) {
154 const v = Ke[15 * (5 * a + da) + (5 * bb + db)];
155 if (v === 0) continue;
156 if (db === 2) {
157 cooI.push(i); cooJ.push(nb); cooV.push(v);
158 } else {
159 const j = 4 * pb + oSlot[db];
160 if (i >= j) Koo.data[i * (b + 1) + (j - i + b)] += v;
161 }
162 }
163 }
164 }
165 }
166 if ((e & 1023) === 0) progress("assembly", e / E);
167 }
168 // pin u0, v0, v1 (in-plane rigid motions), as in condensedStiffness
169 const pins = [4 * inv[0] + 0, 4 * inv[0] + 1, 4 * inv[1] + 1];
170 const w1 = b + 1;
171 for (const p of pins) {
172 for (let j = Math.max(0, p - b); j <= p; j++) Koo.data[p * w1 + (j - p + b)] = 0;
173 for (let i = p + 1; i < Math.min(No, p + b + 1); i++) Koo.data[i * w1 + (p - i + b)] = 0;
174 Koo.data[p * w1 + b] = 1;
175 }
176 const pinned = new Uint8Array(No);
177 for (const p of pins) pinned[p] = 1;
178 for (let k = 0; k < cooI.length; k++) if (pinned[cooI[k]]) cooV[k] = 0;
179 progress("factorisation", 0);
180 if (!bandedCholesky(Koo, (f) => progress("factorisation", f))) throw new Error("condensation: the in-plane/rotation stiffness is not positive definite");
181 progress("done", 1);
182 return {
183 complement(w: Float64Array): Float64Array {
184 const rhs = new Float64Array(No);
185 for (let k = 0; k < cooI.length; k++) rhs[cooI[k]] -= cooV[k] * w[cooJ[k]];
186 bandedSolve(Koo, rhs, 1);
187 const o = new Float64Array(4 * N);
188 for (let n = 0; n < N; n++) for (let s = 0; s < 4; s++) o[4 * n + s] = rhs[4 * inv[n] + s];
189 return o;
190 },
191 };
192}