Solves an elastic–plastic stress study in load steps, with contact between pieces when asked, and reports the peak stress, plastic strain, displacement, safety factor and reactions.
How the precision is kept
Strains, stresses, the tangent and the Newton iteration are 64-bit (Float64Array throughout nonlinear.js). The accuracy of the solution is set by the mesh, the load steps and the stopping tolerances. The per-element von Mises stress and plastic strain are then stored in Float32Array (elementMaxima), and the reported maximum stress, maximum plastic strain and minimum safety factor are read from those 32-bit values. The maximum displacement is taken from a Float64Array.
Open item. Return the per-element von Mises stress and plastic strain in a Float64Array, so the reported maxima and the safety factor are not read from 32-bit values.
Calculations under this policy
11 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.
nonlinear.jsPartCraft 0.8.1 32-bit in part5 calculationsReviewed
Tests tools/partcraft-test/test-sim08.mjs
Bilinear material constants
64-bit
Shear modulus, bulk modulus and hardening modulus from E, ν, yield stress and tangent modulus.
G = E ÷ (2(1 + ν)); K = E ÷ (3(1 − 2ν)); H = E × E_t ÷ (E − E_t)
ν must be below 0.5 and E_t below E, so neither division is by zero; close to those limits 1 − 2ν and E − E_t subtract nearly equal numbers and K or H become large.
nonlinear.js
lines 23 to 29
PartCraft 0.8.1
/** A bilinear (or linear, when sy is missing) material from { E, nu, yield, Et }: the numbers the return map needs. */
export function bilinearMaterial({ E, nu, yield: sy = Infinity, Et = 0, density = 0 }) {
if (!(E > 0) || !(nu >= 0 && nu < 0.5)) throw new Error('bilinear material: E > 0 and 0 ≤ ν < 0.5 are required'); // i18n-ignore: internal check, the UI validates first
if (!(Et >= 0 && Et < E)) throw new Error('bilinear material: 0 ≤ Et < E is required'); // i18n-ignore: internal check
const G = E / (2 * (1 + nu)), K = E / (3 * (1 - 2 * nu));
return { E, nu, sy, Et, H: E * Et / (E - Et), G, K, density, D: elasticD(nu).map((v) => v * E) };
}
Stress update at a Gauss point
64-bit
The stress, the new plastic strain and the algorithmic tangent for a given total strain.
q = √(1.5 s:s); f = q − (σ_y + H α); Δγ = f ÷ (3G + H); σ = p 1 + (1 − 3GΔγ ÷ q) s
The yield test is f > 0 with no tolerance. The division by q happens only on the plastic branch, where q exceeds the current yield stress. The deviatoric norm is formed from squares, so nothing cancels there; the trial stress is a 6 × 6 product in 64-bit numbers.
nonlinear.js
lines 33 to 67
PartCraft 0.8.1
/**
* Stress update at one Gauss point. eps: total strain (6, engineering shear); st: { ep: Float64Array(6) (engineering shear),
* a } is read and the new state written to out.ep / out.a. Returns the stress in sig (6) and the algorithmic tangent in Dep
* (36, row-major). Pure function of its arguments.
*/
export function returnMap(mat, eps, ep0, a0, sig, Dep, out) {
const { G, K, H, sy, D } = mat;
const ee = [0, 0, 0, 0, 0, 0];
for (let r = 0; r < 6; r++) ee[r] = eps[r] - ep0[r];
const st = [0, 0, 0, 0, 0, 0];
for (let r = 0; r < 6; r++) { let s = 0; for (let k = 0; k < 6; k++) s += D[r * 6 + k] * ee[k]; st[r] = s; }
const p = (st[0] + st[1] + st[2]) / 3, s = [st[0] - p, st[1] - p, st[2] - p, st[3], st[4], st[5]];
const s2 = s[0] * s[0] + s[1] * s[1] + s[2] * s[2] + 2 * (s[3] * s[3] + s[4] * s[4] + s[5] * s[5]);
const q = Math.sqrt(1.5 * s2), f = q - (sy + H * a0);
for (let r = 0; r < 6; r++) out.ep[r] = ep0[r];
out.a = a0;
if (!(f > 0)) { // elastic
for (let r = 0; r < 6; r++) sig[r] = st[r];
Dep.set(D);
return false;
}
const dg = f / (3 * G + H), th = 1 - 3 * G * dg / q, nrm = Math.sqrt(s2);
for (let r = 0; r < 6; r++) sig[r] = (r < 3 ? p : 0) + th * s[r];
// plastic strain increment: Δεp = Δγ · 3/2 · s / q (tensor); the engineering shear components are twice the tensor's
for (let r = 0; r < 6; r++) out.ep[r] += dg * 1.5 * s[r] / q * (r < 3 ? 1 : 2);
out.a = a0 + dg;
const thb = 3 * G / (3 * G + H) - (1 - th), nv = s.map((v) => v / nrm);
for (let i = 0; i < 6; i++) for (let j = 0; j < 6; j++) {
let v = K * (i < 3 && j < 3 ? 1 : 0);
const id = (i === j ? (i < 3 ? 1 : 0.5) : 0) - (i < 3 && j < 3 ? 1 / 3 : 0); // I_dev on engineering strains
v += 2 * G * th * id - 2 * G * thb * nv[i] * nv[j];
Dep[i * 6 + j] = v;
}
return true;
}
Element tangent and internal force
64-bit
Integrates each element's tangent stiffness and internal force over its four Gauss points.
K_e = Σ w Bᵀ D_ep B; f_e = Σ w Bᵀ σ; w = weight × det J
Float64Array work arrays; strains are sums of gradient times displacement over the 10 nodes. An element whose Jacobian determinant is not positive stops the solve with an error.
nonlinear.js
lines 69 to 111
PartCraft 0.8.1
/** Element tangent stiffness Ke (30×30) and internal force fe (30) from the strain at the 4 Gauss points. */
function elementResponse(X, el, ue, mat, ep, al, epN, alN, Ke, fe, sigOut, g, sig, Dep, tmp, DB) {
Ke.fill(0); fe.fill(0);
let plastic = 0;
for (let p = 0; p < 4; p++) {
const det = jac(X, el, SH4[p].dN, g);
if (!(det > 0)) throw new Error('inverted element'); // i18n-ignore: internal; the mesher never gives one
let exx = 0, eyy = 0, ezz = 0, gxy = 0, gyz = 0, gzx = 0;
for (let a = 0; a < 10; a++) {
const dx = g[3 * a], dy = g[3 * a + 1], dz = g[3 * a + 2], ux = ue[3 * a], uy = ue[3 * a + 1], uz = ue[3 * a + 2];
exx += dx * ux; eyy += dy * uy; ezz += dz * uz; gxy += dy * ux + dx * uy; gyz += dz * uy + dy * uz; gzx += dz * ux + dx * uz;
}
const eps = tmp.eps; eps[0] = exx; eps[1] = eyy; eps[2] = ezz; eps[3] = gxy; eps[4] = gyz; eps[5] = gzx;
const o = tmp.out, base = 6 * p;
const pl = returnMap(mat, eps, ep.subarray(base, base + 6), al[p], sig, Dep, o);
if (pl) plastic++;
epN.set(o.ep, base); alN[p] = o.a;
if (sigOut) sigOut.set(sig, base);
const w = GP4[p][3] * det;
// DB = Dep·B, then Ke += w BᵀDB and fe += w Bᵀσ
for (let a = 0; a < 10; a++) {
const dx = g[3 * a], dy = g[3 * a + 1], dz = g[3 * a + 2];
for (let r = 0; r < 6; r++) {
DB[r * 30 + 3 * a] = Dep[r * 6] * dx + Dep[r * 6 + 3] * dy + Dep[r * 6 + 5] * dz;
DB[r * 30 + 3 * a + 1] = Dep[r * 6 + 1] * dy + Dep[r * 6 + 3] * dx + Dep[r * 6 + 4] * dz;
DB[r * 30 + 3 * a + 2] = Dep[r * 6 + 2] * dz + Dep[r * 6 + 4] * dy + Dep[r * 6 + 5] * dx;
}
fe[3 * a] += w * (dx * sig[0] + dy * sig[3] + dz * sig[5]);
fe[3 * a + 1] += w * (dy * sig[1] + dx * sig[3] + dz * sig[4]);
fe[3 * a + 2] += w * (dz * sig[2] + dy * sig[4] + dx * sig[5]);
}
for (let a = 0; a < 10; a++) {
const dx = g[3 * a], dy = g[3 * a + 1], dz = g[3 * a + 2];
for (let j = 0; j < 30; j++) {
const c0 = DB[j], c1 = DB[30 + j], c2 = DB[60 + j], c3 = DB[90 + j], c4 = DB[120 + j], c5 = DB[150 + j];
Ke[(3 * a) * 30 + j] += w * (dx * c0 + dy * c3 + dz * c5);
Ke[(3 * a + 1) * 30 + j] += w * (dy * c1 + dx * c3 + dz * c4);
Ke[(3 * a + 2) * 30 + j] += w * (dz * c2 + dy * c4 + dx * c5);
}
}
}
return plastic;
}
Load-stepped Newton solve
64-bit
Raises the load in steps and iterates each step to equilibrium.
r = λ f − f_int(u) (− contact); converged when |r| ÷ (λ |f|) < tol
The residual is relative to the applied load at that step (tolerance 1e-6 by default, 25 iterations). The linear solves stop at 1e-10, and an unconverged one is still used if its residual is under 1e-6. A step that fails is halved, down to 1/64 of a step, after which the result is returned as not converged at the last good load.
nonlinear.js
lines 142 to 195
PartCraft 0.8.1
/**
* Load-stepped Newton solve. model: makeTetModel; mats/elemMat: materials per element; bc: { f (reference load, N), cons
* (makeConstraints) }; opts { steps (10), tol (1e-6, on the residual relative to the load), maxIter (25), contact, prec,
* onProgress, cancelled, trace (called with { lambda, iteration, residual, contact } each Newton iteration) }.
* → { u, state, steps: [{ lambda, iterations, plastic, residual, contact? }], converged, failedAt? , fint }
*/
export function solveNonlinear(model, mats, elemMat, bc, opts = {}) {
const { steps = 10, tol = 1e-6, maxIter = 25, contact = null, prec = undefined, onProgress = null, cancelled = null, tolPCG = 1e-10, trace = null, maxOuter = 12 } = opts;
const N = model.N, f = bc.f, fn = norm(f) || 1;
let u = new Float64Array(N), state = newState(model), lam = 0, dl = 1 / steps, minDl = 1 / (steps * 64);
const log = [];
const newton = (lamNext, u0) => {
const uT = Float64Array.from(u0); let it = 0, res = Infinity, plastic = 0, ctxInfo = null;
for (; it < maxIter; it++) {
if (cancelled && cancelled()) throw new Error('Cancelled');
const nl = internal(model, mats, elemMat, uT, state, { tangent: true });
const r = new Float64Array(N);
for (let i = 0; i < N; i++) r[i] = lamNext * f[i] - nl.fint[i];
if (contact) ctxInfo = contact.apply(model, uT, r, nl.K); // adds the contact forces to r and the penalty stiffness to K
bc.cons && projectResidual(model, bc.cons, r);
res = norm(r) / (lamNext * fn);
plastic = nl.plastic;
if (trace) trace({ lambda: lamNext, iteration: it, residual: res, contact: ctxInfo });
if (res < tol || (res < tol * 1e3 && it > 0 && contact && contact.relaxed)) return { ok: true, u: uT, it, res, plastic, ctxInfo };
const sys = system(model, nl.K, bc.cons, prec ? { prec } : {});
const sol = tetPCG(model, sys, r, { tol: tolPCG, cancelled });
if (!sol.converged && !(sol.residual < 1e-6)) break;
for (let i = 0; i < N; i++) uT[i] += sol.u[i];
}
return { ok: false, u: uT, it, res, plastic, ctxInfo };
};
while (lam < 1 - 1e-12) {
const lamNext = Math.min(1, lam + dl);
let cur = newton(lamNext, u), outer = 0, total = cur.it;
// contact: re-pair on the converged shape and solve again until the pairing stops moving
while (contact && cur.ok && outer < maxOuter && contact.update(model, cur.u)) { outer++; cur = newton(lamNext, cur.u); total += cur.it; }
if (cur.ok) {
u = cur.u; lam = lamNext;
// commit the state at the converged displacement
const nl = internal(model, mats, elemMat, u, state, { tangent: false });
state = { ep: nl.epN, a: nl.aN };
log.push({ lambda: lam, iterations: total, outer, plastic: cur.plastic, residual: cur.res, contact: cur.ctxInfo });
if (contact) contact.commit(model, u);
if (onProgress) onProgress(lam);
dl = Math.min(1 / steps, dl * 1.5);
} else {
if (contact) contact.update(model, u);
dl /= 2;
if (dl < minDl) return { u, state, steps: log, converged: false, failedAt: lamNext, lambda: lam };
}
}
const fin = internal(model, mats, elemMat, u, state, { tangent: false, stress: true });
return { u, state, steps: log, converged: true, lambda: 1, fint: fin.fint, sig: fin.sig };
}
Von Mises stress and plastic strain per element
32-bit in part
The largest von Mises stress and the largest equivalent plastic strain among each element's Gauss points.
vm_e = max over 4 points of vonMises(σ); α_e = max over 4 points of α
Both results are written into Float32Array. The maximum stress, maximum plastic strain and safety factor the study reports are read from these arrays, so they carry about 7 significant digits.
nonlinear.js
lines 207 to 216
PartCraft 0.8.1
/** Von Mises of each element's Gauss-point stresses (max per element) and the equivalent plastic strain (max per element). */
export function elementMaxima(model, sig, state) {
const m = model.tets10.length / 10, vm = new Float32Array(m), ap = new Float32Array(m);
for (let e = 0; e < m; e++) {
let v = 0, a = 0;
for (let p = 0; p < 4; p++) { v = Math.max(v, vonMises(sig.subarray(24 * e + 6 * p, 24 * e + 6 * p + 6))); a = Math.max(a, state.a[4 * e + p]); }
vm[e] = v; ap[e] = a;
}
return { vm, plasticStrain: ap };
}
nl-worker.jsPartCraft 0.8.1 32-bit in part6 calculationsReviewed
Tests tools/partcraft-test/test-nl-mount.mjs
Reported maximum stress, displacement and safety factor
32-bit in part
The peak von Mises stress and where it is, the peak displacement, and the smallest safety factor against yield.
max σ_vm = largest element value; max |u| = largest nodal √(u_x² + u_y² + u_z²); safety = yield ÷ max σ_vm
The stress maximum is ex.vm[ev], an entry of the Float32Array from elementMaxima, and the safety factor divides the yield stress by it. The displacement maximum comes from a Float64Array and keeps 64-bit precision. With CalculiX results the stress maximum is taken from a Float64Array of nodal values instead.
The total force the supports carry, overall and per support group.
R = Σ over constrained nodes of (f_int − f_applied)
A sum of 64-bit nodal forces in node order. Reactions of opposite sign on different nodes cancel in the total, so a small net reaction can carry the rounding of much larger nodal forces.
nl-worker.js
lines 87 to 88
PartCraft 0.8.1
/** The support reaction on a set of nodes: internal force minus the load applied there (a loaded node on a symmetry plane is in both). */
const sumNodes = (fint, f, nodes) => { const s = [0, 0, 0]; for (const n of nodes) for (let a = 0; a < 3; a++) s[a] += fint[3 * n + a] - f[3 * n + a]; return s; };
Element values averaged to nodes
64-bit
The nodal stress and plastic strain fields the plot colors by.
node value = mean of the values of the elements that share the node
Sums and counts are Float64Array. The averaged field is then copied into a Float32Array for the plot (the helper named at), which is display only; the reported maxima are not taken from it.
nl-worker.js
lines 81 to 86
PartCraft 0.8.1
function nodalFromElements(model, perElem) {
const sum = new Float64Array(model.n), cnt = new Float64Array(model.n), m = model.tets10.length / 10;
for (let e = 0; e < m; e++) for (let a = 0; a < 10; a++) { const n = model.tets10[10 * e + a]; sum[n] += perElem[e]; cnt[n]++; }
for (let i = 0; i < model.n; i++) if (cnt[i]) sum[i] /= cnt[i];
return sum;
}
Separating the pieces of a part
32-bit in part
Splits a part's surface mesh into its connected pieces so each can be meshed on its own for a contact study.
vertices are the same when round(x × 10⁴), round(y × 10⁴), round(z × 10⁴) agree
Vertices are welded on a grid of 1e-4 model units for the connectivity test only; positions are not moved. Each piece's positions are written to a new Float32Array, as the incoming surface mesh is.
nl-worker.js
lines 35 to 54
PartCraft 0.8.1
/** The connected pieces of a surface mesh (vertices welded by position), each as its own mesh; one piece → the mesh itself. */
export function splitMesh(mesh) {
const P = mesh.positions, I = mesh.indices, nv = P.length / 3, weld = new Map(), rep = new Int32Array(nv);
for (let v = 0; v < nv; v++) { const k = Math.round(P[3 * v] * 1e4) + ',' + Math.round(P[3 * v + 1] * 1e4) + ',' + Math.round(P[3 * v + 2] * 1e4); if (!weld.has(k)) weld.set(k, v); rep[v] = weld.get(k); }
const parent = Int32Array.from({ length: nv }, (_, i) => i), find = (a) => { while (parent[a] !== a) { parent[a] = parent[parent[a]]; a = parent[a]; } return a; };
for (let t = 0; t < I.length / 3; t++) { const a = find(rep[I[3 * t]]), b = find(rep[I[3 * t + 1]]), c = find(rep[I[3 * t + 2]]); parent[b] = a; parent[find(c)] = a; }
const root = (v) => find(rep[v]), groups = new Map();
for (let t = 0; t < I.length / 3; t++) { const r = root(I[3 * t]); (groups.get(r) || groups.set(r, []).get(r)).push(t); }
if (groups.size <= 1) return [mesh];
const keyOf = (x, y, z) => Math.round(x * 1e4) + ',' + Math.round(y * 1e4) + ',' + Math.round(z * 1e4);
const compOfKey = new Map(); for (let v = 0; v < nv; v++) compOfKey.set(keyOf(P[3 * v], P[3 * v + 1], P[3 * v + 2]), root(v));
return [...groups.entries()].map(([r, tris]) => {
const map = new Map(), pos = [], idx = [], tf = [];
const vid = (v) => { let k = map.get(v); if (k === undefined) { k = map.size; map.set(v, k); pos.push(P[3 * v], P[3 * v + 1], P[3 * v + 2]); } return k; };
for (const t of tris) { for (let k = 0; k < 3; k++) idx.push(vid(I[3 * t + k])); if (mesh.triFace) tf.push(mesh.triFace[t]); }
const out = { positions: Float32Array.from(pos), indices: Uint32Array.from(idx), triFace: mesh.triFace ? Int32Array.from(tf) : null };
if (mesh.edgeSegs) { const s = []; for (let k = 0; k < mesh.edgeSegs.length / 6; k++) if (compOfKey.get(keyOf(mesh.edgeSegs[6 * k], mesh.edgeSegs[6 * k + 1], mesh.edgeSegs[6 * k + 2])) === r) for (let c = 0; c < 6; c++) s.push(mesh.edgeSegs[6 * k + c]); out.edgeSegs = Float32Array.from(s); }
return out;
});
}
Density unit for the solver
64-bit
Converts the material's density from kg/m³ to the solver's tonnes per mm³.
ρ_t/mm³ = ρ_kg/m³ × 1e-12
One multiplication by 1e-12, which is not exactly representable; the product can differ from the correctly rounded value in its last digit.
nl-worker.js
line 29
PartCraft 0.8.1
const T_PER_KGM3 = 1e-12; // kg/m³ → tonnes/mm³ (the units of tetfea.js)
Von Mises stress from CalculiX results
64-bit
The von Mises stress at each node from the six stress components of a results file.
The components are the numbers parsed from the file by calculix.js, which this review did not cover; the formula is one 64-bit expression of squares under a square root and the result is kept in a Float64Array.
nl-worker.js
line 151
PartCraft 0.8.1
if (stressBlk) for (const [id, v] of stressBlk.values) { const n = id - 1; if (n >= 0 && n < P.model.n && v.length >= 6) { const [a, b, c, d, e, f] = v; nodeVm[n] = Math.sqrt(0.5 * ((a - b) ** 2 + (b - c) ** 2 + (c - a) ** 2) + 3 * (d * d + e * e + f * f)); } }
In detail
Method
Small-strain J2 plasticity with a bilinear curve and isotropic hardening; radial return and the consistent tangent (the file header cites Simo and Hughes, Computational Inelasticity, 1998, section 3.3 and box 3.1) on 10-node tetrahedra with a 4-point rule; load-stepped Newton–Raphson with each linear system solved by preconditioned conjugate gradients from tetfea.js.
Number types
Float64Array for displacements, forces, stresses and state; Float32Array for the per-element maxima and for the plot's vertex fields.
Tolerances
Newton: residual norm ÷ (λ × load norm) under 1e-6 by default, at most 25 iterations; with relaxed contact 1000 × that after the first iteration. Linear solve tolerance 1e-10; an unconverged linear solve is accepted when its residual is under 1e-6. A failed step is halved down to 1/64 of a step. Up to 12 contact re-pairings per step. Contact search cutoff max(2h, 1e-3).
Limits
Small displacement and strain, monotone loading, isotropic materials, as the header says. Pieces of a part are separated by welding vertices on a 1e-4 grid.
Approach: 32-bit step, flagged
One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.