Resources Numerical accuracy Calculation policy
Works out the gap between touching parts and the contact forces and stiffness that keep them from passing through each other.
A penalty method in 64-bit: parts are allowed to overlap slightly and are pushed apart in proportion to the overlap, so the answer depends on the penalty stiffness (20 × E ÷ h by default) as well as on the mesh. The closest point on a curved master face is searched on four flat sub-triangles. Node coordinates are Float64Array.
8 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.
Tests tools/partcraft-test/test-sim08.mjs
The contact stiffness used when none is given.
kn = 20 × E ÷ h
One expression. The factor 20 is a choice in the code; a larger factor leaves less penetration and a stiffer system to solve.
export const autoPenalty = (E, h, factor = 20) => factor * E / h;The point of a flat triangle nearest a given point, as weights of the three corners.
region tests on d₁ … d₆ = dot products of edges with corner-to-point vectors; inside: v = v_b ÷ (v_a + v_b + v_c)
The region is chosen by signs of dot products and their 2 × 2 combinations in ordinary 64-bit arithmetic, not by exact predicates; near a region border either branch returns nearly the same point. The interior case divides by v_a + v_b + v_c without a guard, which is zero only for a triangle of no area.
/** Closest point of triangle abc to p (Ericson, Real-Time Collision Detection § 5.1.5) → { d2, w: [wa, wb, wc] }. */
function closestTri(p, a, b, c) {
const ab = [b[0] - a[0], b[1] - a[1], b[2] - a[2]], ac = [c[0] - a[0], c[1] - a[1], c[2] - a[2]], ap = [p[0] - a[0], p[1] - a[1], p[2] - a[2]];
const dot = (x, y) => x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
const d1 = dot(ab, ap), d2 = dot(ac, ap);
let w;
if (d1 <= 0 && d2 <= 0) w = [1, 0, 0];
else {
const bp = [p[0] - b[0], p[1] - b[1], p[2] - b[2]], d3 = dot(ab, bp), d4 = dot(ac, bp);
if (d3 >= 0 && d4 <= d3) w = [0, 1, 0];
else {
const vc = d1 * d4 - d3 * d2;
if (vc <= 0 && d1 >= 0 && d3 <= 0) { const v = d1 / (d1 - d3); w = [1 - v, v, 0]; }
else {
const cp = [p[0] - c[0], p[1] - c[1], p[2] - c[2]], d5 = dot(ab, cp), d6 = dot(ac, cp);
if (d6 >= 0 && d5 <= d6) w = [0, 0, 1];
else {
const vb = d5 * d2 - d1 * d6;
if (vb <= 0 && d2 >= 0 && d6 <= 0) { const v = d2 / (d2 - d6); w = [1 - v, 0, v]; }
else {
const va = d3 * d6 - d5 * d4;
if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) { const v = (d4 - d3) / ((d4 - d3) + (d5 - d6)); w = [0, 1 - v, v]; }
else { const den = 1 / (va + vb + vc), v = vb * den, ww = vc * den; w = [1 - v - ww, v, ww]; }
}
}
}
}
}
const q = [0, 1, 2].map((i) => w[0] * a[i] + w[1] * b[i] + w[2] * c[i]);
return { d2: (p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2 + (p[2] - q[2]) ** 2, w };
}The nearest point of a six-node master face, in the face's own coordinates.
best of 4 flat sub-triangles (corners and midsides); (L₁, L₂) = weights × reference coordinates
The search is done on four flat triangles through the six nodes, and the weights found there are carried to the curved face. For a curved face the returned point is therefore near, but not exactly, the closest point of the quadratic surface.
/** Closest point on boundary T6 triangle t (current coordinates Y) to p → { d2, L1, L2 } (reference coordinates of the point). */
function closestT6(Y, v, p) {
let best = null;
const P = Array.from(v, (n) => [Y[3 * n], Y[3 * n + 1], Y[3 * n + 2]]);
for (const [i, j, k] of SUB) {
const r = closestTri(p, P[i], P[j], P[k]);
if (!best || r.d2 < best.d2) best = { d2: r.d2, L1: r.w[0] * REF[i][0] + r.w[1] * REF[j][0] + r.w[2] * REF[k][0], L2: r.w[0] * REF[i][1] + r.w[1] * REF[j][1] + r.w[2] * REF[k][1] };
}
return best;
}The position and outward unit normal of a six-node face at given face coordinates.
x = Σ N_k x_k; n = unit(Σ ∂N_k/∂L₁ x_k × Σ ∂N_k/∂L₂ x_k)
Sums of six terms and a cross product. A normal of zero length is divided by 1 instead, so the result is a zero vector and not NaN.
/** Position, T6 shape values and outward unit normal of boundary triangle v at (L1, L2) on coordinates Y. */
function surfacePoint(Y, v, L1, L2) {
const { N, d1, d2 } = t6(L1, L2), x = [0, 0, 0], a = [0, 0, 0], b = [0, 0, 0];
for (let k = 0; k < 6; k++) for (let c = 0; c < 3; c++) { const y = Y[3 * v[k] + c]; x[c] += N[k] * y; a[c] += d1[k] * y; b[c] += d2[k] * y; }
const 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]], l = Math.hypot(...n) || 1;
return { N, x, n: n.map((q) => q / l) };
}The area each contact point on a slave face stands for.
w = w_q × |∂x/∂L₁ × ∂x/∂L₂| ÷ 2
The rule's weight times the local area factor of the curved face, from the undeformed node positions. The weights are computed once when the contact is planned.
for (const [, L1, L2, wq] of TRI6) {
const { N, d1, d2 } = t6(L1, L2), a = [0, 0, 0], b = [0, 0, 0];
for (let k = 0; k < 6; k++) for (let c = 0; c < 3; c++) { a[c] += d1[k] * X[3 * sv[k] + c]; b[c] += d2[k] * X[3 * sv[k] + c]; }
items.push({ sv, N, w: wq * Math.hypot(a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]) / 2, cands });
}Finds the master point each slave point faces, and reports whether the pairing changed.
paired when d² ≤ cutoff²; moved when the face changes or |ΔL| > 2e-3
The threshold of 2e-3 is in the face's reference coordinates (0 to 1), so it scales with the face size. A bonded point keeps its first pairing. The outer loop in nonlinear.js repeats the solve until this reports no movement.
function pairAll(Y) {
let moved = false;
for (const pl of P) for (const it of pl.items) {
const p = [0, 0, 0]; for (let k = 0; k < 6; k++) for (let a = 0; a < 3; a++) p[a] += it.N[k] * Y[3 * it.sv[k] + a];
if (pl.type === 'bonded' && it.pair) continue; // a bonded point keeps the pairing it starts with
let best = null;
for (const t of it.cands) { const rr = closestT6(Y, bt.subarray(6 * t, 6 * t + 6), p); if (!best || rr.d2 < best.d2) best = { ...rr, t }; }
if (!best || !(best.d2 <= pl.cut * pl.cut)) { if (it.pair) moved = true; it.pair = null; continue; }
const { N, n } = surfacePoint(Y, bt.subarray(6 * best.t, 6 * best.t + 6), best.L1, best.L2), old = it.pair;
if (!old || old.t !== best.t || Math.abs(old.L1 - best.L1) > 2e-3 || Math.abs(old.L2 - best.L2) > 2e-3) moved = true;
it.pair = { t: best.t, L1: best.L1, L2: best.L2, N, n };
}
return moved;
}The normal gap at a contact point and the force and stiffness it adds.
g = n · (x_s − x_m); force = kn × w × (−g) when g < 0; stiffness = kn × w × c cᵀ ⊗ n nᵀ; pressure = force ÷ w
The gap is the difference of two nearly equal positions, each a sum of six shape-weighted node positions (rest position plus displacement), so its absolute error is at the rounding level of the coordinates, not of the gap. A point with gap above 1e-6 mm adds no force and 1e-4 of the stiffness.
// gap on the current shape with the frozen foot point and normal: n·(x_s − Σ N_i x_i)
const n = pr.n, xs = [0, 0, 0]; let g = 0;
for (let a = 0; a < 3; a++) { let d = 0; for (let q = 0; q < 12; q++) d += c[q] * (X0[3 * nodes[q] + a] + u[3 * nodes[q] + a]); g += d * n[a]; if (a >= 0) for (let q = 0; q < 6; q++) xs[a] += it.N[q] * (X0[3 * nodes[q] + a] + u[3 * nodes[q] + a]); }
const gtol = pl.gtol ?? 1e-6;
if (!(g <= gtol)) {
// an open point carries a token stiffness (1e-4 of the penalty) so a part that only the contact holds is not a
// singular system before it closes; no force
const kr = REGULAR * kn * w;
for (let q = 0; q < 12; q++) for (let q2 = 0; q2 < 12; q2++) addBlock(K, nodes[q], nodes[q2], (i, j) => kr * c[q] * c[q2] * n[i] * n[j]);
continue;
}
const kw = kn * w, f = g < 0 ? kw * -g : 0;
if (f > 0) { active++; fn += f; if (-g > maxPen) maxPen = -g; state.last.push({ gap: g, force: f, w, pressure: f / w, at: xs }); }
for (let q = 0; q < 12; q++) for (let a = 0; a < 3; a++) r[3 * nodes[q] + a] += f * c[q] * n[a];
for (let q = 0; q < 12; q++) for (let q2 = 0; q2 < 12; q2++) addBlock(K, nodes[q], nodes[q2], (i, j) => kw * c[q] * c[q2] * n[i] * n[j]);Holds a slave point to the master point it started nearest to.
rel = (x_s − x_m) − (x_s − x_m) at the first call; force = −kt × w × rel
The initial offset between the two points is stored on first use and subtracted afterwards, so an initial gap between the meshes does not create force. The tie has the penalty's compliance: the faces separate by stress ÷ kt.
if (bonded) {
const rel = [0, 0, 0];
for (let a = 0; a < 3; a++) for (let q = 0; q < 12; q++) rel[a] += c[q] * (X0[3 * nodes[q] + a] + u[3 * nodes[q] + a]);
const kw = kn * w;
for (let a = 0; a < 3; a++) rel[a] -= it.rest ? it.rest[a] : 0;
if (!it.rest) { it.rest = rel.slice(); for (let a = 0; a < 3; a++) rel[a] = 0; }
for (let q = 0; q < 12; q++) for (let a = 0; a < 3; a++) r[3 * nodes[q] + a] -= kw * c[q] * rel[a];
for (let q = 0; q < 12; q++) for (let q2 = 0; q2 < 12; q2++) addBlock(K, nodes[q], nodes[q2], (i, j) => (i === j ? kw * c[q] * c[q2] : 0));
active++;
continue;
}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.