Voxel stress analysis: elements, loads, reactions and modes
Builds the stiffness of the voxel elements, applies supports and loads, and works out reactions, element stresses, natural frequencies and the mesh's volume and mass.
How the precision is kept
Element matrices, the matrix-free stiffness product, loads, reactions and the modal reduction are 64-bit in Float64Array. The answer's accuracy is set by the method: the voxel staircase mesh, the conjugate gradient tolerance (1e-8 static, 1e-10 inside the modal loop) and the modal stop at a relative change of 1e-6. Small dense problems use Gaussian elimination and a cyclic Jacobi eigenvalue routine with absolute thresholds.
Open item. Give the modal result a converged flag, and scale the Jacobi stopping threshold to the size of the matrix.
Calculations under this policy
15 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Closed form in 64-bit. The factor 1 − 2ν goes to zero as Poisson's ratio approaches 0.5 and λ grows without bound; this function does not check the ratio.
fea.js
lines 40 to 46
PartCraft 0.8.1
/** Isotropic D (6×6, engineering shear strains) for E = 1. */
export function elasticD(nu) {
const l = nu / ((1 + nu) * (1 - 2 * nu)), m = 1 / (2 * (1 + nu)), D = new Float64Array(36);
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) D[i * 6 + j] = l + (i === j ? 2 * m : 0);
for (let i = 3; i < 6; i++) D[i * 6 + i] = m;
return D;
}
Hexahedron stiffness and stress operators
64-bit
The 24 × 24 stiffness of a unit cube and the operators that give stress at its corners and center.
K = Σ_gauss Bᵀ D B × ⅛; K_e = K_cc − K_ci K_ii⁻¹ K_ic; σ = (E ÷ h) S u
Integrated at the eight Gauss points ±1 ÷ √3 with Math.sqrt(3). The incompatible modes are condensed with a small dense solve, and the condensed matrix is averaged with its transpose to remove rounding asymmetry. Computed once per Poisson's ratio and element type and cached.
fea.js
lines 91 to 124
PartCraft 0.8.1
/**
* The element of a unit cube with E = 1 → { Ke: Float64Array(576), S: Float64Array(9·6·24) } (S: stress
* operators at the 8 corners then the centre). For a cube of side h and modulus E: K = E·h·Ke, σ = (E/h)·S·u.
*/
export function hexElement(nu, type = 'h8i') {
const key = nu + type;
if (hexCache.has(key)) return hexCache.get(key);
const D = elasticD(nu), g = 1 / Math.sqrt(3), detw = 1 / 8;
const Kcc = new Float64Array(576), Kci = new Float64Array(24 * 9), Kii = new Float64Array(81);
for (const x of [-g, g]) for (const y of [-g, g]) for (const z of [-g, g]) {
const Bc = bRows(dNc(x, y, z), 8);
atDb(Bc, 24, D, Bc, 24, detw, Kcc);
if (type === 'h8i') { const Bi = bRows(dNi(x, y, z), 3); atDb(Bc, 24, D, Bi, 9, detw, Kci); atDb(Bi, 9, D, Bi, 9, detw, Kii); }
}
let Ke = Kcc, A = null;
if (type === 'h8i') {
// α = A u, A = −Kii⁻¹ Kic ; Ke = Kcc + Kci A
const Kic = new Float64Array(9 * 24); for (let i = 0; i < 24; i++) for (let j = 0; j < 9; j++) Kic[j * 24 + i] = Kci[i * 9 + j];
A = solveSmall(Kii, 9, Kic, 24).map((v) => -v);
Ke = new Float64Array(576);
for (let i = 0; i < 24; i++) for (let j = 0; j < 24; j++) { let s = Kcc[i * 24 + j]; for (let k = 0; k < 9; k++) s += Kci[i * 9 + k] * A[k * 24 + j]; Ke[i * 24 + j] = s; }
for (let i = 0; i < 24; i++) for (let j = 0; j < i; j++) { const m = (Ke[i * 24 + j] + Ke[j * 24 + i]) / 2; Ke[i * 24 + j] = Ke[j * 24 + i] = m; }
}
const S = new Float64Array(9 * 144);
const pts = XI.map((_, a) => [XI[a], ETA[a], ZETA[a]]).concat([[0, 0, 0]]);
pts.forEach(([x, y, z], p) => {
const B = bRows(dNc(x, y, z), 8);
if (A) { const Bi = bRows(dNi(x, y, z), 3); for (let r = 0; r < 6; r++) for (let j = 0; j < 24; j++) { let s = 0; for (let k = 0; k < 9; k++) s += Bi[r * 9 + k] * A[k * 24 + j]; B[r * 24 + j] += s; } }
for (let r = 0; r < 6; r++) for (let j = 0; j < 24; j++) { let s = 0; for (let k = 0; k < 6; k++) s += D[r * 6 + k] * B[k * 24 + j]; S[p * 144 + r * 24 + j] = s; }
});
const out = { Ke, S, type, nu };
hexCache.set(key, out);
return out;
}
Small dense solve
64-bit
Solves a small system by Gaussian elimination with row exchanges.
Partial pivoting on Float64Array copies. The matrix is called singular only when the best pivot is under 1e-300 in absolute value, so a nearly singular system is solved without a warning.
fea.js
lines 70 to 88
PartCraft 0.8.1
/** Dense solve of A X = B (A n×n, B n×m), Gaussian elimination with partial pivoting (small systems). */
export function solveSmall(A, n, B, m) {
A = Float64Array.from(A); B = Float64Array.from(B);
for (let c = 0; c < n; c++) {
let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(A[r * n + c]) > Math.abs(A[p * n + c])) p = r;
if (Math.abs(A[p * n + c]) < 1e-300) throw new Error('singular matrix'); // i18n-ignore: internal (a dense solve on a matrix the callers keep regular)
if (p !== c) { for (let k = 0; k < n; k++) [A[c * n + k], A[p * n + k]] = [A[p * n + k], A[c * n + k]]; for (let k = 0; k < m; k++) [B[c * m + k], B[p * m + k]] = [B[p * m + k], B[c * m + k]]; }
for (let r = c + 1; r < n; r++) {
const f = A[r * n + c] / A[c * n + c]; if (!f) continue;
for (let k = c; k < n; k++) A[r * n + k] -= f * A[c * n + k];
for (let k = 0; k < m; k++) B[r * m + k] -= f * B[c * m + k];
}
}
for (let c = n - 1; c >= 0; c--) for (let k = 0; k < m; k++) {
let s = B[c * m + k]; for (let j = c + 1; j < n; j++) s -= A[c * n + j] * B[j * m + k];
B[c * m + k] = s / A[c * n + c];
}
return B;
}
Stiffness times displacement
64-bit
Applies the whole model's stiffness to a vector, element by element.
y = Σ_e s_e × K_e x_e, s_e = E h × scale_e
Each element contributes a 24-term 64-bit sum per row, added into the result in element order. No global matrix is stored.
fea.js
lines 145 to 162
PartCraft 0.8.1
/** y = K x (scale: per-model-element factor or null). */
export function applyK(model, x, y, scale = null) {
const { node0, off, Ke, kScale } = model, ne = node0.length;
const xe = new Float64Array(24), dof = new Int32Array(24);
y.fill(0);
for (let m = 0; m < ne; m++) {
const s = scale ? scale[m] * kScale : kScale;
if (s === 0) continue;
const n0 = node0[m];
for (let c = 0; c < 8; c++) { const d = 3 * (n0 + off[c]); dof[3 * c] = d; dof[3 * c + 1] = d + 1; dof[3 * c + 2] = d + 2; xe[3 * c] = x[d]; xe[3 * c + 1] = x[d + 1]; xe[3 * c + 2] = x[d + 2]; }
for (let i = 0, r = 0; i < 24; i++, r += 24) {
let t = 0;
for (let j = 0; j < 24; j++) t += Ke[r + j] * xe[j];
y[dof[i]] += s * t;
}
}
return y;
}
Node-block preconditioner
64-bit
Inverts each node's 3 × 3 diagonal block of the stiffness.
A⁻¹ = adj(A) ÷ det(A)
The inverse is by cofactors and the determinant. A block whose determinant is zero or not a number is set to zero. It is used only as a preconditioner, so it changes the iteration count, not the tolerance the solution meets.
fea.js
lines 164 to 182
PartCraft 0.8.1
/** Node-block diagonal of K, inverted → Float64Array(9·nn) (inactive nodes 0). */
export function blockJacobi(model, scale = null) {
const { node0, off, Ke, kScale, nn } = model, B = new Float64Array(9 * nn);
for (let m = 0; m < node0.length; m++) {
const s = scale ? scale[m] * kScale : kScale;
for (let c = 0; c < 8; c++) {
const b = 9 * (node0[m] + off[c]);
for (let a = 0; a < 3; a++) for (let d = 0; d < 3; d++) B[b + a * 3 + d] += s * Ke[(3 * c + a) * 24 + 3 * c + d];
}
}
for (let n = 0; n < nn; n++) {
const b = 9 * n, a = B.subarray(b, b + 9);
const det = a[0] * (a[4] * a[8] - a[5] * a[7]) - a[1] * (a[3] * a[8] - a[5] * a[6]) + a[2] * (a[3] * a[7] - a[4] * a[6]);
if (!(Math.abs(det) > 0)) { a.fill(0); continue; }
const inv = [a[4] * a[8] - a[5] * a[7], a[2] * a[7] - a[1] * a[8], a[1] * a[5] - a[2] * a[4], a[5] * a[6] - a[3] * a[8], a[0] * a[8] - a[2] * a[6], a[2] * a[3] - a[0] * a[5], a[3] * a[7] - a[4] * a[6], a[1] * a[6] - a[0] * a[7], a[0] * a[4] - a[1] * a[3]];
for (let k = 0; k < 9; k++) a[k] = inv[k] / det;
}
return B;
}
Constraint projections
64-bit
Turns each node's held directions into a projection that removes them.
P = I − Σ q qᵀ over the orthonormal held directions q
Directions are made orthonormal by Gram-Schmidt; a direction with less than 1e-6 left after removing the earlier ones is dropped. The threshold is absolute and assumes directions of about unit length.
fea.js
lines 186 to 203
PartCraft 0.8.1
/** Constraint set → { nodes: Int32Array, P: Float64Array(9 each), dirs: [[unit vectors]] } (P = I − QQᵀ). */
export function makeConstraints(nodeDirs) {
const nodes = [], P = [], dirs = [];
for (const [n, list] of nodeDirs) {
const q = [];
for (const d of list) { // Gram–Schmidt
let v = [...d];
for (const u of q) { const t = v[0] * u[0] + v[1] * u[1] + v[2] * u[2]; v = [v[0] - t * u[0], v[1] - t * u[1], v[2] - t * u[2]]; }
const l = Math.hypot(...v);
if (l > 1e-6) q.push(v.map((c) => c / l));
if (q.length === 3) break;
}
const p = [1, 0, 0, 0, 1, 0, 0, 0, 1];
for (const u of q) for (let a = 0; a < 3; a++) for (let b = 0; b < 3; b++) p[a * 3 + b] -= u[a] * u[b];
nodes.push(n); P.push(...p); dirs.push(q);
}
return { nodes: Int32Array.from(nodes), P: Float64Array.from(P), dirs };
}
Eigenvalues of a small symmetric matrix
64-bit
Eigenvalues and vectors by cyclic Jacobi rotations.
The rotation uses the form that avoids subtracting nearly equal numbers. Sweeps stop when the sum of squared off-diagonal entries is under 1e-30 or after 100 sweeps; that threshold is absolute, not scaled to the matrix, and no flag reports which stop was reached.
fea.js
lines 217 to 233
PartCraft 0.8.1
/** Symmetric eigenvalues of a small matrix (cyclic Jacobi) → { values, vectors (column-major n×n) }. */
export function symEig(A, n) {
const a = Float64Array.from(A), v = new Float64Array(n * n);
for (let i = 0; i < n; i++) v[i * n + i] = 1;
for (let sweep = 0; sweep < 100; sweep++) {
let off = 0; for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += a[i * n + j] ** 2;
if (off < 1e-30) break;
for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
const apq = a[p * n + q]; if (Math.abs(apq) < 1e-300) continue;
const th = (a[q * n + q] - a[p * n + p]) / (2 * apq), t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1)), c = 1 / Math.sqrt(t * t + 1), s = t * c;
for (let k = 0; k < n; k++) { const akp = a[k * n + p], akq = a[k * n + q]; a[k * n + p] = c * akp - s * akq; a[k * n + q] = s * akp + c * akq; }
for (let k = 0; k < n; k++) { const apk = a[p * n + k], aqk = a[q * n + k]; a[p * n + k] = c * apk - s * aqk; a[q * n + k] = s * apk + c * aqk; }
for (let k = 0; k < n; k++) { const vkp = v[k * n + p], vkq = v[k * n + q]; v[k * n + p] = c * vkp - s * vkq; v[k * n + q] = s * vkp + c * vkq; }
}
}
return { values: Array.from({ length: n }, (_, i) => a[i * n + i]), vectors: v };
}
Is the body held
64-bit
Checks that the supports stop all six rigid motions.
Builds a 6 × 6 matrix from the held directions, with node positions measured from the grid center and divided by the grid size. A motion is reported free when its eigenvalue is at most 1e-9 of the largest or at most 1e-12.
fea.js
lines 235 to 259
PartCraft 0.8.1
/**
* Is the body held? The six rigid motions (about the grid's centre, lengths scaled by its size) restricted
* to the constrained directions must be independent. → null, or the motions left free in words.
*/
export function rigidCheck(model, cons) {
const g = model.grid, c = [g.origin[0] + g.nx * g.h / 2, g.origin[1] + g.ny * g.h / 2, g.origin[2] + g.nz * g.h / 2], L = Math.max(g.nx, g.ny, g.nz) * g.h;
const G = new Float64Array(36);
for (let m = 0; m < cons.nodes.length; m++) {
const p = VX.nodeXYZ(g, cons.nodes[m]).map((v, a) => (v - c[a]) / L);
const R = [[1, 0, 0, 0, p[2], -p[1]], [0, 1, 0, -p[2], 0, p[0]], [0, 0, 1, p[1], -p[0], 0]]; // rows x,y,z; cols tx ty tz rx ry rz
for (const q of cons.dirs[m]) {
const r = [0, 1, 2, 3, 4, 5].map((k) => q[0] * R[0][k] + q[1] * R[1][k] + q[2] * R[2][k]);
for (let i = 0; i < 6; i++) for (let j = 0; j < 6; j++) G[i * 6 + j] += r[i] * r[j];
}
}
const { values, vectors } = symEig(G, 6);
const big = Math.max(...values, 1e-30), names = [t('study.fea.free.move_x'), t('study.fea.free.move_y'), t('study.fea.free.move_z'), t('study.fea.free.turn_x'), t('study.fea.free.turn_y'), t('study.fea.free.turn_z')];
const free = [];
values.forEach((v, k) => {
if (v > 1e-9 * big && v > 1e-12) return;
let best = 0; for (let i = 1; i < 6; i++) if (Math.abs(vectors[i * 6 + k]) > Math.abs(vectors[best * 6 + k])) best = i;
free.push(names[best]);
});
return free.length ? free : null;
}
Gravity load
64-bit
Shares each element's weight equally to its eight corners.
corner force = scale × ρ g h³ ÷ 8 × unit direction
64-bit; the direction is divided by its length (1 if the length is zero). When the load carries no g, G_MM from materials.js is used, which is 9810 mm/s²; standard gravity is 9806.65 mm/s².
fea.js
lines 359 to 366
PartCraft 0.8.1
if (l.type === 'gravity') {
const d = l.dir || [0, -1, 0], dl = Math.hypot(...d) || 1, w = density * (l.g ?? G_MM) * g.h ** 3 / 8;
for (let m = 0; m < model.node0.length; m++) {
const s = scale ? scale[m] : 1;
for (let c = 0; c < 8; c++) { const n = model.node0[m] + model.off[c]; for (let a = 0; a < 3; a++) { f[3 * n + a] += s * w * d[a] / dl; total[a] += s * w * d[a] / dl; } }
}
groups.loads.push({ type: 'gravity', nodes: new Int32Array(0) });
continue;
Stress at element corners and center
64-bit
The six stress components and von Mises stress at each element's eight corners and center, and the largest center value.
σ = scale × (E ÷ h) × S u_e − scale × E α ΔT ÷ (1 − 2ν) on the direct components
Stresses are 64-bit sums. The largest element value is compared before it is stored in the 32-bit vmElem array, so the reported maximum element stress is 64-bit. Corner values are summed per node in a Float64Array for the nodal average.
fea.js
lines 394 to 410
PartCraft 0.8.1
const thK = thermal ? model.E * thermal.alpha / (1 - 2 * model.nu) : 0;
let elMax = 0, elMaxAt = -1;
for (let m = 0; m < ne; m++) {
const n0 = node0[m], k = (scale ? scale[m] : 1) * sF;
for (let c = 0; c < 8; c++) { const d = 3 * (n0 + off[c]); ue[3 * c] = u[d]; ue[3 * c + 1] = u[d + 1]; ue[3 * c + 2] = u[d + 2]; }
for (let p = 0; p < 9; p++) {
for (let r = 0; r < 6; r++) { let t = 0; const b = p * 144 + r * 24; for (let j = 0; j < 24; j++) t += S[b + j] * ue[j]; s[r] = k * t; }
if (thermal) {
let dt = 0; if (p < 8) dt = thermal.dT[n0 + off[p]]; else { for (let c = 0; c < 8; c++) dt += thermal.dT[n0 + off[c]]; dt /= 8; }
const sc = scale ? scale[m] : 1, st = sc * thK * dt; s[0] -= st; s[1] -= st; s[2] -= st;
const C = thermal.corr; // h8i: the bubble modes' thermal share per unit E·α (thermal.js stressCorrection)
if (C) for (let r = 0; r < 6; r++) { let t = 0; for (let c = 0; c < 8; c++) t += C[p * 48 + r * 8 + c] * thermal.dT[n0 + off[c]]; s[r] += sc * model.E * thermal.alpha * t; }
}
const vm = vonMises(s);
if (p < 8) { vmSum[n0 + off[p]] += vm; cnt[n0 + off[p]]++; } else { vmElem[m] = vm; if (vm > elMax) { elMax = vm; elMaxAt = m; } }
}
}
Reactions and compliance
64-bit
The support reactions and the work of the loads.
R = (K u − f) projected on the held directions; compliance = Σ fᵢ uᵢ
K u is recomputed in 64-bit and f subtracted at the held nodes; the total is a running sum over those nodes. The reactions inherit the solver's residual of 1e-8 of the load.
fea.js
lines 420 to 431
PartCraft 0.8.1
// reactions: (K u − f) in the constrained directions
const ku = applyK(model, u, new Float64Array(3 * nn), scale);
const reaction = [0, 0, 0], nodeR = new Map();
for (let m = 0; m < cons.nodes.length; m++) {
const d = 3 * cons.nodes[m], r = [ku[d] - f[d], ku[d + 1] - f[d + 1], ku[d + 2] - f[d + 2]], rr = [0, 0, 0];
for (const q of cons.dirs[m]) { const t = r[0] * q[0] + r[1] * q[1] + r[2] * q[2]; for (let a = 0; a < 3; a++) rr[a] += t * q[a]; }
for (let a = 0; a < 3; a++) reaction[a] += rr[a];
nodeR.set(cons.nodes[m], rr);
}
// per constraint (a node shared by two constraints counts in both)
const groupReactions = groups ? groups.map((gr) => { const s = [0, 0, 0]; for (const n of gr.nodes) { const rr = nodeR.get(n); if (rr) for (let a = 0; a < 3; a++) s[a] += rr[a]; } return s; }) : null;
let compliance = 0; for (let i = 0; i < f.length; i++) compliance += f[i] * u[i];
Field value at a point
64-bit
Interpolates a nodal result at a point inside the grid.
value = Σ wᵢ vᵢ ÷ Σ wᵢ over the active corners of the cell
Trilinear weights in 64-bit, renormalized over the corners that belong to the solid. Where the point's own cell has no solid corner the search widens up to two cells, which is a nearest-value estimate, not an interpolation.
fea.js
lines 441 to 463
PartCraft 0.8.1
/** A nodal grid field sampled at a point: trilinear over the cell's active corners (renormalised); the
* nearest cell with active corners when the point's own has none (the staircase). comps 1 or 3. */
export function sampleField(model, values, comps, p, out = new Float64Array(comps)) {
const g = model.grid, NX = g.nx + 1, NXY = NX * (g.ny + 1);
const fx = (p[0] - g.origin[0]) / g.h, fy = (p[1] - g.origin[1]) / g.h, fz = (p[2] - g.origin[2]) / g.h;
for (let ring = 0; ring <= 2; ring++) {
out.fill(0); let wsum = 0;
const i0 = Math.min(g.nx - 1, Math.max(0, Math.floor(fx))), j0 = Math.min(g.ny - 1, Math.max(0, Math.floor(fy))), k0 = Math.min(g.nz - 1, Math.max(0, Math.floor(fz)));
for (let dk = -ring; dk <= ring + 1; dk++) for (let dj = -ring; dj <= ring + 1; dj++) for (let di = -ring; di <= ring + 1; di++) {
const i = i0 + di, j = j0 + dj, k = k0 + dk;
if (i < 0 || j < 0 || k < 0 || i > g.nx || j > g.ny || k > g.nz) continue;
const n = i + NX * j + NXY * k;
if (!model.active[n]) continue;
const w = Math.max(0, 1 - Math.abs(fx - i) / (ring + 1)) * Math.max(0, 1 - Math.abs(fy - j) / (ring + 1)) * Math.max(0, 1 - Math.abs(fz - k) / (ring + 1));
if (w <= 0) continue;
for (let c = 0; c < comps; c++) out[c] += w * values[comps * n + c];
wsum += w;
}
if (wsum > 0) { for (let c = 0; c < comps; c++) out[c] /= wsum; return out; }
}
out.fill(NaN);
return out;
}
Lumped mass
64-bit
The mass at each node for the modal analysis.
node mass = Σ over its elements of ρ h³ ÷ 8
One 64-bit value per element corner, added per node.
fea.js
lines 516 to 521
PartCraft 0.8.1
/** Lumped mass per node (t), 3 equal dofs. */
export function lumpedMass(model, density) {
const m = new Float64Array(model.nn), me = density * model.grid.h ** 3 / 8;
for (let e = 0; e < model.node0.length; e++) for (let c = 0; c < 8; c++) m[model.node0[e] + model.off[c]] += me;
return m;
}
Natural frequencies
64-bit
The lowest natural frequencies and mode shapes.
f = √λ ÷ (2π); stop when |λ − λ(previous)| ≤ 1e-6 × λ for every wanted mode
Subspace iteration with inner conjugate gradient solves at 1e-10, whose converged flag is not read. The eigenvalues are accepted at a relative change of 1e-6, so a frequency is settled to about 6 digits, and after 40 iterations the last estimate is returned without a flag. A reduced mass pivot at or below zero is replaced by 1e-300 before the square root.
fea.js
lines 523 to 570
PartCraft 0.8.1
/**
* The lowest `modes` natural frequencies (Hz) by subspace iteration: X ← K⁻¹ M X, Rayleigh–Ritz on
* (XᵀKX, XᵀMX), until every wanted eigenvalue changes < 1e-6 (relative). → { freqs, shapes: [Float64Array] }.
*/
export function modal(model, cons, density, { modes = 4, maxIter = 40, onProgress = null, cancelled = null } = {}) {
const n = 3 * model.nn, q = Math.min(2 * modes, modes + 8), mass = lumpedMass(model, density), proj = projector(model, cons);
const Minv = blockJacobi(model);
// start vectors: the mass diagonal and pseudo-random vectors
let rnd = 12345; const rand = () => ((rnd = (rnd * 1103515245 + 12345) % 2147483648) / 2147483648) - 0.5;
let X = Array.from({ length: q }, (_, k) => { const v = new Float64Array(n); for (let i = 0; i < n; i++) v[i] = k === 0 ? 1 : rand(); return proj(v); });
const Mx = (x) => { const y = new Float64Array(n); for (let i = 0; i < n; i++) y[i] = mass[Math.floor(i / 3)] * x[i]; return y; };
let prev = null, lam = null, Y = X.map(() => null);
for (let it = 0; it < maxIter; it++) {
const Z = X.map((x, k) => pcg(model, cons, Mx(x), { x0: Y[k], tol: 1e-10, Minv, cancelled }).u);
Y = Z;
// Rayleigh–Ritz
const KZ = Z.map((z) => proj(applyK(model, z, new Float64Array(n)))), MZ = Z.map(Mx);
const Kr = new Float64Array(q * q), Mr = new Float64Array(q * q);
for (let a = 0; a < q; a++) for (let b = a; b < q; b++) {
let sk = 0, sm = 0; for (let i = 0; i < n; i++) { sk += Z[a][i] * KZ[b][i]; sm += Z[a][i] * MZ[b][i]; }
Kr[a * q + b] = Kr[b * q + a] = sk; Mr[a * q + b] = Mr[b * q + a] = sm;
}
// Mr = L Lᵀ (Cholesky), C = L⁻¹ Kr L⁻ᵀ, eig C
const L = new Float64Array(q * q);
for (let i = 0; i < q; i++) for (let j = 0; j <= i; j++) {
let s = Mr[i * q + j]; for (let k = 0; k < j; k++) s -= L[i * q + k] * L[j * q + k];
L[i * q + j] = i === j ? Math.sqrt(Math.max(s, 1e-300)) : s / L[j * q + j];
}
const Li = solveSmall(L, q, Float64Array.from({ length: q * q }, (_, k) => (Math.floor(k / q) === k % q ? 1 : 0)), q);
const C = new Float64Array(q * q);
for (let i = 0; i < q; i++) for (let j = 0; j < q; j++) { let s = 0; for (let a = 0; a < q; a++) for (let b = 0; b < q; b++) s += Li[i * q + a] * Kr[a * q + b] * Li[j * q + b]; C[i * q + j] = s; }
const { values, vectors } = symEig(C, q);
const order = values.map((v, i) => i).sort((a, b) => values[a] - values[b]);
lam = order.map((i) => values[i]);
// new X = Z · L⁻ᵀ · vectors (sorted), normalised
X = order.map((col) => {
const w = new Float64Array(q); for (let a = 0; a < q; a++) { let s = 0; for (let b = 0; b < q; b++) s += Li[b * q + a] * vectors[b * q + col]; w[a] = s; }
const x = new Float64Array(n); for (let a = 0; a < q; a++) { const c = w[a]; if (c) for (let i = 0; i < n; i++) x[i] += c * Z[a][i]; }
return x;
});
Y = X.map((x, k) => Float64Array.from(x, (v) => v / lam[k])); // warm start: K⁻¹M x ≈ x/λ for a converged x
if (onProgress) onProgress((it + 1) / maxIter, `Modes: iteration ${it + 1}`); // i18n-ignore: progress label, never shown
if (prev && lam.slice(0, modes).every((l, i) => Math.abs(l - prev[i]) <= 1e-6 * l)) break;
prev = lam.slice();
}
const freqs = lam.slice(0, modes).map((l) => Math.sqrt(Math.max(l, 0)) / (2 * Math.PI));
return { freqs, shapes: X.slice(0, modes) };
}
Mesh volume and mass
64-bit
The volume and mass reported with a study.
volume = elements × h³; mass (kg) = volume × ρ (t/mm³) × 1000
64-bit arithmetic, but the volume is that of the solid voxels, not of the part's exact shape.
8-node hexahedron with 2 × 2 × 2 Gauss integration; the default element adds Wilson-Taylor incompatible modes condensed out once (file header). Constraints as projections. Modal: lumped mass, subspace iteration with Rayleigh-Ritz.
Number types
Float64Array for matrices, loads and displacements; Int32Array for indexes.
Tolerances
Static solve: residual under 1e-8 of the load (the job can change it). Modal: inner solves 1e-10, eigenvalue change 1e-6 relative, at most 40 iterations. Jacobi eigenvalues: sum of squared off-diagonals under 1e-30 (absolute) or 100 sweeps. Small dense solve: pivot under 1e-300 is singular. Rigid motion free: eigenvalue at most 1e-9 of the largest, or at most 1e-12. Constraint directions shorter than 1e-6 after Gram-Schmidt are dropped.
Limits
Curved faces lose up to about half a cell; stress at re-entrant corners and constrained edges depends on the mesh (file header). The reported volume and mass are those of the voxels. The modal routine returns its last estimate when 40 iterations pass without meeting the tolerance. Gravity uses 9810 mm/s² from materials.js when the load gives no value.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.