Assembles the stiffness of the tetrahedral model, applies supports and loads, and solves for displacements, natural frequencies and the pre-stress of a buckling study.
How the precision is kept
Matrices and vectors are Float64Array; the stiffness is kept as its lower triangle in compressed rows with Int32Array indexes. Displacements come from preconditioned conjugate gradients that stop at a relative residual, so the solution carries the iteration tolerance (1e-8 by default for a static study) rather than 15 digits. Natural frequencies come from an iterative eigen-solver that stops at a relative residual of 1e-4.
Calculations under this policy
12 calculations in 2 source files. 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.
Adds every element's stiffness (and mass) into the model's sparse matrices and totals the volume.
K = Σₑ Kₑ; volume = Σₑ volumeₑ
Entries are added element by element in mesh order into Float64Array values of the lower triangle. The volume and the nodal ∫N dV used for gravity are plain running sums in the same order.
tetfea.js
lines 212 to 240
PartCraft 0.8.1
/** Assemble K (lower CSR values) → { K: Float64Array, vol }. Also the node mass (scalar lower graph) when rho given. */
export function assemble(model, E, nu, { rho = 0, mass = false, onProgress = null } = {}) {
const D = elasticD(nu), K = new Float64Array(model.col.length), m = model.tets10.length / 10;
const Ke = new Float64Array(900), el = new Int32Array(10);
const Mv = mass ? new Float64Array(model.nodeAdj.length) : null, Me = new Float64Array(100), Ni = new Float64Array(10), nodeVolN = new Float64Array(model.n);
let vol = 0;
const pos = new Int32Array(100); // node-pair → first scalar entry offsets (per element)
let t0 = Date.now();
for (let e = 0; e < m; e++) {
for (let a = 0; a < 10; a++) el[a] = model.tets10[10 * e + a];
vol += t10Stiffness(model.X, el, E, D, Ke);
for (let a = 0; a < 10; a++) for (let b = 0; b < 10; b++) {
const A = el[a], B = el[b];
if (B > A) continue;
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) {
if (B === A && j > i) continue;
K[entry(model, A, i, B, j)] += Ke[(3 * a + i) * 30 + 3 * b + j];
}
}
if (mass || rho) {
t10Mass(model.X, el, rho || 1, Me, Ni);
for (let a = 0; a < 10; a++) nodeVolN[el[a]] += Ni[a];
if (mass) for (let a = 0; a < 10; a++) for (let b = 0; b < 10; b++) { const A = el[a], B = el[b]; if (B > A) continue; Mv[nodeEntry(model, A, B)] += Me[a * 10 + b]; }
}
if (onProgress && Date.now() - t0 > 200) { t0 = Date.now(); onProgress(e / m); }
}
void pos;
return { K, vol, M: Mv, nodeVolN };
}
Matrix times vector
64-bit
Multiplies the symmetric stiffness, stored as its lower triangle, by a vector.
y = K x
Each row's products are summed in column order and the mirrored upper-triangle terms are added into y as they are met. Plain 64-bit sums, no compensation.
tetfea.js
lines 259 to 269
PartCraft 0.8.1
/** y = A x for a symmetric matrix stored as lower scalar CSR (values V). */
export function symMul(model, V, x, y) {
const { rp, col, N } = model;
y.fill(0);
for (let r = 0; r < N; r++) {
let s = 0; const xr = x[r];
for (let k = rp[r], e = rp[r + 1] - 1; k < e; k++) { const c = col[k], v = V[k]; s += v * x[c]; y[c] += v * xr; }
y[r] += s + V[rp[r + 1] - 1] * xr; // the diagonal is the row's last entry
}
return y;
}
Constrained stiffness
64-bit
Builds the stiffness with supports applied.
K̃ = P K P + α (I − P)
Only the 3 × 3 blocks that touch a supported node are multiplied by the projections, on a copy of the Float64Array values. α is the mean of the diagonal of K, so the supported directions get a stiffness of the same size as the rest.
tetfea.js
lines 299 to 316
PartCraft 0.8.1
/** K̃ = PKP + α(I − P) as lower CSR values (a copy). */
function projectedK(model, K, pj, alpha) {
const Kt = Float64Array.from(K), { has, P } = pj, { n, nodePtr, nodeAdj } = model;
const blk = new Float64Array(9), t = new Float64Array(9), I3 = [1, 0, 0, 0, 1, 0, 0, 0, 1];
const Pm = (a) => (has[a] ? P.subarray(9 * a, 9 * a + 9) : I3);
for (let a = 0; a < n; a++) for (let k = nodePtr[a]; k < nodePtr[a + 1]; k++) {
const b = nodeAdj[k];
if (!has[a] && !has[b]) continue;
const diag = a === b;
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) blk[3 * i + j] = diag && j > i ? Kt[entry(model, a, j, a, i)] : Kt[entry(model, a, i, b, j)];
const Pa = Pm(a), Pb = Pm(b);
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) { let s = 0; for (let k2 = 0; k2 < 3; k2++) s += blk[3 * i + k2] * Pb[3 * k2 + j]; t[3 * i + j] = s; }
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) { let s = 0; for (let k2 = 0; k2 < 3; k2++) s += Pa[3 * i + k2] * t[3 * k2 + j]; blk[3 * i + j] = s; }
if (diag) for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) blk[3 * i + j] += alpha * ((i === j ? 1 : 0) - Pa[3 * i + j]);
for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) { if (diag && j > i) continue; Kt[entry(model, a, i, b, j)] = blk[3 * i + j]; }
}
return Kt;
}
Scale of the supported directions
64-bit
The value α put on the diagonal for supported directions.
α = Σ Kᵢᵢ ÷ N
One running sum of the N diagonal entries, then one division.
tetfea.js
lines 362 to 364
PartCraft 0.8.1
let dsum = 0; for (let r = 0; r < model.N; r++) dsum += K[model.rp[r + 1] - 1];
const alpha = dsum / model.N;
const Kt = projectedK(model, K, pj, alpha);
Incomplete Cholesky preconditioner
64-bit
Factors the constrained stiffness on its own sparsity pattern to speed up the iteration.
L Lᵀ ≈ K̃ on the pattern of K̃
A pivot must be greater than 1e-14 × the matching diagonal entry of the matrix; otherwise the whole factorization restarts with the diagonal multiplied by 1 + shift (1e-3, then × 4, up to 12 attempts). The shift changes only the preconditioner, not the system that is solved.
tetfea.js
lines 318 to 347
PartCraft 0.8.1
/** IC(0) of a lower-CSR SPD matrix (in place on a copy) with a diagonal shift raised on breakdown → { L, shift }. */
export function ic0(model, A) {
const { rp, col, N } = model, diagIdx = new Int32Array(N);
for (let r = 0; r < N; r++) diagIdx[r] = rp[r + 1] - 1;
const w = new Float64Array(N);
for (let shift = 0, attempt = 0; attempt < 12; attempt++, shift = shift ? shift * 4 : 1e-3) {
const L = Float64Array.from(A);
if (shift) for (let r = 0; r < N; r++) L[diagIdx[r]] *= 1 + shift;
let ok = true;
for (let i = 0; i < N && ok; i++) {
const s0 = rp[i], s1 = rp[i + 1] - 1;
for (let k = s0; k < s1; k++) w[col[k]] = L[k];
for (let k = s0; k < s1; k++) {
const c = col[k];
// L_ic = (A_ic − Σ_{j<c} L_ij L_cj) / L_cc, row c of L holds L_cj
let s = w[c];
for (let q = rp[c], qe = rp[c + 1] - 1; q < qe; q++) s -= w[col[q]] * L[q];
s /= L[diagIdx[c]];
w[c] = s; L[k] = s;
}
let d = L[s1];
for (let k = s0; k < s1; k++) d -= L[k] * L[k];
for (let k = s0; k < s1; k++) w[col[k]] = 0;
if (!(d > 1e-14 * Math.abs(A[s1]))) { ok = false; break; }
L[s1] = Math.sqrt(d);
}
if (ok) return { L, shift };
}
throw new Error(t('study.tet.precond_failed'));
}
Conjugate-gradient solve for displacements
64-bit
Solves the constrained system for the displacement of every node.
stop when ‖r‖ ÷ ‖P f‖ < tol
The residual is updated recursively (r −= a·q) and not recomputed from K̃u, and the stop test uses that recursive residual relative to the norm of the projected load. Default tolerance 1e-9; the cap is max(1000, 40 √N) iterations. The loop also leaves when p·K̃p is not positive.
tetfea.js
lines 388 to 419
PartCraft 0.8.1
/** Preconditioned CG on K̃ u = P f → { u, iterations, residual, converged }. */
export function tetPCG(model, sys, f, { tol = 1e-9, maxIter = 0, x0 = null, onProgress = null, cancelled = null, label = 'Solving' } = {}) {
const N = model.N, u = x0 ? Float64Array.from(x0) : new Float64Array(N);
const b = sys.pj.proj(Float64Array.from(f));
sys.pj.proj(u);
let bn = 0; for (let i = 0; i < N; i++) bn += b[i] * b[i]; bn = Math.sqrt(bn);
if (bn === 0) return { u: new Float64Array(N), iterations: 0, residual: 0, converged: true };
const r = new Float64Array(N), z = new Float64Array(N), p = new Float64Array(N), q = new Float64Array(N);
sys.apply(u, q); for (let i = 0; i < N; i++) r[i] = b[i] - q[i];
sys.pj.proj(r);
sys.prec(r, z); p.set(z);
let rz = 0; for (let i = 0; i < N; i++) rz += r[i] * z[i];
const cap = maxIter || Math.max(1000, 4 * Math.sqrt(N) * 10);
let it = 0, res = 1, t0 = Date.now();
for (; it < cap; it++) {
let rn = 0; for (let i = 0; i < N; i++) rn += r[i] * r[i];
res = Math.sqrt(rn) / bn;
if (res < tol) break;
if (onProgress && Date.now() - t0 > 200) { t0 = Date.now(); onProgress(Math.min(0.99, Math.log(res) / Math.log(tol)), `${label}: iteration ${it}, residual ${res.toExponential(1)}`); } // i18n-ignore: progress label, never shown (the workers pass on only the fraction)
if (cancelled && (it & 15) === 0 && cancelled()) throw new Error('Cancelled');
sys.apply(p, q); sys.pj.proj(q);
let pq = 0; for (let i = 0; i < N; i++) pq += p[i] * q[i];
if (!(pq > 0)) break;
const a = rz / pq;
for (let i = 0; i < N; i++) { u[i] += a * p[i]; r[i] -= a * q[i]; }
sys.prec(r, z);
let rz2 = 0; for (let i = 0; i < N; i++) rz2 += r[i] * z[i];
const beta = rz2 / rz; rz = rz2;
for (let i = 0; i < N; i++) p[i] = z[i] + beta * p[i];
}
return { u, iterations: it, residual: res, converged: res < tol };
}
Loads on faces and gravity
64-bit
Turns forces, pressures and gravity into nodal loads and totals them.
force: t = total ÷ A; pressure: t = −p n; f = Σ N t dA; gravity: f = ρ g ∫N dV
A is the meshed area of the picked faces from the same 6-point rule, so a force's traction is the total divided by that area. Each rule point's normal is normalized with Math.hypot, with 1 used when the length is zero. Gravity takes g from the load or the constant G_MM imported from materials.js.
tetfea.js
lines 499 to 521
PartCraft 0.8.1
const total = [0, 0, 0];
for (const l of bcs.loads || []) {
if (l.type === 'gravity') {
const d = l.dir || [0, -1, 0], dl = Math.hypot(...d) || 1, g = l.g ?? G_MM;
for (let a = 0; a < model.n; a++) for (let c = 0; c < 3; c++) { const v = density * g * nodeVolN[a] * d[c] / dl; f[3 * a + c] += v; total[c] += v; }
groups.loads.push({ type: 'gravity', nodes: new Int32Array(0) });
continue;
}
const tris = l.nodes ? [] : trisOf(l.faces);
if (!tris.length) throw new Error(t('study.tet.load_no_faces', { type: typeWord(l.type) }));
const A = meshArea(tris), nodes = new Set();
let traction; // (N, n, p) → vector per unit area at the point
if (l.type === 'force') {
const vec = Array.isArray(l.dir) ? l.dir.map((v) => v / (Math.hypot(...l.dir) || 1)) : null;
traction = (n) => { const ln = Math.hypot(...n) || 1; return vec ? vec.map((v) => v * l.total / A) : n.map((v) => -v / ln * l.total / A); };
} else {
traction = (n) => { const ln = Math.hypot(...n) || 1; return n.map((v) => -v / ln * l.p); };
}
for (const t of tris) faceIntegrate(model, t, (Nn, n, p, v) => {
const tr = traction(n), da = Math.hypot(...n);
for (let k = 0; k < 6; k++) { nodes.add(v[k]); for (let c = 0; c < 3; c++) { const val = Nn[k] * tr[c] * da; f[3 * v[k] + c] += val; total[c] += val; } }
});
groups.loads.push({ type: l.type, nodes: Int32Array.from(nodes), area: A, exactArea: exactArea && l.faces ? l.faces.reduce((s, fi) => s + (exactArea[fi] || 0), 0) : null });
Support directions on faces
64-bit
Finds which nodes of a supported face are held and the normal direction used for a sliding support.
held when distance to the exact face ≤ 0.002 × mean edge length
The normal comes from the exact plane, cylinder or surface record where the face has one, with its sign taken from the element's own normal; otherwise from the element. The on-face test is relative to the triangle's mean edge length (2e-3).
tetfea.js
lines 470 to 486
PartCraft 0.8.1
let hl = 0; for (let k = 0; k < 3; k++) { const a = v[k], b = v[(k + 1) % 3]; hl += Math.hypot(model.X[3 * a] - model.X[3 * b], model.X[3 * a + 1] - model.X[3 * b + 1], model.X[3 * a + 2] - model.X[3 * b + 2]) / 3; }
// the normal at each of the 6 nodes (corners at L = vertices, midsides at edge middles); a node of a triangle
// that straddles a sharp edge (the mesher's chamfer) is held only if it lies on the face itself
[[0, 0], [1, 0], [0, 1], [0.5, 0], [0.5, 0.5], [0, 0.5]].forEach(([L1, L2], k) => {
const p = [model.X[3 * v[k]], model.X[3 * v[k] + 1], model.X[3 * v[k] + 2]], dist = exactDist(p, fi);
if (dist !== null && dist > 2e-3 * hl) return;
const { d1, d2 } = t6(L1, L2), a = [0, 0, 0], b = [0, 0, 0];
for (let q = 0; q < 6; q++) for (let cc = 0; cc < 3; cc++) { const x = model.X[3 * v[q] + cc]; a[cc] += d1[q] * x; b[cc] += d2[q] * x; }
let n = [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
const r = faces && faces[fi], nl = Math.hypot(...n) || 1;
// exact normals where the face has one (sign from the element's outward side)
if (r && r.type === 'plane' && r.n) { const s = Math.sign(n[0] * r.n[0] + n[1] * r.n[1] + n[2] * r.n[2]) || 1; n = r.n.map((x) => s * x * nl); }
else if (r && r.type === 'cylinder' && r.axis) { const o = r.axis.o, ad = r.axis.d, q = [p[0] - o[0], p[1] - o[1], p[2] - o[2]], t = q[0] * ad[0] + q[1] * ad[1] + q[2] * ad[2], rad = [q[0] - t * ad[0], q[1] - t * ad[1], q[2] - t * ad[2]], rl = Math.hypot(...rad) || 1, s = Math.sign(n[0] * rad[0] + n[1] * rad[1] + n[2] * rad[2]) || 1; n = rad.map((x) => s * x / rl * nl); }
else if (r && r.surf) { const sn = closestOnSurface(r.surf, p).n, s = Math.sign(n[0] * sn[0] + n[1] * sn[1] + n[2] * sn[2]) || 1; n = sn.map((x) => s * x * nl); }
let m = normals.get(v[k]); if (!m) normals.set(v[k], (m = new Map()));
const sm = m.get(fi) || [0, 0, 0]; m.set(fi, [sm[0] + n[0], sm[1] + n[1], sm[2] + n[2]]);
});
Is the part held
64-bit
Checks that the supports stop all six rigid-body motions.
free when an eigenvalue of Σ rᵀr is ≤ 10⁻⁹ × the largest, or ≤ 10⁻¹²
Node positions are centered on the bounding box and divided by its largest side before the 6 × 6 matrix is built, so the test does not depend on the part's size or position. The second threshold (1e-12) is absolute on that scaled matrix.
tetfea.js
lines 526 to 539
PartCraft 0.8.1
/** Is the body held? (fea.js rigidCheck on node coordinates) → null or the free motions in words. */
export function tetRigidCheck(model, cons) {
const X = model.X, lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
for (let i = 0; i < model.n; i++) for (let a = 0; a < 3; a++) { lo[a] = Math.min(lo[a], X[3 * i + a]); hi[a] = Math.max(hi[a], X[3 * i + a]); }
const c = lo.map((v, a) => (v + hi[a]) / 2), L = Math.max(...hi.map((v, a) => v - lo[a])) || 1, G = new Float64Array(36);
for (let m = 0; m < cons.nodes.length; m++) {
const nd = cons.nodes[m], p = [0, 1, 2].map((a) => (X[3 * nd + a] - 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]];
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), 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')], 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;
}
Natural frequencies
64-bit
The lowest natural frequencies and mode shapes of the supported part.
K̃ x = λ M x; f = √max(λ, 0) ÷ 2π
Block iteration with Rayleigh–Ritz in 64-bit; it stops when every wanted mode's residual ‖K̃x − λMx‖ ÷ (|λ| ‖Mx‖) is below 1e-4, or after 3000 iterations. Directions of the Gram matrix below 1e-10 × its largest eigenvalue are dropped. A negative eigenvalue is set to 0 before the square root.
tetfea.js
lines 641 to 699
PartCraft 0.8.1
/**
* The lowest `modes` eigenpairs of K̃ x = λ M x on the admissible space (LOBPCG, block size modes + 2, IC(0)
* preconditioner) → { freqs (Hz), shapes, iterations, residuals, converged }.
*/
export function tetModal(model, sys, Mv, { modes = 4, tol = 1e-4, maxIter = 3000, onProgress = null, cancelled = null } = {}) {
const N = model.N, q = modes + Math.max(3, Math.ceil(modes / 2)), proj = sys.pj.proj;
const Kx = (x) => proj(sys.apply(x, new Float64Array(N))), Mx = (x) => massMul(model, Mv, x, new Float64Array(N));
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 ? (i % 3 === 1 ? 1 : 0.1) : rand(); return proj(v); });
// the preconditioned start: a few smoothing steps push the random vectors towards the low modes
X = X.map((x) => { const z = new Float64Array(N); sys.prec(Mx(x), z); return z; });
const dot = (a, b) => { let s = 0; for (let i = 0; i < N; i++) s += a[i] * b[i]; return s; };
const comb = (vs, C, cols, ld) => cols.map((j) => { const out = new Float64Array(N); vs.forEach((v, i) => { const c = C[i * ld + j]; if (c) for (let k = 0; k < N; k++) out[k] += c * v[k]; }); return out; });
// Rayleigh–Ritz on basis S (with K S, M S) → { vals, C (ld = S.length), cols kept }
const rr = (S, KS, MS) => {
const s = S.length, GA2 = new Float64Array(s * s), GB2 = new Float64Array(s * s);
for (let i = 0; i < s; i++) for (let j = i; j < s; j++) { GA2[i * s + j] = GA2[j * s + i] = dot(S[i], KS[j]); GB2[i * s + j] = GB2[j * s + i] = dot(S[i], MS[j]); }
// GB = V Λ Vᵀ; drop near-dependent directions; T = V Λ^-½
const { values: bv, vectors: bV } = symEig(GB2, s), bmax = Math.max(...bv);
const keep = bv.map((v, i) => i).filter((i) => bv[i] > 1e-10 * bmax);
const k = keep.length, T = new Float64Array(s * k);
keep.forEach((c, j) => { const sc = 1 / Math.sqrt(bv[c]); for (let i = 0; i < s; i++) T[i * k + j] = bV[i * s + c] * sc; });
const Cm = new Float64Array(k * k);
for (let a = 0; a < k; a++) for (let b = a; b < k; b++) { let t = 0; for (let i = 0; i < s; i++) { let u2 = 0; for (let j = 0; j < s; j++) u2 += GA2[i * s + j] * T[j * k + b]; t += T[i * k + a] * u2; } Cm[a * k + b] = Cm[b * k + a] = t; }
const { values, vectors } = symEig(Cm, k), order = values.map((v, i) => i).sort((a, b) => values[a] - values[b]);
const C = new Float64Array(s * k);
order.forEach((c, j) => { for (let i = 0; i < s; i++) { let t = 0; for (let l = 0; l < k; l++) t += T[i * k + l] * vectors[l * k + c]; C[i * k + j] = t; } });
return { vals: order.map((c) => values[c]), C, k };
};
let KX = X.map(Kx), MX = X.map(Mx);
let r0 = rr(X, KX, MX), lam = r0.vals.slice(0, q);
X = comb(X, r0.C, [...Array(q).keys()], r0.k); KX = X.map(Kx); MX = X.map(Mx);
let P = null, KP = null, MP = null, it = 0, res = [];
for (; it < maxIter; it++) {
if (cancelled && cancelled()) throw new Error('Cancelled');
const R = X.map((x, j) => { const r = new Float64Array(N); for (let i = 0; i < N; i++) r[i] = KX[j][i] - lam[j] * MX[j][i]; return proj(r); });
res = R.map((r, j) => Math.sqrt(dot(r, r)) / (Math.abs(lam[j]) * Math.sqrt(dot(MX[j], MX[j])) || 1));
if (onProgress) onProgress(Math.min(0.99, (it + 1) / maxIter + 0.5 * Math.max(0, Math.log(Math.max(...res.slice(0, modes))) / Math.log(tol))), `Modes: iteration ${it + 1}, residual ${Math.max(...res.slice(0, modes)).toExponential(1)}`); // i18n-ignore: progress label, never shown
if (res.slice(0, modes).every((x) => x < tol)) break;
const W = R.map((r) => { const z = new Float64Array(N); sys.prec(r, z); return z; });
// M-normalise W (the Gram matrix drops near-dependent directions by size, so every block must be on one scale)
const MW = W.map(Mx);
W.forEach((w, j) => { const l = Math.sqrt(Math.max(dot(w, MW[j]), 1e-300)); for (let i = 0; i < N; i++) { w[i] /= l; MW[j][i] /= l; } });
const KW = W.map(Kx);
const S = [...X, ...W, ...(P || [])], KS = [...KX, ...KW, ...(KP || [])], MS = [...MX, ...MW, ...(MP || [])];
const r = rr(S, KS, MS);
if (r.k < q) { P = null; KP = null; MP = null; continue; }
const cols = [...Array(q).keys()];
lam = r.vals.slice(0, q);
// P = the W and P parts of the new X
const Cp = Float64Array.from(r.C); for (let i = 0; i < q; i++) for (let j = 0; j < r.k; j++) Cp[i * r.k + j] = 0;
P = comb(S, Cp, cols, r.k); KP = comb(KS, Cp, cols, r.k); MP = comb(MS, Cp, cols, r.k);
X = comb(S, r.C, cols, r.k); KX = comb(KS, r.C, cols, r.k); MX = comb(MS, r.C, cols, r.k);
// P normalised for conditioning
for (let j = 0; j < q; j++) { const l = Math.sqrt(Math.max(dot(P[j], MP[j]), 1e-300)); for (let i = 0; i < N; i++) { P[j][i] /= l; KP[j][i] /= l; MP[j][i] /= l; } }
}
const converged = res.slice(0, modes).every((x) => x < tol);
return { freqs: lam.slice(0, modes).map((l) => Math.sqrt(Math.max(l, 0)) / (2 * Math.PI)), shapes: X.slice(0, modes), iterations: it, residuals: res.slice(0, modes), converged };
}
The temperature rise at every node over the stress-free temperature, and the load it adds.
ΔT = T − T_ref (20 when not given); f += thermal load
The rise is a single subtraction per node into a Float64Array, from a uniform value, given temperatures, a steady heat solve or one step of a transient run. The thermal load vector is added to the mechanical loads entry by entry.
tet-elements.js
lines 226 to 237
PartCraft 0.8.1
const dT = new Float64Array(model.n);
if (tl.mode === 'uniform') dT.fill(tl.dT);
else if (job.temperatures) for (let n = 0; n < model.n; n++) dT[n] = job.temperatures[n] - (tl.Tref ?? 20); // given nodal temperatures (tests)
else if (tl.mode === 'transient') { // round 5 (builder C): the transient run's temperatures at the step nearest tl.time (buckling over time)
if (!job.transient) throw new Error(t('study.thermal.history.needs_transient'));
onProgress(0.35, 'Temperatures over time'); // i18n-ignore: progress label, never shown
at = stepAt(job.transient, tl.time);
heat = heatOnMesh(job, mesh, { frameSteps: [at.step], cancelled });
for (let n = 0; n < model.n; n++) dT[n] = heat.res.frames[0][n] - (tl.Tref ?? 20);
} else { onProgress(0.35, 'Temperatures'); heat = heatOnTets(job, model, null, cancelled, { steadyOnly: true }); for (let n = 0; n < model.n; n++) dT[n] = heat.res.T[n] - (tl.Tref ?? 20); }
const ft = thermalLoadVector(mesh, mat, dT);
for (let i = 0; i < ft.length; i++) bc.f[i] += ft[i];
Buckling pre-stress and geometric stiffness
64-bit
Solves the pre-stress and assembles the stress stiffness that the buckling eigen-solve uses.
K u = f (tol 10⁻¹⁰); G = −Σₑ Gₑ; each inner solve tol 10⁻⁸
The pre-stress solve uses a relative residual of 1e-10 and each solve inside the eigen-iteration 1e-8; the worst unconverged inner residual is kept and returned. The load factors themselves are computed by bucklingEigen in buckling.js.
tet-elements.js
lines 267 to 284
PartCraft 0.8.1
const sol = tetPCG(model, sys, bc.f, { tol: 1e-10, cancelled });
if (!sol.converged) throw new Error(t('study.buckling.prestress_no_converge', { res: sol.residual.toExponential(1), n: sol.iterations }));
const sig = stressAtPoints(mesh, sol.u, mat, s.thermal);
// Kσ assembled once on the node graph (the element's npe×npe block ⊗ I₃), applied as −Kσ x
const Gv = new Float64Array(model.nodeAdj.length), X = new Float64Array(30), Ge = new Float64Array(100);
for (let e = 0; e < mesh.ne; e++) {
mesh.elemX(e, X); TET10.geometricStiffness(X, sig.subarray(e * 24, e * 24 + 24), Ge);
for (let a = 0; a < 10; a++) for (let b = 0; b < 10; b++) { const A = model.tets10[10 * e + a], B = model.tets10[10 * e + b]; if (B <= A) Gv[nodeEntry(model, A, B)] -= Ge[a * 10 + b]; }
}
const applyG = (x, y) => massMul(model, Gv, x, y);
void geometricOperator;
let solves = 0, worst = 0;
const eig = bucklingEigen({
n: model.N, modes: job.modes || 4, cancelled, onProgress: (f, t) => onProgress(0.5 + 0.5 * f, t),
applyK: (x, y) => sys.pj.proj(sys.apply(x, y)),
applyG: (x, y) => sys.pj.proj(applyG(x, y)),
solve: (b, x0) => { const r = tetPCG(model, sys, b, { x0, tol: 1e-8, cancelled }); solves++; if (!r.converged) worst = Math.max(worst, r.residual); return r.u; },
});
In detail
Method
Constraints as per-node projections P, solved as K̃ = P K P + α (I − P) with α the mean diagonal of K. Preconditioner: incomplete Cholesky IC(0) up to 5,000 unknowns, the multigrid of amg.js above that. Preconditioned conjugate gradients for displacements; LOBPCG for modes (the file header cites Knyazev 2001). Face loads by the 6-point triangle rule; gravity from ρ g ∫N dV.
Number types
64-bit numbers and Float64Array for matrices, vectors and loads; Int32Array for the sparse pattern.
Tolerances
Conjugate gradients: ‖r‖ ÷ ‖P f‖ below 1e-9 by default, 1e-8 or the job's value for a static study, 1e-10 for the buckling pre-stress, 1e-8 inside the buckling eigen-solve; at most max(1000, 40 √N) iterations. IC(0): a pivot must exceed 1e-14 × the matrix diagonal, else the diagonal is scaled by 1 + shift, shift 1e-3 then × 4, up to 12 attempts. Modes: relative residual below 1e-4, at most 3000 iterations; Gram directions below 1e-10 × the largest eigenvalue are dropped. Held check: a support eigenvalue counts as free when it is not above 1e-9 × the largest or not above 1e-12. A node is held by a face support when it lies within 2e-3 × the triangle's mean edge length of the face's exact surface.
Limits
Gravity uses G_MM from materials.js when the load gives no value. A force on faces is spread as the total divided by the meshed area. The conjugate-gradient loop also stops when p·K̃p is not positive; the caller then reports no convergence.
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.