Resources Numerical accuracy Calculation policy

Non-uniform scale of an exact body

Scales a solid by different factors along different axes while keeping it an exact boundary model.

How the precision is kept

Lines, planes and B-spline or Bézier geometry are mapped exactly by moving their poles in 64-bit. Circles, cylinders, cones, spheres, tori and surfaces of revolution cannot stay analytic under a non-uniform scale, so they are replaced by cubic B-splines that interpolate the true shape at equal angle steps. The step count is chosen so the deviation is at most 4e-7 mm, up to a cap of 720 steps; the deviation reached is returned and the shape's tolerances are raised to cover it.

Open item. Confirm how many digits the OpenCascade BRep text writer emits for this build, since every coordinate of a scaled body passes through that text.

Calculations under this policy

10 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.

brep-gtrsf.js PartCraft 0.8.1 64-bit 10 calculations Reviewed

Tests tools/partcraft-test/test-bodies.mjs

Number of steps round a circle

64-bit

How many interpolation points a round direction gets, from its effective radius.

h = (384 × DEV ÷ R)^¼; n = ceil(2π ÷ h) rounded up to a multiple of 4, between 32 and 720

Solves h⁴ R ÷ 384 = DEV for the step. The result is clamped to 720, so for a large enough radius the deviation exceeds DEV; the caller computes the deviation actually reached from the final n.

  • brep-gtrsf.js
  • lines 58 to 64
  • PartCraft 0.8.1
/** Periodic uniform cubic interpolation of the unit circle: n poles, pole i (0-based) at angle (i − 1)·h (OCCT's periodic alignment: the value at knot j·h comes from poles j … j+2). */
function circleN(rEff) {
  const h = Math.pow(384 * DEV / Math.max(rEff, 1e-9), 0.25);
  let n = Math.ceil(2 * Math.PI / h);
  n = Math.max(N_MIN, Math.min(N_MAX, Math.ceil(n / 4) * 4));
  return n;
}

Cubic spline through the unit circle

64-bit

The poles of a periodic cubic B-spline that passes through the circle at every knot.

pole i = k × (cos((i − 1)h), sin((i − 1)h)), k = 6 ÷ (4 + 2 cos h), err = h⁴ ÷ 384

With this k the spline value at a knot, (P₋₁ + 4P₀ + P₁) ÷ 6, lies on the unit circle, so the curve interpolates the circle at the knots. The error figure h⁴ ÷ 384 is the leading term of the spline's radial deviation midway between knots; it is a formula, not a measurement of the built curve.

  • brep-gtrsf.js
  • lines 65 to 70
  • PartCraft 0.8.1
function unitCircle(n) {
  const h = 2 * Math.PI / n, k = 6 / (4 + 2 * Math.cos(h));
  const cs = [];
  for (let i = 0; i < n; i++) { const a = (i - 1) * h; cs.push([k * Math.cos(a), k * Math.sin(a)]); }
  return { n, h, cs, err: Math.pow(h, 4) / 384 };
}

Largest stretch of the map

64-bit

The largest factor by which the scale lengthens any direction, used to size the deviation.

σ = |M x| after 60 power iterations of x ← unit(Mᵀ M x)

A fixed 60 iterations from the start vector (1, 0.7, 0.3) with no convergence test. Power iteration approaches the largest singular value from below, so the estimate can be slightly low when two stretch factors are close.

  • brep-gtrsf.js
  • lines 51 to 56
  • PartCraft 0.8.1
/** Largest singular value of a 3×3 (row-major in a 12-array) by power iteration on mᵀm. */
function sigmaMax(m) {
  let x = [1, 0.7, 0.3];
  for (let k = 0; k < 60; k++) { const y = av(m, x); const z = [m[0] * y[0] + m[4] * y[1] + m[8] * y[2], m[1] * y[0] + m[5] * y[1] + m[9] * y[2], m[2] * y[0] + m[6] * y[1] + m[10] * y[2]]; x = norm(z); }
  return len(av(m, x));
}

Image of a curve

64-bit

A line, circle, ellipse, Bézier or B-spline curve after the scale, as a B-spline in the same parameter.

line, Bézier, B-spline: poles ↦ T × pole (exact); circle, ellipse: T × (c + r₁ x cos + r₂ y sin) at the spline's poles

Poles are mapped by one 3 × 4 matrix product each. Weights and knots are passed through unchanged. For circles and ellipses the deviation bound is the unit-circle error times the larger radius times the largest stretch.

  • brep-gtrsf.js
  • lines 162 to 183
  • PartCraft 0.8.1
/**
 * A curve (local coordinates) under the affine map T (12-array, local → image) as a B-spline in the same parameter.
 * range(dir, origin) → [a, b]: the parameter range a straight direction must cover. → { rec: {deg, periodic, poles, w,
 * knots}, err } (err = the interpolation bound in mm, 0 when exact).
 */
function mapCurve(c, T, range) {
  const s = sigmaMax(T);
  if (c.type === 'line') {
    const d = c.d, [a, b] = range(d, c.p);
    return { rec: { deg: 1, periodic: false, poles: [ap(T, add(c.p, mul(d, a))), ap(T, add(c.p, mul(d, b)))], w: null, knots: [[a, 2], [b, 2]] }, err: 0 };
  }
  if (c.type === 'circle' || c.type === 'ellipse') {
    const r1 = c.type === 'circle' ? c.r : c.r1, r2 = c.type === 'circle' ? c.r : c.r2;
    const U = unitCircle(circleN(Math.max(r1, r2) * s));
    const poles = U.cs.map(([x, y]) => ap(T, add(c.p, add(mul(c.x, r1 * x), mul(c.y, r2 * y)))));
    return { rec: { deg: 3, periodic: true, poles, w: null, knots: periodicDir(U.n).knots }, err: U.err * Math.max(r1, r2) * s };
  }
  if (c.type === 'bezier') return { rec: { deg: c.deg, periodic: false, poles: c.poles.map((p) => ap(T, p)), w: c.w, knots: [[0, c.deg + 1], [1, c.deg + 1]] }, err: 0 };
  if (c.type === 'bspline') return { rec: { deg: c.deg, periodic: c.periodic, poles: c.poles.map((p) => ap(T, p)), w: c.w, knots: c.knots }, err: 0 };
  if (c.type === 'trimmed') return mapCurve(c.basis, T, range);
  throw new Error(`a ${c.type} curve`);   // i18n-ignore: internal, caught and reworded below
}

Image of a surface

64-bit

A plane, cylinder, cone, sphere, torus, extrusion, revolution or spline surface after the scale.

plane ↦ plane with frame from T x, T x × T y; cylinder, cone: ring(u) × (R + v sin α) + v cos α Z; sphere, torus: ring(u) × (R₁ + R₂ cos v) + R₂ sin v Z

Planes and spline surfaces are mapped exactly. Round surfaces become tensor-product B-splines: cubic interpolation round the axis, degree 1 along straight directions. For spheres and tori the stated bound adds the errors of the two round directions and multiplies by 1.5, a factor fixed in the code.

  • brep-gtrsf.js
  • lines 198 to 253
  • PartCraft 0.8.1
/** A surface under T as a B-spline in the same parameters → { rec | plane, err }. */
function mapSurface(S, T, range) {
  const s = sigmaMax(T);
  if (S.type === 'plane') {
    // the image of a plane is a plane: an orthonormal frame with the same handedness (pcurves on it are recomputed)
    const hand = Math.sign(dot(cross(S.x, S.y), S.n)) || 1;
    const X = norm(av(T, S.x)), Z = norm(cross(av(T, S.x), av(T, S.y)));
    const Y = cross(Z, X);
    return { plane: { p: ap(T, S.p), n: mul(Z, hand), x: X, y: Y }, err: 0 };
  }
  if (S.type === 'rtrimmed') return mapSurface(S.basis, T, range);
  if (S.type === 'bezsurf') {
    const dir = (d) => ({ deg: d, periodic: false, knots: [[0, d + 1], [1, d + 1]] });
    return { rec: { u: dir(S.ud), v: dir(S.vd), poles: S.poles.map((r) => r.map((p) => ap(T, p))), w: S.w }, err: 0 };
  }
  if (S.type === 'bsurf') return { rec: { u: { deg: S.ud, periodic: S.up, knots: S.uk }, v: { deg: S.vd, periodic: S.vp, knots: S.vk }, poles: S.poles.map((r) => r.map((p) => ap(T, p))), w: S.w }, err: 0 };
  const tensor = (udir, vdir, f, w = null) => ({ u: udir, v: vdir, poles: Array.from({ length: udir.n }, (_, i) => Array.from({ length: vdir.n }, (__, j) => ap(T, f(i, j)))), w });
  const { p: O, x: X, y: Y, n: Z } = S;
  const ring = (c) => add(mul(X, c[0]), mul(Y, c[1]));
  if (S.type === 'cylinder' || S.type === 'cone') {
    // O + (R + v·sin a)(cos u X + sin u Y) + v·cos a Z: periodic in u, straight (degree 1, exact) in v
    const sa = S.type === 'cone' ? Math.sin(S.ang) : 0, ca = S.type === 'cone' ? Math.cos(S.ang) : 1;
    const [a, b] = range(mul(Z, ca), O);
    const rmax = Math.max(Math.abs(S.r + a * sa), Math.abs(S.r + b * sa));
    const U = unitCircle(circleN(rmax * s));
    const vs = [a, b];
    const rec = tensor({ ...periodicDir(U.n), n: U.n }, { ...linearDir(a, b), n: 2 }, (i, j) => add(O, add(mul(ring(U.cs[i]), S.r + vs[j] * sa), mul(Z, vs[j] * ca))));
    return { rec, err: U.err * rmax * s };
  }
  if (S.type === 'sphere' || S.type === 'torus') {
    // both directions periodic: (cos v, sin v) interpolated over a full turn (the face only uses v in [−π/2, π/2] on a sphere)
    const R1 = S.type === 'sphere' ? 0 : S.r1, R2 = S.type === 'sphere' ? S.r : S.r2;
    const U = unitCircle(circleN((R1 + R2) * s)), W = unitCircle(circleN(R2 * s));
    const rec = tensor({ ...periodicDir(U.n), n: U.n }, { ...periodicDir(W.n), n: W.n }, (i, j) => add(O, add(mul(ring(U.cs[i]), R1 + R2 * W.cs[j][0]), mul(Z, R2 * W.cs[j][1]))));
    return { rec, err: (U.err * (R1 + R2) + W.err * R2) * s * 1.5 };
  }
  if (S.type === 'extrusion') {
    // C(u) + v·D
    const cd = resolveDir(curveDir(S.basis, range), s);
    const [a, b] = range(S.d, cd.poles[0], { spread: Math.max(...cd.poles.map((q) => len(sub(q, cd.poles[0])))) });
    const vs = [a, b];
    const rec = tensor({ ...cd.dir, n: cd.poles.length }, { ...linearDir(a, b), n: 2 }, (i, j) => add(cd.poles[i], mul(S.d, vs[j])), cd.w ? cd.poles.map((_, i) => [cd.w[i], cd.w[i]]) : null);
    return { rec, err: (cd.err || 0) * s };
  }
  if (S.type === 'revolution') {
    // P + Rot(D, u)(C(v) − P): each pole of the basis curve swept round a periodic cubic circle
    const D = norm(S.d), P = S.p;
    const cd = resolveDir(curveDir(S.basis, (d, o) => range(d, o, { revolve: P })), s);
    const rmax = Math.max(1e-9, ...cd.poles.map((q) => len(sub(sub(q, P), mul(D, dot(sub(q, P), D))))));
    const U = unitCircle(circleN(rmax * s));
    const rot = (q, c, sn) => { const w = sub(q, P), along = mul(D, dot(w, D)), rad = sub(w, along), perp = cross(D, rad); return add(P, add(along, add(mul(rad, c), mul(perp, sn)))); };
    const rec = tensor({ ...periodicDir(U.n), n: U.n }, { ...cd.dir, n: cd.poles.length }, (i, j) => rot(cd.poles[j], U.cs[i][0], U.cs[i][1]), cd.w ? Array.from({ length: U.n }, () => [...cd.w]) : null);
    return { rec, err: (U.err * rmax + (cd.err || 0)) * s };
  }
  throw new Error(`a ${S.type} surface`);   // i18n-ignore: internal, caught and reworded below
}

Inverse of a placement

64-bit

Inverts a 3 × 4 placement matrix found in the shape.

M⁻¹ = adj(M) ÷ det(M); t⁻¹ = −M⁻¹ t

Cofactor formula with a division by the determinant and no guard for a zero determinant. The determinant is a sum of products that can cancel for a nearly singular matrix; placements inside a shape are normally rigid motions.

  • brep-gtrsf.js
  • lines 32 to 38
  • PartCraft 0.8.1
const invM = (m) => {
  const [a, b, c, , d, e, f, , g, h, i] = m;
  const det = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
  const r = [(e * i - f * h) / det, (c * h - b * i) / det, (b * f - c * e) / det, 0, (f * g - d * i) / det, (a * i - c * g) / det, (c * d - a * f) / det, 0, (d * h - e * g) / det, (b * g - a * h) / det, (a * e - b * d) / det, 0];
  for (let k = 0; k < 3; k++) r[k * 4 + 3] = -(r[k * 4] * m[3] + r[k * 4 + 1] * m[7] + r[k * 4 + 2] * m[11]);
  return r;
};

Parameter range of straight directions

64-bit

How far a line or the straight direction of a surface must extend to cover the body.

t₀ = ((center − o) · d) ÷ |d|²; range = t₀ ± (1.02 × D ÷ |d| + 1)

D is the bounding-box diagonal, enlarged for revolutions and extrusions. The range is deliberately generous (2 % plus 1); the mapped line or ruling is exact over it, so the margin does not affect accuracy.

  • brep-gtrsf.js
  • lines 320 to 330
  • PartCraft 0.8.1
const rangeFor = (M) => {
  const centre = ap(invM(M), mid);
  return (d, o, extra) => {
    const L = len(d), t0 = dot(sub(centre, o), d) / (L * L);
    let D = diag;
    if (extra?.revolve) D = diag + 2 * len(sub(centre, extra.revolve));
    if (extra?.spread !== undefined) D = 2 * diag + extra.spread;
    D /= L;
    return [t0 - 1.02 * D - 1, t0 + 1.02 * D + 1];
  };
};

Raised shape tolerances

64-bit

The tolerance written on every vertex, edge and face of the scaled body.

tol = max(old, max(1e-7, 3 × largest deviation + 1e-9))

The tolerance is never lowered. With no approximated geometry the floor is 1e-7 mm.

  • brep-gtrsf.js
  • lines 401 to 406
  • PartCraft 0.8.1
// tolerances: raised to cover the interpolation (vertices, edges, faces)
const tl = tolUp();
for (const o of out) {
  if (o.type === 'Ve') o.geo[0] = fmt(Math.max(Number(o.geo[0]), Number(tl)));
  if (o.type === 'Ed') { const h = o.geo[0].trim().split(/\s+/); h[0] = fmt(Math.max(Number(h[0]), Number(tl))); o.geo[0] = ' ' + h.join(' '); }
  if (o.type === 'Fa') { const t = o.geo[0].split(/\s+/).filter(Boolean); t[1] = fmt(Math.max(Number(t[1]), Number(tl))); o.geo[0] = `${t[0]}  ${t[1]} ${t[2]} ${t[3]}`; }

Numbers written to the shape text

64-bit

Turns each coordinate, weight and knot into the text OpenCascade reads back.

text = String(v); −0 → 0

String(v) gives the shortest decimal that parses back to the same 64-bit number, so nothing is rounded on the way out. Negative zero is written as 0.

  • brep-gtrsf.js
  • line 50
  • PartCraft 0.8.1
const fmt = (v) => (Object.is(v, -0) ? '0' : String(v));

Matching equal placements

64-bit

Decides when two uses of the same curve or surface share one mapped copy.

key = entries rounded to 1e-12, joined

Matrix entries are rounded to 12 decimals for the lookup key only. Two placements that differ by less than that share the geometry mapped with the first one's matrix.

  • brep-gtrsf.js
  • line 49
  • PartCraft 0.8.1
const key = (m) => m.map((v) => Math.round(v * 1e12) / 1e12).join(',');

In detail

Method
Rewrites the shape's OpenCascade BRep text: every curve and surface is replaced by its image in the same parameters. Round directions use a periodic uniform cubic B-spline through the circle at 2π ÷ n steps, pole radius factor 6 ÷ (4 + 2 cos h); the header gives the deviation as h⁴ R ÷ 384.
Number types
64-bit numbers parsed from and written to decimal text. Written values use JavaScript's shortest text that reads back to the same 64-bit number.
Tolerances
Target deviation DEV = 4e-7 mm. Steps n between 32 and 720, a multiple of 4. Shape tolerances raised to max(1e-7, 3 × largest deviation + 1e-9). Placements that agree to 1e-12 per matrix entry share one mapped copy.
Limits
At the cap of 720 steps the bound h⁴ R ÷ 384 equals 4e-7 mm at an effective radius of about 26 m; larger round faces deviate more, and the larger value is what is returned. The digits of the incoming geometry are those the OpenCascade text writer emits, which this file does not control. Offset curves and surfaces, parabolas and hyperbolas are refused.

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.

Every policy that takes this approach