Finds the lowest load factors at which a part buckles and their shapes, and sets up thermal loads for a stress study on the voxel model.
How the precision is kept
Everything is Float64Array. The buckling factors come from an iteration that stops when each wanted value changes by less than 1e-6 (relative) between two rounds, with inner linear solves stopped at 1e-9 and the pre-stress solve at 1e-10. The model is a voxel grid of the part, so the grid spacing and the iteration tolerances set the accuracy, not the arithmetic.
Calculations under this policy
7 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.
The six stress components at every integration point of every element, from the displacements.
σ = D (ε − α ΔT I) per point, by the element's own routine
Gathers each element's nodal displacements into Float64Array work arrays and calls the element's stress routine; results are stored in a Float64Array. The arithmetic of the stress itself is in the element code, not in this file.
buckling.js
lines 25 to 36
PartCraft 0.8.1
/** σ at every element's integration points (6·nip per element) from a displacement field (3 per node). th: null | { alpha, dT (nn) } */
export function stressAtPoints(mesh, u, { E, nu }, th = null, scale = null) {
const { ne, npe, element: el } = mesh, nip = el.nip;
const out = new Float64Array(ne * nip * 6), X = new Float64Array(3 * npe), ue = new Float64Array(3 * npe), te = new Float64Array(npe), sg = new Float64Array(6 * nip);
for (let e = 0; e < ne; e++) {
mesh.elemX(e, X);
for (let a = 0; a < npe; a++) { const n = mesh.conn[e * npe + a]; ue[3 * a] = u[3 * n]; ue[3 * a + 1] = u[3 * n + 1]; ue[3 * a + 2] = u[3 * n + 2]; if (th) te[a] = th.dT[n]; }
el.stressPoints(X, E * (scale ? scale[e] : 1), nu, ue, sg, th ? { alpha: th.alpha, dT: te } : null);
out.set(sg, e * nip * 6);
}
return out;
}
Geometric stiffness times a vector
64-bit
Applies the stress-dependent stiffness Kσ to a vector, element by element.
y = s × Σ_elements G_e x_e on each of the three directions
Each element's contribution is summed over its nodes in order and added into the result; contributions from neighboring elements accumulate in element order. No 32-bit storage.
buckling.js
lines 38 to 56
PartCraft 0.8.1
/** y = s · Kσ x (3 dofs per node), element by element; s = −1 gives the −Kσ of the eigenproblem. */
export function geometricOperator(mesh, sig, s = 1) {
const { ne, npe, element: el } = mesh, nip = el.nip;
const X = new Float64Array(3 * npe), G = new Float64Array(npe * npe);
return (x, y) => {
y.fill(0);
for (let e = 0; e < ne; e++) {
mesh.elemX(e, X);
el.geometricStiffness(X, sig.subarray(e * nip * 6, (e + 1) * nip * 6), G);
for (let a = 0; a < npe; a++) {
const na = 3 * mesh.conn[e * npe + a];
let t0 = 0, t1 = 0, t2 = 0;
for (let b = 0; b < npe; b++) { const g = G[a * npe + b], nb = 3 * mesh.conn[e * npe + b]; t0 += g * x[nb]; t1 += g * x[nb + 1]; t2 += g * x[nb + 2]; }
y[na] += s * t0; y[na + 1] += s * t1; y[na + 2] += s * t2;
}
}
return y;
};
}
Buckling factors by subspace iteration
64-bit
The lowest positive buckling load factors and their mode shapes.
Z = K⁻¹ G X; reduce (Zᵀ G Z, Zᵀ K Z); Kr = L Lᵀ; eig(L⁻¹ Gr L⁻ᵀ) = μ; λ = 1 ÷ μ
Stops when every wanted μ changes by at most 1e-6 × |μ| from the previous round, or after 80 rounds; the result carries a flag named converged. The projected matrices are plain dot products summed in order over all unknowns. A pivot of the small Cholesky factor is floored at 1e-300, so a dependent trial vector does not stop the run. The start vectors come from a fixed-seed generator, so runs repeat.
buckling.js
lines 65 to 111
PartCraft 0.8.1
/**
* The largest μ of K⁻¹ G φ = μ φ (G = −Kσ) by subspace iteration with Rayleigh–Ritz → { factors λ = 1/μ (lowest positive first),
* shapes, iterations, converged, negative (the smallest reverse-load factor |λ| found, or null) }.
* op { n, solve(b, x0) → x (admissible, K⁻¹ b), applyK(x, y), applyG(x, y), modes, maxIter 80, tol 1e-6, onProgress, cancelled }
*/
export function bucklingEigen({ n, solve, applyK, applyG, modes = 4, maxIter = 80, tol = 1e-6, onProgress = null, cancelled = null }) {
const q = Math.max(modes + 4, Math.min(2 * modes, modes + 8));
let rnd = 4321; const rand = () => ((rnd = (rnd * 1103515245 + 12345) % 2147483648) / 2147483648) - 0.5;
let X = Array.from({ length: q }, () => Float64Array.from({ length: n }, rand)), Y = X.map(() => null);
let mu = null, prev = null, it = 0, converged = false;
const g = new Float64Array(n);
for (; it < maxIter; it++) {
if (cancelled && cancelled()) throw new Error('Cancelled');
const Z = X.map((x, k) => { applyG(x, g); return solve(g, Y[k]); });
const KZ = Z.map((z) => applyK(z, new Float64Array(n))), GZ = Z.map((z) => applyG(z, new Float64Array(n)));
const Kr = new Float64Array(q * q), Gr = new Float64Array(q * q);
for (let a = 0; a < q; a++) for (let b = a; b < q; b++) {
let sk = 0, sg = 0; for (let i = 0; i < n; i++) { sk += Z[a][i] * KZ[b][i]; sg += Z[a][i] * GZ[b][i]; }
Kr[a * q + b] = Kr[b * q + a] = sk; Gr[a * q + b] = Gr[b * q + a] = sg;
}
// Kr = L Lᵀ, C = L⁻¹ Gr L⁻ᵀ, eig C, y = L⁻ᵀ w
const L = new Float64Array(q * q);
for (let i = 0; i < q; i++) for (let j = 0; j <= i; j++) {
let s = Kr[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 = FE.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] * Gr[a * q + b] * Li[j * q + b]; C[i * q + j] = s; }
const { values, vectors } = FE.symEig(C, q);
const order = values.map((v, i) => i).sort((a, b) => values[b] - values[a]); // largest μ (smallest positive λ) first
mu = order.map((i) => values[i]);
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 * mu[k])); // warm start: K⁻¹G x ≈ μ x
if (onProgress) onProgress((it + 1) / maxIter, `Buckling: iteration ${it + 1}`); // i18n-ignore: progress label, never shown
const want = mu.slice(0, modes);
if (prev && want.every((m, i) => Math.abs(m - prev[i]) <= tol * Math.abs(m))) { converged = true; it++; break; }
prev = want.slice();
}
const pos = mu.map((m, i) => ({ m, i })).filter((p) => p.m > 0).slice(0, modes);
const negs = mu.filter((m) => m < 0);
return { factors: pos.map((p) => 1 / p.m), shapes: pos.map((p) => X[p.i]), iterations: it, converged, negative: negs.length ? 1 / Math.min(...negs) : null };
}
Solver tolerances of a buckling study
64-bit
How tightly the pre-stress solve and each inner solve are converged.
pre-stress: 1e-10; each K⁻¹ solve: 1e-9
An inner solve that does not converge does not stop the study; its residual is recorded and the largest is returned as solveResidual.
The temperature change that drives thermal stress: a uniform rise, or the temperature field minus the stress-free temperature.
ΔT = T − T_ref, T_ref = 20 when not given
One subtraction per node in a Float64Array. When the study gives no reference temperature the code uses 20.
buckling.js
lines 139 to 153
PartCraft 0.8.1
const dT = new Float64Array(model.nn);
if (tl.mode === 'uniform') { for (let n = 0; n < model.nn; n++) if (model.active[n]) dT[n] = tl.dT; }
else if (job.temperatures) { for (let n = 0; n < model.nn; n++) if (model.active[n]) dT[n] = job.temperatures[n] - (tl.Tref ?? 20); } // given nodal temperatures (tests)
else {
onProgress(0, 'Temperatures');
// 'field': the steady temperatures; 'transient': the run's temperatures at the step nearest tl.time
if (tl.mode === 'transient' && !job.transient) throw new Error(t('study.thermal.history.needs_transient'));
at = tl.mode === 'transient' ? stepAt(job.transient, tl.time) : null;
heat = TH.heatOnGrid(job, grid, null, cancelled, at ? { frameSteps: [at.step] } : { steadyOnly: true });
const T = at ? heat.res.frames[0] : heat.res.T;
for (let n = 0; n < model.nn; n++) if (model.active[n]) dT[n] = T[n] - (tl.Tref ?? 20);
}
const ft = TH.thermalLoadVector(mesh, mat, dT);
for (let i = 0; i < ft.length; i++) bc.f[i] += ft[i];
thermal = { dT, alpha: mat.alpha, corr: TH.stressCorrection(type, mat.nu), at };
Time step for a transient thermal load
64-bit
Which saved time step of a transient run is nearest the requested time.
step = round(time ÷ (t ÷ steps)), held in [0, steps]; time = step × dt
The requested time is snapped to the nearest step, so the temperatures used can be up to half a step away from the time asked for; the snapped time is returned.
buckling.js
lines 115 to 119
PartCraft 0.8.1
/** The time step a transient thermal load uses: the step nearest `time` (default: the end) → { step, time }. */
export function stepAt(tr, time) {
const dt = tr.t / tr.steps, step = Math.max(0, Math.min(tr.steps, Math.round((time ?? tr.t) / dt)));
return { step, time: step * dt };
}
Voxel model volume and mass
64-bit
The volume and mass of the analysis model listed with the study.
volume = elements × h³; mass = volume × density (tonne/mm³) × 1000
These are the volume and mass of the voxel grid, a count of cubes, not of the exact part. The factor 1000 converts tonnes to kilograms and is exact.
Linear buckling as K⁻¹(−Kσ)φ = μφ, λ = 1 ÷ μ, by subspace iteration with Rayleigh–Ritz on the projected pair, as described in the file's header. Stresses at the element integration points give the geometric stiffness Kσ.
Number types
64-bit numbers and Float64Array throughout.
Tolerances
Eigenvalue change ≤ 1e-6 × |μ| between iterations, at most 80 iterations. Inner conjugate-gradient solves to 1e-9, pre-stress to 1e-10, thermal static solve to the study's tolerance or 1e-8. Projected Cholesky pivots are floored at 1e-300.
Checked by
tools/partcraft-test/test-thermal.mjs calls solveThermoStatic and solveBuckling on bars and columns. It loads the development build's copy of this file.
Limits
The header states the limits: linear, no imperfections, elastic; real structures usually buckle below these loads. The stopping test compares successive iterates and is not an error bound. Reported volume and mass are those of the voxel model.
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.