Moves the points and radii of a sketch until every relation and dimension is met, and reports how many degrees of freedom are left and which relations conflict or repeat.
How the precision is kept
The solver is iterative in 64-bit arithmetic and stops when every residual is within 1e-9 × max(1, sketch extent). That tolerance, not the number format, bounds how closely a solved sketch meets its dimensions. Lengths and the dimensionless residuals (the sine for parallel, the cosine for perpendicular, radians for an angle) are tested against the same number.
Open item. Give each kind of residual its own scale (lengths against the sketch size, sines and angles against a fixed small number) so that a large sketch does not loosen the angular relations.
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.
The derivative is a difference of two nearby function values, so it carries roughly ten digits rather than fifteen. It is used only to choose the step; the stopping test reads the residuals themselves.
solver.js
lines 52 to 67
PartCraft 0.8.1
/** Sparse Jacobian rows: [Map(varIndex → ∂f/∂x)] */
function jacobian(eqs, x) {
return eqs.map((q) => {
const row = new Map();
for (const k of q.deps) {
if (row.has(k)) continue;
const x0 = x[k], h = 1e-6 * Math.max(1, Math.abs(x0));
x[k] = x0 + h; const fp = q.f();
x[k] = x0 - h; const fm = q.f();
x[k] = x0;
const d = (fp - fm) / (2 * h);
if (d !== 0) row.set(k, d);
}
return row;
});
}
Cholesky solve
64-bit
Solves the symmetric system of each Newton step.
A = L·Lᵀ; L·y = b; Lᵀ·z = y
Plain 64-bit Cholesky with no pivoting. A diagonal value that is not positive returns null, and the caller then raises the damping and tries again.
solver.js
lines 76 to 97
PartCraft 0.8.1
/** Solve the symmetric positive (semi)definite m×m system by Cholesky; null if it breaks down. */
function cholSolve(A, b) {
const n = b.length;
const L = A.map((r) => r.slice());
for (let j = 0; j < n; j++) {
let s = L[j][j];
for (let k = 0; k < j; k++) s -= L[j][k] * L[j][k];
if (!(s > 0)) return null;
const d = Math.sqrt(s);
L[j][j] = d;
for (let i = j + 1; i < n; i++) {
let t = L[i][j];
for (let k = 0; k < j; k++) t -= L[i][k] * L[j][k];
L[i][j] = t / d;
}
}
const y = new Array(n);
for (let i = 0; i < n; i++) { let s = b[i]; for (let k = 0; k < i; k++) s -= L[i][k] * y[k]; y[i] = s / L[i][i]; }
const z = new Array(n);
for (let i = n - 1; i >= 0; i--) { let s = y[i]; for (let k = i + 1; k < n; k++) s -= L[k][i] * z[k]; z[i] = s / L[i][i]; }
return z;
}
Damped Newton iteration
64-bit
Repeats least-norm steps until every residual is within the tolerance.
dx = Jᵀ · (J·Jᵀ + μ·max(scale, 1e-12)·I)⁻¹ · (−r); x ← x + α·dx, α halved until Σr² falls
It returns ok only when every |residual| ≤ tol. A step is accepted only if the sum of squared residuals falls; when no step helps after the damping has been raised (up to 12 tries, or past 1e6) the loop ends with the best point found. The damping term added to the diagonal is at least 1e-300.
solver.js
lines 99 to 136
PartCraft 0.8.1
function newton(eqs, x, tol, maxIter = 100) {
let r = residuals(eqs);
let f = norm2(r);
let mu = 1e-12;
let it = 0;
for (; it < maxIter; it++) {
if (r.every((v) => Math.abs(v) <= tol)) return { ok: true, it };
const J = jacobian(eqs, x);
const m = J.length;
const A = [];
let scale = 0;
for (let i = 0; i < m; i++) {
A.push(new Array(m));
for (let j = 0; j <= i; j++) { const v = sparseDot(J[i], J[j]); A[i][j] = v; if (j < i) A[j][i] = v; }
scale = Math.max(scale, A[i][i]);
}
let improved = false;
for (let tries = 0; tries < 12 && !improved; tries++) {
const B = A.map((row, i) => row.map((v, j) => (i === j ? v + mu * Math.max(scale, 1e-12) + 1e-300 : v)));
const y = cholSolve(B, r.map((v) => -v));
if (!y) { mu *= 100; continue; }
const dx = new Array(x.length).fill(0);
J.forEach((row, i) => { for (const [k, v] of row) dx[k] += v * y[i]; });
const x0 = x.slice();
let alpha = 1;
for (let ls = 0; ls < 20; ls++) {
for (let k = 0; k < x.length; k++) x[k] = x0[k] + alpha * dx[k];
const r2 = residuals(eqs), f2 = norm2(r2);
if (f2 < f) { r = r2; f = f2; improved = true; break; }
alpha /= 2;
}
if (!improved) { for (let k = 0; k < x.length; k++) x[k] = x0[k]; mu *= 100; if (mu > 1e6) break; }
else mu = Math.max(1e-12, mu / 10);
}
if (!improved) return { ok: r.every((v) => Math.abs(v) <= tol), it };
}
return { ok: r.every((v) => Math.abs(v) <= tol), it };
}
Solving a sketch
64-bit
Solves the sketch in place and classifies it as under, fully or over defined.
tol = 1e-9 × max(1, extent)
The tolerance scales with the largest coordinate in the sketch and has a floor of 1e-9. If the solve fails, every point, radius and spline angle is restored from a copy taken at the start. A dragged radius is held at 1e-9 or more.
solver.js
lines 186 to 220
PartCraft 0.8.1
/**
* Solve the sketch in place. Pinned points/radii are held at the given values (dragging);
* when that cannot be met the drag is dropped and the sketch solved where it was. When the
* relations cannot all be met the geometry is left exactly as it was and the conflicts named.
*/
export function solveSketch(sk, { pin = null, pinRadius = null } = {}) {
const snap = sk.entities.map((e) => (e.type === 'point' ? [e.x, e.y] : e.type === 'circle' ? e.r : e.type === 'spline' ? [e.a0, e.a1] : null));
const restore = () => sk.entities.forEach((e, i) => {
if (e.type === 'point') { e.x = snap[i][0]; e.y = snap[i][1]; } else if (e.type === 'circle') e.r = snap[i];
else if (e.type === 'spline') for (const [j, k] of [[0, 'a0'], [1, 'a1']]) if (snap[i][j] !== undefined) e[k] = snap[i][j];
});
const tolScale = Math.max(1, extent(sk));
const tol = 1e-9 * tolScale;
let res = { ok: false, it: 0 };
let pinnedOk = false;
if ((pin && pin.size) || (pinRadius && pinRadius.size)) {
const pinned = new Set([...(pin?.keys() || []), ...(pinRadius?.keys() || [])]);
for (const e of sk.entities) {
if (pin?.has(e.id)) { [e.x, e.y] = pin.get(e.id); }
if (pinRadius?.has(e.id)) e.r = Math.max(1e-9, pinRadius.get(e.id));
}
const { vars, x } = buildVars(sk, pinned);
res = newton(equations(sk, x, vars), x, tol, 60);
if (res.ok) { writeBack(sk, vars, x); pinnedOk = true; } else restore();
}
if (!pinnedOk) {
const { vars, x } = buildVars(sk);
res = newton(equations(sk, x, vars), x, tol, 100);
if (res.ok) writeBack(sk, vars, x); else restore();
}
const a = analyze(sk, tolScale);
const geometry = sk.entities.some((e) => e.type !== 'point') || sk.entities.length > 0;
let status = !geometry ? 'empty' : a.conflicts.length || a.redundant.length ? 'over' : a.dof === 0 ? 'full' : 'under';
return { ok: res.ok, dragged: pinnedOk, status, dof: a.dof, conflicts: a.conflicts, redundant: a.redundant, determined: a.determined, iterations: res.it };
}
Writing the solution back
64-bit
Copies the solved variables to the sketch's points, radii and spline end angles.
angle ← atan2(sin a, cos a)
Coordinates and radii are copied unchanged. Spline end angles are passed through sin, cos and atan2 to bring them into (−π, π], which can change the last bits of an angle already in range.
solver.js
lines 37 to 47
PartCraft 0.8.1
function writeBack(sk, vars, x) {
for (const e of sk.entities) {
const v = vars.get(e.id);
if (v === undefined) continue;
if (e.type === 'point') { e.x = x[v[0]]; e.y = x[v[1]]; } else if (e.type === 'circle') e.r = x[v];
}
for (const e of sk.entities) {
if (e.type !== 'spline') continue;
for (const k of ['a0', 'a1']) { const v = vars.get(e.id + ':' + k); if (v !== undefined) e[k] = Math.atan2(Math.sin(x[v]), Math.cos(x[v])); }
}
}
Rank, redundancy and determined geometry
64-bit
Counts the degrees of freedom left, names relations that repeat or conflict, and marks which entities are fully defined.
dof = variables − rank; row redundant when its remainder after projection is under 1e-7; variable determined when Σ q_k² > 1 − 1e-6
Each Jacobian row is normalized and orthogonalized against the earlier rows twice, which limits the loss of orthogonality of a single pass. The decisions are threshold tests on 64-bit values, not exact rank computations, so a nearly dependent relation can be called either way.
solver.js
lines 138 to 184
PartCraft 0.8.1
/** Rank / redundancy / determined variables at the current geometry. */
export function analyze(sk, tolScale) {
const { vars, x } = buildVars(sk);
const eqs = equations(sk, x, vars);
const tol = 1e-8 * (tolScale ?? Math.max(1, extent(sk)));
const J = jacobian(eqs, x);
const n = x.length;
const Q = [];
const redundant = new Set();
const conflicts = new Set();
eqs.forEach((q, i) => {
const r = q.f();
if (Math.abs(r) > tol * 100) conflicts.add(q.cid);
const v = new Float64Array(n);
let vn = 0;
for (const [k, d] of J[i]) { v[k] = d; vn += d * d; }
vn = Math.sqrt(vn);
if (vn < 1e-12) { redundant.add(q.cid); return; }
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 !== 0) for (let k = 0; k < n; k++) v[k] -= d * b[k];
}
}
let wn = 0;
for (let k = 0; k < n; k++) wn += v[k] * v[k];
wn = Math.sqrt(wn);
if (wn < 1e-7) { redundant.add(q.cid); return; }
for (let k = 0; k < n; k++) v[k] /= wn;
Q.push(v);
});
const fixedVar = new Array(n).fill(0);
for (const b of Q) for (let k = 0; k < n; k++) fixedVar[k] += b[k] * b[k];
const det = (k) => fixedVar[k] > 1 - 1e-6;
const determined = new Set([ORIGIN]);
const ptDet = (id) => id === ORIGIN || (vars.has(id) && vars.get(id).every(det));
for (const e of sk.entities) {
let ok = pointsOf(e).every(ptDet);
if (e.type === 'circle') ok = ok && det(vars.get(e.id));
if (e.type === 'spline') for (const k of [':a0', ':a1']) if (vars.has(e.id + k)) ok = ok && det(vars.get(e.id + k));
if (ok) determined.add(e.id);
}
const strip = (s) => [...s].filter((c) => !c.startsWith('arc:'));
return { dof: n - Q.length, rank: Q.length, vars: n, redundant: strip(redundant), conflicts: strip(conflicts), determined };
}
The equation that keeps an arc's end point at the same distance from the center as its start point.
f = |e − c| − |s − c|
A difference of two Math.hypot values, in length units. At a solved sketch the two radii agree within the solver tolerance, not bit for bit.
sketch.js
lines 539 to 543
PartCraft 0.8.1
// arcs: end point on the same radius as the start point
for (const e of sk.entities) {
if (e.type !== 'arc') continue;
push('arc:' + e.id, depsE(e.id), () => { const c = P(e.c), s = P(e.s), q = P(e.e); return Math.hypot(q[0] - c[0], q[1] - c[1]) - Math.hypot(s[0] - c[0], s[1] - c[1]); });
}
Residuals: parallel and perpendicular
64-bit
The equations for two lines being parallel or square to each other.
parallel: (u × v) ÷ (|u|·|v|); perpendicular: (u · v) ÷ (|u|·|v|)
Both are divided by the two lengths, so they are the sine and the cosine of the angle and do not depend on the size of the lines. A zero length divides by 1e-300 in place of zero.
sketch.js
lines 565 to 566
PartCraft 0.8.1
case 'parallel': push(c.id, deps, () => { const u = dir(r[0]), v = dir(r[1]); return (u[0] * v[1] - u[1] * v[0]) / ((Math.hypot(...u) * Math.hypot(...v)) || 1e-300); }); break;
case 'perpendicular': push(c.id, deps, () => { const u = dir(r[0]), v = dir(r[1]); return (u[0] * v[0] + u[1] * v[1]) / ((Math.hypot(...u) * Math.hypot(...v)) || 1e-300); }); break;
Residuals: tangent
64-bit
The equations for a line touching a circle or arc, and for two circles or arcs touching.
shared end, line and arc: (u · (p − c)) ÷ (|u|·r); shared end, two arcs: ((p − a) × (p − b)) ÷ (rₐ·r_b); else: ±dist(c, line) − r, or d − (rₐ + r_b), or d − ±(rₐ − r_b)
Where the curves share an end point the code uses a form whose gradient is not zero at the solution (its comment explains why); those two forms are dimensionless. The general forms are differences of lengths, so two nearly equal lengths are subtracted at the solution.
sketch.js
lines 581 to 595
PartCraft 0.8.1
const j = sharedEnd(sk, E(r[0]), E(r[1]));
if (j && sig === 'LR') {
push(c.id, deps, () => { const u = dir(r[0]), ce = P(E(r[1]).c), p = P(j); return (u[0] * (p[0] - ce[0]) + u[1] * (p[1] - ce[1])) / ((Math.hypot(...u) * R(E(r[1]))) || 1e-300); });
break;
}
if (j) {
push(c.id, deps, () => { const p = P(j), a = P(E(r[0]).c), b = P(E(r[1]).c); return ((p[0] - a[0]) * (p[1] - b[1]) - (p[1] - a[1]) * (p[0] - b[0])) / ((R(E(r[0])) * R(E(r[1]))) || 1e-300); });
break;
}
if (sig === 'LR') push(c.id, deps, () => (c.sign || 1) * sd(P(E(r[1]).c), r[0]) - R(E(r[1])));
else push(c.id, deps, () => {
const a = E(r[0]), b = E(r[1]), ca = P(a.c), cb = P(b.c);
const d = Math.hypot(cb[0] - ca[0], cb[1] - ca[1]);
return c.kind === 'int' ? d - (c.sign || 1) * (R(a) - R(b)) : d - (R(a) + R(b));
});
Residuals: distance dimensions
64-bit
The equations that hold a distance between two points, the length of a line, or the distance from a point or a line to a line.
f = |b − a| − value; f = sign × signed distance − value
The dimension's value enters as the stored 64-bit number with no conversion. The residual is in sketch length units and is compared with the solver tolerance.
The equation that holds the angle between two lines at a value in degrees.
f = atan2(u × v, u · v) − sign × value × π ÷ 180, wrapped to (−π, π]
The value is converted to radians with Math.PI inside the residual; the stored dimension stays in degrees. The residual is in radians and is tested against the same tolerance as the lengths.
sketch.js
lines 651 to 654
PartCraft 0.8.1
let d = Math.atan2(u[0] * v[1] - u[1] * v[0], u[0] * v[0] + u[1] * v[1]) - (c.sign || 1) * c.value * Math.PI / 180;
while (d > Math.PI) d -= 2 * Math.PI;
while (d <= -Math.PI) d += 2 * Math.PI;
return d;
In detail
Method
Damped Newton with the least-norm step dx = Jᵀ(JJᵀ + μI)⁻¹(−r), a halving line search and growth of μ when a step does not help (the file header describes it); Jacobian by central differences; Cholesky solve; rank and redundancy by Gram–Schmidt done twice.
Number types
64-bit numbers in plain arrays; Float64Array for the orthonormal rows of the analysis.
Tolerances
Residual: 1e-9 × max(1, extent). Iterations: at most 100, or 60 while dragging. Difference step: 1e-6 × max(1, |x|). μ starts at 1e-12, × 100 on a failed step, stops past 1e6; up to 20 halvings in the line search. Analysis: conflict when |residual| > 1e-6 × max(1, extent); a row is redundant when its length is under 1e-12 or under 1e-7 after projection; a variable is determined when Σq² > 1 − 1e-6.
Limits
A solved dimension is met to within the tolerance, which for a sketch whose largest coordinate is 1000 is 1e-6. JJᵀ is formed explicitly, which squares the conditioning of the Jacobian. When the relations cannot all be met the geometry is put back exactly as it was.
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.