Finds joint values that make every closed chain of joints in an assembly meet, and reports how many values remain free.
How the precision is kept
An iterative solve in 64-bit: damped Newton with a least-norm step, a Jacobian by central differences, and a backtracking line search. It stops when the largest residual is below 1e-9 × max(1, L), where L is half the largest part bounding-box side, or after 60 steps. The answer's accuracy is set by that tolerance and by the finite-difference Jacobian.
Calculations under this policy
6 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.
Turns the rotation part of a mismatch into an axis times an angle in radians.
θ = acos(clamp((trace − 1) ÷ 2, −1, 1)); w = (R − Rᵀ)ᵛ × θ ÷ (2 sin θ)
The acos argument is clamped. Below θ = 1e-6 the vector is taken straight from the antisymmetric part, which avoids acos near 1 for the small rotations that matter at convergence. Within 1e-6 of a half turn the axis comes from the diagonal, with the square root guarded against a negative argument.
assembly.js
lines 295 to 309
PartCraft 0.8.1
/** The rotation of a 3 × 3 (in a column-major 4 × 4) as a rotation vector (axis × angle in radians). */
function rotVec(m) {
const c = Math.max(-1, Math.min(1, (m[0] + m[5] + m[10] - 1) / 2)), th = Math.acos(c);
const w = [m[6] - m[9], m[8] - m[2], m[1] - m[4]]; // 2 sin θ · axis
if (th < 1e-6) return w.map((x) => x / 2);
if (Math.PI - th < 1e-6) { // half a turn: the axis from the symmetric part
const d = [m[0], m[5], m[10]].map((x) => Math.sqrt(Math.max(0, (x + 1) / 2)));
const k = d.indexOf(Math.max(...d));
const ax = [0, 0, 0]; ax[k] = d[k];
const off = [[0, m[4], m[8]], [m[1], 0, m[9]], [m[2], m[6], 0]];
for (let q = 0; q < 3; q++) if (q !== k) ax[q] = off[k][q] / (2 * d[k]);
return ax.map((x) => x * th);
}
return w.map((x) => (x * th) / (2 * Math.sin(th)));
}
Loop residuals
64-bit
For each loop-closing joint, how far its B end is from where the loop puts it.
D = have⁻¹ × want; r = (D_x, D_y, D_z, w_x L, w_y L, w_z L)
The mismatch is formed as one relative transform, so the position terms are differences of nearly equal placements once the loop is close to closed. Rotation is scaled by L (half the largest bounding-box side, at least 5) so that radians and millimeters weigh alike.
assembly.js
lines 356 to 365
PartCraft 0.8.1
const resid = (xx, k = Infinity) => {
const { closing } = walk(asm, joints, frames, valsOf(xx));
const r = [];
for (const c of closing.slice(0, k)) {
const D = V.matMul(rigidInverse(c.have), c.want);
const w = rotVec(D);
r.push(D[12], D[13], D[14], w[0] * L, w[1] * L, w[2] * L);
}
return { r, closing };
};
Jacobian by central differences
64-bit
The rate of change of every residual with every unknown joint value.
J_qc = (r_q(x + h e_c) − r_q(x − h e_c)) ÷ 2h, h = 1e-6 × max(1, L)
Each entry subtracts two nearly equal residuals and divides by 2h, which multiplies their rounding error by 1 ÷ 2h. This affects the step, not the stopping test.
Moves the joint values toward closing the loops while changing them as little as possible.
dx = Jᵀ (J Jᵀ + μ D)⁻¹ (−r); accept x + α dx when Σr² falls, α = 1, ½, ¼ …
Stops when the largest absolute residual is at most 1e-9 × max(1, L), or after the step cap, or when 12 damping tries with 16 halvings each find no decrease. The normal matrix J Jᵀ squares the condition number of J; the damping adds μ × (diagonal + 1e-9 × largest diagonal) + 1e-15. Values are clamped to their limits before each trial.
assembly.js
lines 385 to 424
PartCraft 0.8.1
const newton = (x0, k = Infinity, cap = LOOP_MAX_ITER, cols = ALL) => {
const xs = x0.slice();
let cur = resid(xs, k), it = 0, mu = 1e-9;
const n = cols.length;
while (it < cap && n && worst(cur.r) > tol) {
it++;
const J = jac(xs, cur.r, k, cols), m = cur.r.length;
// the least-norm step dx = Jᵀ (J Jᵀ + μI)⁻¹ (−r) = (JᵀJ + μI)⁻¹ Jᵀ (−r): through whichever square is smaller (many loops
// give more equations than unknowns, most of them met twice over)
const wide = m > n, sz = wide ? n : m;
const A = Array.from({ length: sz }, () => new Float64Array(sz));
for (let a = 0; a < sz; a++) for (let b = 0; b <= a; b++) {
let s2 = 0;
if (wide) for (let q = 0; q < m; q++) s2 += J[q][a] * J[q][b];
else for (let c = 0; c < n; c++) s2 += J[a][c] * J[b][c];
A[a][b] = s2; A[b][a] = s2;
}
const g = wide ? Array.from({ length: n }, (_, c) => { let s2 = 0; for (let q = 0; q < m; q++) s2 -= J[q][c] * cur.r[q]; return s2; }) : cur.r.map((v) => -v);
let peak = 0;
for (let a = 0; a < sz; a++) peak = Math.max(peak, A[a][a]);
let improved = false;
for (let tries = 0; tries < 12 && !improved; tries++) {
const Am = Array.from({ length: sz }, (_, a) => Array.from(A[a], (v, b) => v + (a === b ? mu * (A[a][a] + 1e-9 * peak) + 1e-15 : 0)));
const y = V.solveDense(Am, g);
if (!y || !y.every(Number.isFinite)) { mu = Math.max(mu * 100, 1e-9); continue; }
let dx = y;
if (!wide) { dx = new Array(n).fill(0); for (let c = 0; c < n; c++) { let s2 = 0; for (let a = 0; a < m; a++) s2 += J[a][c] * y[a]; dx[c] = s2; } }
for (let alpha = 1, ls = 0; ls < 16; ls++, alpha /= 2) {
const xn = xs.slice();
cols.forEach((i, c) => { xn[i] += alpha * dx[c]; });
clamp(xn);
const rn = resid(xn, k);
if (norm2(rn.r) < norm2(cur.r)) { for (const i of cols) xs[i] = xn[i]; cur = rn; improved = true; break; }
}
if (improved) mu = Math.max(mu / 10, 1e-12); else mu = Math.max(mu * 100, 1e-9);
}
if (!improved) break;
}
return { x: xs, cur, it };
};
Angle on the nearest turn
64-bit
Returns each solved angle on the turn closest to where it started, inside its limits.
w = s₀ + ((val − s₀ + 180) mod 360) − 180
A modulo on degrees in 64-bit. Adding and subtracting 180 and s₀ around the modulo can change the last bits of the solved value; the size of that change is at the rounding level of the terms.
assembly.js
lines 446 to 456
PartCraft 0.8.1
// an angle comes back on the turn nearest where it started (the equations do not tell 50° from 410°), inside its limits
vars.forEach((v, c) => {
let a = best.x[c];
if (VALUE_UNIT[v.k] === 'deg') {
const s0 = x[c] / v.sc, val = a / v.sc;
let w = s0 + ((((val - s0 + 180) % 360) + 360) % 360) - 180;
if (w < v.lo || w > v.hi) { const up = w + 360, dn = w - 360; w = up >= v.lo && up <= v.hi ? up : dn >= v.lo && dn <= v.hi ? dn : val; }
a = w * v.sc;
}
x[c] = a;
});
Rank, mobility and conflicts
64-bit
Counts the independent conditions at the answer, how many values could still be set, and which loops did not close.
mobility = unknowns − rank; conflict when max |r| of a loop > max(1000 × tol, 1e-6)
Rank is decided by Gram–Schmidt with two passes on normalized Jacobian rows: a row is dropped when its norm is not above 1e-9 or its remainder is below 1e-6. Both thresholds are fixed numbers applied to rows of a finite-difference Jacobian, so a nearly dependent condition can be counted either way. A value is called free when its component in the row space is below 1 − 1e-6.
assembly.js
lines 458 to 486
PartCraft 0.8.1
// analysis at the answer
const r = cur.r, m = r.length;
const J = n ? jac(x, r) : [];
const Q = [];
let redundant = 0;
for (let a = 0; a < m; a++) {
const v = Float64Array.from(J[a] || []);
const vn = Math.hypot(...v);
if (!(vn > 1e-9)) { redundant++; continue; }
for (let k = 0; k < n; k++) v[k] /= vn;
for (let pass = 0; pass < 2; pass++) for (const b of Q) { let d = 0; for (let k = 0; k < n; k++) d += b[k] * v[k]; if (d) for (let k = 0; k < n; k++) v[k] -= d * b[k]; }
const wn = Math.hypot(...v);
if (wn < 1e-6) { redundant++; continue; }
for (let k = 0; k < n; k++) v[k] /= wn;
Q.push(v);
}
const det = new Float64Array(n);
for (const b of Q) for (let k = 0; k < n; k++) det[k] += b[k] * b[k];
const conflicts = [];
cur.closing.forEach((c, i) => { if (Math.max(...r.slice(i * 6, i * 6 + 6).map(Math.abs)) > Math.max(tol * 1e3, 1e-6)) conflicts.push(c.joint.id); });
const badGroups = new Set(conflicts.map((id) => group.get(id)));
const states = new Map();
for (const j of joints) {
if (!onLoop.has(j.id)) continue;
if (badGroups.has(group.get(j.id))) { states.set(j.id, 'over'); continue; }
const free = vars.some((v, i) => v.j === j && det[i] < 1 - 1e-6);
states.set(j.id, free ? 'under' : 'full');
}
const mobility = n - Q.length;
In detail
Method
Six equations per loop-closing joint (three of position in mm, three of the rotation vector times L). Step dx = Jᵀ(J Jᵀ + μI)⁻¹(−r) or (JᵀJ + μI)⁻¹Jᵀ(−r), through whichever square is smaller, with Levenberg–Marquardt damping and step halving. Rank and mobility by Gram–Schmidt on the Jacobian rows.
Number types
64-bit numbers; Float64Array for the Jacobian and the normal matrix.
Tolerances
Residual tolerance 1e-9 × max(1, L). Finite-difference step 1e-6 × max(1, L). Damping μ from 1e-9, at least 1e-12. A loop is a conflict when a residual exceeds max(1000 × tolerance, 1e-6). Jacobian row treated as empty below norm 1e-9 and as dependent below 1e-6 after projection. Caps: 12 steps on the first try, 60 per solve, 12 damping tries and 16 halvings per step.
Checked by
tools/partcraft-test/test-asm-loops.mjs calls placeAll and loopJoints on a four-bar linkage and compares the closing gap with 1e-6 mm. It loads the development build's copy of this file.
Limits
A central-difference Jacobian with step 1e-6 × max(1, L) carries truncation and rounding error, which slows convergence near the answer but does not change the stopping test. A start far from a closed position can end in a local minimum, reported as a conflict.
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.