Builds the tetrahedral analysis mesh of a solid from its surface triangles: element sizes, inside and outside, cut points on the surface, boundary vertices on the true faces, and the mesh quality and volume that are reported.
How the precision is kept
Node coordinates are Float64Array and every distance, volume and angle is 64-bit. Inside and outside, on-surface and ray hits are decided with floating-point tests and fixed tolerances scaled to the part's bounding-box diagonal, not with exact predicates. The mesh follows the tessellation handed in for faces that have no exact surface record, so the mesh and its size field set the accuracy of the analysis.
Calculations under this policy
20 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.
The nearest point of a surface triangle to a given point, and the squared distance.
Region tests on dot products decide vertex, edge or interior; the edge and interior cases each use one division by a sum of those products. The triangle's edge vectors are differences of absolute coordinates. The squared distance is returned, so no square root is taken here.
tetmesh.js
lines 79 to 93
PartCraft 0.8.1
// closest point on triangle t to (x,y,z) → squared distance, point in q (Ericson 5.1.5)
const triClosest = (t, x, y, z) => {
const a = 3 * I[3 * t], b = 3 * I[3 * t + 1], c = 3 * I[3 * t + 2];
const ax = P[a], ay = P[a + 1], az = P[a + 2], abx = P[b] - ax, aby = P[b + 1] - ay, abz = P[b + 2] - az, acx = P[c] - ax, acy = P[c + 1] - ay, acz = P[c + 2] - az;
const apx = x - ax, apy = y - ay, apz = z - az, d1 = abx * apx + aby * apy + abz * apz, d2 = acx * apx + acy * apy + acz * apz;
let u, v;
if (d1 <= 0 && d2 <= 0) { u = 0; v = 0; } else {
const bpx = x - P[b], bpy = y - P[b + 1], bpz = z - P[b + 2], d3 = abx * bpx + aby * bpy + abz * bpz, d4 = acx * bpx + acy * bpy + acz * bpz;
const cpx = x - P[c], cpy = y - P[c + 1], cpz = z - P[c + 2], d5 = abx * cpx + aby * cpy + abz * cpz, d6 = acx * cpx + acy * cpy + acz * cpz;
const vc = d1 * d4 - d3 * d2, vb = d5 * d2 - d1 * d6, va = d3 * d6 - d5 * d4;
if (d3 >= 0 && d4 <= d3) { u = 1; v = 0; } else if (d6 >= 0 && d5 <= d6) { u = 0; v = 1; } else if (vc <= 0 && d1 >= 0 && d3 <= 0) { u = d1 / (d1 - d3); v = 0; } else if (vb <= 0 && d2 >= 0 && d6 <= 0) { u = 0; v = d2 / (d2 - d6); } else if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) { const w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); u = 1 - w; v = w; } else { const den = 1 / (va + vb + vc); u = vb * den; v = vc * den; }
}
q[0] = ax + abx * u + acx * v; q[1] = ay + aby * u + acy * v; q[2] = az + abz * u + acz * v;
return (x - q[0]) ** 2 + (y - q[1]) ** 2 + (z - q[2]) ** 2;
};
Ray against a triangle
64-bit
Where a ray crosses a surface triangle, used for inside tests, wall thickness and cut points.
t = (e₂·q) ÷ det, det = e₁·(d × e₂)
A ray is treated as parallel when |det| ≤ 1e-14 × |e₁| |e₂| |d|, a tolerance scaled to the three lengths. The barycentric bounds have a slack of 1e-12, widened to 1e-7 for first-hit queries, and a hit within 1e-9 of an edge is flagged ambiguous rather than counted.
tetmesh.js
lines 133 to 145
PartCraft 0.8.1
const rayTri = (t, ox, oy, oz, dx, dy, dz) => {
const a = 3 * I[3 * t], b = 3 * I[3 * t + 1], c = 3 * I[3 * t + 2];
const e1x = P[b] - P[a], e1y = P[b + 1] - P[a + 1], e1z = P[b + 2] - P[a + 2], e2x = P[c] - P[a], e2y = P[c + 1] - P[a + 1], e2z = P[c + 2] - P[a + 2];
const px = dy * e2z - dz * e2y, py = dz * e2x - dx * e2z, pz = dx * e2y - dy * e2x, det = e1x * px + e1y * py + e1z * pz;
const scale = Math.hypot(e1x, e1y, e1z) * Math.hypot(e2x, e2y, e2z) * Math.hypot(dx, dy, dz);
if (Math.abs(det) <= 1e-14 * scale) return 0; // parallel: 0 = no hit
const inv = 1 / det, sx = ox - P[a], sy = oy - P[a + 1], sz = oz - P[a + 2], u = (sx * px + sy * py + sz * pz) * inv;
if (u < -tolB || u > 1 + tolB) return 0;
const qx = sy * e1z - sz * e1y, qy = sz * e1x - sx * e1z, qz = sx * e1y - sy * e1x, v = (dx * qx + dy * qy + dz * qz) * inv;
if (v < -tolB || u + v > 1 + tolB) return 0;
hit[0] = (e2x * qx + e2y * qy + e2z * qz) * inv; hit[1] = u; hit[2] = v;
return u < 1e-9 || v < 1e-9 || u + v > 1 - 1e-9 ? 2 : 1; // 2 = on an edge (ambiguous)
};
Is a point inside the part
64-bit
Decides whether a point lies inside the closed surface.
inside when most clean rays cross the surface an odd number of times
Up to three fixed skewed directions are tried; a ray that grazes a triangle edge is discarded, and the answer is the majority of the rest (outside when none is clean). This is a floating-point parity test, not an exact predicate.
tetmesh.js
lines 164 to 176
PartCraft 0.8.1
/** Inside the closed mesh? Parity of crossings along skewed rays; an edge-grazing ray is retried; majority of the clean ones. */
const inside = (p) => {
let votes = 0, clean = 0;
for (const D of DIRS) {
let cnt = 0, amb = false;
const L = 4 * diag + 1, d = [D[0] * L, D[1] * L, D[2] * L];
walk(p, d, 0, 1, (t, k) => { if (k === 2) { amb = true; return true; } cnt++; return false; });
if (amb) continue;
clean++; if (cnt & 1) votes++;
if (clean === 2 && (votes === 0 || votes === 2)) break;
}
return clean ? votes * 2 > clean : false;
};
Projection onto a face
64-bit
Moves a point onto a plane, a cylinder, a recorded surface or, failing those, the face's triangles.
plane: p − ((p − o)·n) n; cylinder: axis point + r × radial ÷ |radial|
The plane case assumes a unit normal and is one dot product. The cylinder case divides by the radial distance and is skipped when that distance is not above 1e-12 (absolute, in model units). Other surfaces go through closestOnSurface in surfaces.js.
tetmesh.js
lines 191 to 202
PartCraft 0.8.1
/** Project p onto B-rep face record f (plane / cylinder / surfaces.js record exactly; else the face's tessellation via bvh) → [x,y,z]. */
export function projectToFace(p, f, fi, bvh, triFace) {
if (f && f.type === 'plane' && f.n && f.o) { const d = (p[0] - f.o[0]) * f.n[0] + (p[1] - f.o[1]) * f.n[1] + (p[2] - f.o[2]) * f.n[2]; return [p[0] - d * f.n[0], p[1] - d * f.n[1], p[2] - d * f.n[2]]; }
if (f && f.type === 'cylinder' && f.axis && f.r > 0) {
const o = f.axis.o, a = f.axis.d, qx = p[0] - o[0], qy = p[1] - o[1], qz = p[2] - o[2], t = qx * a[0] + qy * a[1] + qz * a[2];
const rx = qx - t * a[0], ry = qy - t * a[1], rz = qz - t * a[2], rl = Math.hypot(rx, ry, rz);
if (rl > 1e-12) return [o[0] + t * a[0] + f.r * rx / rl, o[1] + t * a[1] + f.r * ry / rl, o[2] + t * a[2] + f.r * rz / rl];
}
if (f && f.surf) return closestOnSurface(f.surf, p).p;
const c = bvh.closest(p, triFace ? (t) => triFace[t] === fi : null);
return c.tri >= 0 ? c.p.slice() : p.slice();
}
Projection onto the meeting of several faces
64-bit
Moves a point onto the edge or corner where two or three faces meet.
repeat: project onto each face in turn, until the move < 10⁻¹² × diagonal
Alternating projections, at most 12 rounds, stopping when a full round moves the point less than 1e-12 × the bounding-box diagonal. The result is the last iterate, so for faces meeting at a shallow angle it may stop at the round limit before that tolerance.
tetmesh.js
lines 203 to 209
PartCraft 0.8.1
/** Onto the intersection of two or more faces: alternating projections (converges for faces meeting at an angle). */
export function projectToFaces(p, list, faces, bvh, triFace, iters = 12) {
let x = p.slice();
if (list.length === 1) return projectToFace(x, faces[list[0]], list[0], bvh, triFace);
for (let k = 0; k < iters; k++) { const x0 = x; for (const fi of list) x = projectToFace(x, faces[fi], fi, bvh, triFace); if (Math.hypot(x[0] - x0[0], x[1] - x0[1], x[2] - x0[2]) < 1e-12 * (bvh.diag || 1)) break; }
return x;
}
Signed volume of a tetrahedron
64-bit
The volume of one element, positive for the expected vertex order.
V = (u · (v × w)) ÷ 6
The three edge vectors are differences from the first vertex, then a 3 × 3 determinant by cofactors. Its sign is used directly to test for inverted elements, with no tolerance.
tetmesh.js
lines 215 to 220
PartCraft 0.8.1
/** Signed volume of tet (a,b,c,d) of coordinates X. */
export function tetVolume(X, a, b, c, d) {
const ax = X[3 * a], ay = X[3 * a + 1], az = X[3 * a + 2];
const ux = X[3 * b] - ax, uy = X[3 * b + 1] - ay, uz = X[3 * b + 2] - az, vx = X[3 * c] - ax, vy = X[3 * c + 1] - ay, vz = X[3 * c + 2] - az, wx = X[3 * d] - ax, wy = X[3 * d + 1] - ay, wz = X[3 * d + 2] - az;
return (ux * (vy * wz - vz * wy) - uy * (vx * wz - vz * wx) + uz * (vx * wy - vy * wx)) / 6;
}
Dihedral angles and aspect ratio of a tetrahedron
64-bit
The smallest and largest dihedral angle and the aspect ratio of one element.
Angles come from atan2 of a sine term and a cosine term, so there is no acos of a value near ±1 and a flat element reads 0° or 180°. An element with a zero-area face gets an angle of 0; a zero inradius gives an infinite aspect ratio.
tetmesh.js
lines 222 to 243
PartCraft 0.8.1
/** Min / max dihedral angle (degrees) and aspect ratio (longest edge ÷ (2√6 × inradius): 1 for a regular tet) of one tet.
* Angles from θ = atan2(3V·|e| / (2·A₁A₂), −n₁·n₂) with face normals from the vertex order, so a flat tet reads 0° / 180°. */
export function tetQuality(X, v) {
const p = v.map((i) => [X[3 * i], X[3 * i + 1], X[3 * i + 2]]);
const sub = (a, b) => [a[0] - b[0], a[1] - b[1], a[2] - b[2]], cross = (a, b) => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
const vol = Math.abs(tetVolume(X, ...v));
const sgn = tetVolume(X, ...v) >= 0 ? 1 : -1;
// outward normals (×2 area) of the face opposite each vertex, for a positively ordered tet
const F = [[1, 2, 3], [0, 3, 2], [0, 1, 3], [0, 2, 1]];
const fn = F.map(([a, b, c]) => cross(sub(p[b], p[a]), sub(p[c], p[a])).map((x) => x * sgn));
const area = fn.map((n) => Math.hypot(...n) / 2);
let mn = 180, mx = 0;
for (const [a, b, c, d] of DIH_EDGES) {
const n1 = fn[c], n2 = fn[d], A1 = area[c], A2 = area[d], L = Math.hypot(...sub(p[a], p[b]));
const cs = -(n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2]) / 4, sn = 1.5 * vol * L; // (A1A2 cos θ, A1A2 sin θ)
const ang = A1 > 0 && A2 > 0 ? Math.atan2(sn, cs) * 180 / Math.PI : 0;
if (ang < mn) mn = ang; if (ang > mx) mx = ang;
}
const rin = 3 * vol / (area[0] + area[1] + area[2] + area[3] || 1);
let lmax = 0; for (const [a, b] of DIH_EDGES) lmax = Math.max(lmax, Math.hypot(...sub(p[a], p[b])));
return { minDih: mn, maxDih: mx, aspect: rin > 0 ? lmax / (2 * Math.sqrt(6) * rin) : Infinity };
}
Mesh quality summary
64-bit
The smallest and largest dihedral angle, the worst and mean aspect ratio and a histogram, as shown with the study.
Minimum, maximum and a running sum over all elements in 64-bit; the histogram bins the minimum angle by Math.floor(angle ÷ 10), which only chooses a bin.
tetmesh.js
lines 245 to 258
PartCraft 0.8.1
/** Quality summary of a mesh: min/max dihedral, worst/mean aspect ratio, a histogram of min dihedral angles (10° bins). */
export function meshQuality(nodes, tets) {
const m = tets.length / 4, hist = new Array(10).fill(0);
let mn = 180, mx = 0, amax = 0, asum = 0, bad = 0;
const v = [0, 0, 0, 0];
for (let e = 0; e < m; e++) {
v[0] = tets[4 * e]; v[1] = tets[4 * e + 1]; v[2] = tets[4 * e + 2]; v[3] = tets[4 * e + 3];
const q = tetQuality(nodes, v);
if (q.minDih < mn) mn = q.minDih; if (q.maxDih > mx) mx = q.maxDih; if (q.aspect > amax) amax = q.aspect; asum += q.aspect;
hist[Math.min(9, Math.floor(q.minDih / 10))]++;
if (q.minDih < 5) bad++;
}
return { minDihedral: m ? mn : 0, maxDihedral: m ? mx : 0, maxAspect: amax, meanAspect: m ? asum / m : 0, under5: bad, hist };
}
Element size from surface curvature
64-bit
Estimates the radius of curvature between neighboring surface triangles of one face and the element size it allows.
angle = acos(n₁·n₂); r ≈ |c₁ − c₂| ÷ angle; size = r × curv
The acos argument is clamped to [−1, 1] and pairs turning by less than 1e-3 rad are skipped, so the division is guarded. The radius is an estimate from the tessellation; the code comment says it reads about a third of the true radius on doubly curved faces, which have their own exact path.
tetmesh.js
lines 295 to 298
PartCraft 0.8.1
const n1 = bvh.triN(o), n2 = bvh.triN(t), ang = Math.acos(Math.max(-1, Math.min(1, n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2])));
if (ang < 1e-3) continue;
const c1 = cen(o), c2 = cen(t), r = Math.hypot(c1[0] - c2[0], c1[1] - c2[1], c1[2] - c2[2]) / ang, s = r * curv;
if (s < h) samples.push([(P[3 * a] + P[3 * b]) / 2, (P[3 * a + 1] + P[3 * b + 1]) / 2, (P[3 * a + 2] + P[3 * b + 2]) / 2, Math.max(floor, s)]);
Element size from wall thickness and gaps
64-bit
Measures the wall behind and the gap in front of each surface triangle and the element size they allow.
A ray starts 1e-7 × the box diagonal off the triangle's centroid along its normal and runs 2 × h × thick; only triangles whose normal has a dot product below −0.5 with this triangle's normal count. The measured distance is the ray parameter times the ray length.
tetmesh.js
lines 305 to 312
PartCraft 0.8.1
const c = cen(t), n = bvh.triN(t), eps = 1e-7 * bvh.diag;
for (const sgn of [-1, 1]) {
const o = [c[0] + sgn * eps * n[0], c[1] + sgn * eps * n[1], c[2] + sgn * eps * n[2]], d = [sgn * n[0] * h * 2 * thick, sgn * n[1] * h * 2 * thick, sgn * n[2] * h * 2 * thick];
let best = Infinity;
bvh.walk(o, d, 0, 1, (tt, k, th) => { if (tt === t) return false; const m = bvh.triN(tt); if (m[0] * n[0] + m[1] * n[1] + m[2] * n[2] < -0.5 && th < best) best = th; return false; });
if (best === Infinity) continue;
const dist = best * h * 2 * thick, s = sgn < 0 ? dist / thick : dist / 1.5;
if (s < h) samples.push([c[0], c[1], c[2], Math.max(floor, s)]);
Background lattice spacing
64-bit
Chooses the cube size of the background lattice so the part's box is a whole number of cubes where possible.
s = L_min ÷ ceil(L_min ÷ h − 10⁻⁹); per axis: L ÷ n when within 8 % of s
The 1e-9 inside the ceiling keeps a side that is an exact multiple of h from gaining a cube through rounding. Lattice points are computed as box minimum + index × spacing, not by repeated addition, so they do not drift.
tetmesh.js
lines 363 to 388
PartCraft 0.8.1
/**
* BCC lattice over box (padded one cube) → { pts, tets: number[] (flat 4s), spacing, cubes }. The spacing is ≤ h: the
* box's thinnest side is a whole number of cubes, and each other side too when that stretches the cube ≤ 8 %; a side
* that would need more stays unaligned at its far end (its faces are cut, not lattice planes).
*/
export function bccLattice(bbox, h) {
const L = [0, 1, 2].map((a) => bbox.max[a] - bbox.min[a]), lmin = Math.max(Math.min(...L), 1e-9 * Math.max(...L));
const s0 = lmin / Math.ceil(lmin / h - 1e-9), s = s0 >= h / 2 ? s0 : h; // a side under h/2 is left to the size field
const n = [0, 1, 2].map((a) => { const k = Math.max(1, Math.round(L[a] / s)); return Math.abs(L[a] / (k * s) - 1) <= 0.08 ? k : Math.ceil(L[a] / s); });
const sp = [0, 1, 2].map((a) => (Math.abs(L[a] / (n[a] * s) - 1) <= 0.08 ? L[a] / n[a] : s) || s);
const pts = new Pts(), NC = n.map((v) => v + 3), CC = n.map((v) => v + 2); // corners i = −1 … n+1, centres −1 … n
const corner = (i, j, k) => (i + 1) + NC[0] * ((j + 1) + NC[1] * (k + 1));
for (let k = -1; k <= n[2] + 1; k++) for (let j = -1; j <= n[1] + 1; j++) for (let i = -1; i <= n[0] + 1; i++) pts.push(bbox.min[0] + i * sp[0], bbox.min[1] + j * sp[1], bbox.min[2] + k * sp[2]);
const c0 = pts.n, centre = (i, j, k) => c0 + (i + 1) + CC[0] * ((j + 1) + CC[1] * (k + 1));
for (let k = -1; k <= n[2]; k++) for (let j = -1; j <= n[1]; j++) for (let i = -1; i <= n[0]; i++) pts.push(bbox.min[0] + (i + 0.5) * sp[0], bbox.min[1] + (j + 0.5) * sp[1], bbox.min[2] + (k + 0.5) * sp[2]);
const tets = [];
const add = (a, b, c, d) => { if (tetVolume(pts.a, a, b, c, d) < 0) tets.push(a, b, d, c); else tets.push(a, b, c, d); };
for (let k = -1; k <= n[2]; k++) for (let j = -1; j <= n[1]; j++) for (let i = -1; i <= n[0]; i++) {
const c1 = centre(i, j, k);
// the faces towards +x, +y, +z neighbours: 4 tets (c1, c2, face edge)
if (i < n[0]) { const c2 = centre(i + 1, j, k), q = [corner(i + 1, j, k), corner(i + 1, j + 1, k), corner(i + 1, j + 1, k + 1), corner(i + 1, j, k + 1)]; for (let m = 0; m < 4; m++) add(c1, c2, q[m], q[(m + 1) % 4]); }
if (j < n[1]) { const c2 = centre(i, j + 1, k), q = [corner(i, j + 1, k), corner(i + 1, j + 1, k), corner(i + 1, j + 1, k + 1), corner(i, j + 1, k + 1)]; for (let m = 0; m < 4; m++) add(c1, c2, q[m], q[(m + 1) % 4]); }
if (k < n[2]) { const c2 = centre(i, j, k + 1), q = [corner(i, j, k + 1), corner(i + 1, j, k + 1), corner(i + 1, j + 1, k + 1), corner(i, j + 1, k + 1)]; for (let m = 0; m < 4; m++) add(c1, c2, q[m], q[(m + 1) % 4]); }
}
return { pts, tets, spacing: sp, cubes: n };
}
Longest edge chosen for bisection
64-bit
Picks the edge of an element that is split, the same way in every element that shares it.
compare round(length² ÷ unit × 10⁷); ties by vertex numbers
Squared lengths are divided by a reference squared length and rounded to an integer of 1e-7 units before comparing, so edges whose squared lengths differ by less than that are treated as equal and the tie is broken by vertex ids. The new vertex is the plain midpoint (a + b) ÷ 2.
tetmesh.js
lines 401 to 408
PartCraft 0.8.1
const longest = (T, e) => { // → index into PAIRS
let best = -1, bl = -1, bk = -1;
for (let p = 0; p < 6; p++) {
const a = T[4 * e + PAIRS[p][0]], b = T[4 * e + PAIRS[p][1]], l = Math.round(q2(a, b) / unit * 1e7), k = ekey(a, b);
if (l > bl || (l === bl && k < bk)) { bl = l; bk = k; best = p; }
}
return best;
};
Where a lattice edge crosses the surface
64-bit
The point on a lattice edge where it leaves the part.
p = a + t (b − a), t from the first ray hit; else 40 bisections of the inside test
The hit parameter comes from the ray and triangle test. When the ray reports no hit although the ends are classified in and out, 40 bisections on the classification give t to about 1e-12 of the edge. An edge that starts on the surface is probed 1e-4 of its length along.
tetmesh.js
lines 550 to 566
PartCraft 0.8.1
const cutPoint = (a, b) => { // a inside, b outside
const pa = pts.get(a), pb = pts.get(b), d = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]];
let tmin = 0;
if (onSurf[a]) {
const dl = 1e-4, q = [pa[0] + dl * d[0], pa[1] + dl * d[1], pa[2] + dl * d[2]];
if (!bvh.near(q, epsOn) && !bvh.inside(q)) return { t: 0, p: pa };
tmin = dl;
}
const hh = bvh.firstHit(pa, d, 1, tmin);
let t = hh ? hh.t : 0.5;
if (!hh) { // classification and the ray disagree (a grazing edge): bisect on the classification (on = in)
let lo = tmin, hi = 1;
for (let it = 0; it < 40; it++) { const m = (lo + hi) / 2, p = [pa[0] + m * d[0], pa[1] + m * d[1], pa[2] + m * d[2]]; if (bvh.near(p, epsOn) || bvh.inside(p)) lo = m; else hi = m; }
t = (lo + hi) / 2;
}
return { t, p: [pa[0] + t * d[0], pa[1] + t * d[1], pa[2] + t * d[2]] };
};
Lattice vertices warped onto the surface
64-bit
Moves a lattice vertex onto a nearby cut point so no very short edge is created.
warp when fraction of the edge < α (0.24999, 0.41189 or 0.33)
The fraction is the cut parameter itself and the edge class is chosen by comparing the edge's axis components scaled by the lattice spacing (thresholds 0.95 and 0.45). A vertex with several candidates takes the nearest cut point.
tetmesh.js
lines 575 to 578
PartCraft 0.8.1
const L = Math.hypot(...[0, 1, 2].map((x) => pts.a[3 * a + x] - pts.a[3 * b + x]));
for (const [v, f] of [[inn, c.t], [out, 1 - c.t]]) if (f < al) { const d = f * L, w = warpTo.get(v); if (!w || d < w[0]) warpTo.set(v, [d, c.p]); }
}
for (const [v, [, p]] of warpTo) { pts.a[3 * v] = p[0]; pts.a[3 * v + 1] = p[1]; pts.a[3 * v + 2] = p[2]; sign[v] = 0; }
Closest point on the part's edges
64-bit
Finds the nearest point on the part's edge polylines, used to snap boundary vertices onto sharp edges.
Every edge segment is tested in turn; the parameter is clamped to [0, 1] and a zero-length segment uses t = 0. One square root is taken of the smallest squared distance. The edge array is read as handed in.
tetmesh.js
lines 708 to 717
PartCraft 0.8.1
const segClosest = (p) => {
let best = Infinity, bp = null;
for (let s = 0; s < ns; s++) {
const ax = S[6 * s], ay = S[6 * s + 1], az = S[6 * s + 2], dx = S[6 * s + 3] - ax, dy = S[6 * s + 4] - ay, dz = S[6 * s + 5] - az, l2 = dx * dx + dy * dy + dz * dz;
let t = l2 > 0 ? ((p[0] - ax) * dx + (p[1] - ay) * dy + (p[2] - az) * dz) / l2 : 0; t = Math.max(0, Math.min(1, t));
const x = ax + t * dx, y = ay + t * dy, z = az + t * dz, d = (p[0] - x) ** 2 + (p[1] - y) ** 2 + (p[2] - z) ** 2;
if (d < best) { best = d; bp = [x, y, z]; }
}
return { d: Math.sqrt(best), p: bp };
};
Boundary vertices onto exact faces
64-bit
Moves each boundary vertex onto the true surface of its faces when all of them have an exact record.
move when 10⁻¹² × diagonal ≤ distance ≤ 0.25 h and the tets around keep min dihedral ≥ min(10°, before)
A move below 1e-12 × the box diagonal is skipped as already on the surface, one above a quarter of the element size is refused. The move is undone when any tet around the vertex inverts or falls below the angle limit.
tetmesh.js
lines 745 to 748
PartCraft 0.8.1
const p = [nodes[3 * v], nodes[3 * v + 1], nodes[3 * v + 2]], q = projectToFaces(p, list, faces, bvh, tf);
const d = Math.hypot(q[0] - p[0], q[1] - p[1], q[2] - p[2]);
if (d < 1e-12 * bvh.diag || d > 0.25 * h) continue;
if (moveTo(v, q, 10)) projected++;
Mesh volume
64-bit
The volume of the linear tetrahedral mesh, kept in its statistics.
V = Σₑ volume of tet e
A running 64-bit sum of the signed element volumes in element order.
tetmesh.js
line 754
PartCraft 0.8.1
let vol = 0; for (let e = 0; e < m; e++) vol += tetVolume(nodes, tets[4 * e], tets[4 * e + 1], tets[4 * e + 2], tets[4 * e + 3]);
Interior vertices moved to improve angles
64-bit
Moves interior vertices of poor elements when that raises the worst dihedral angle around them.
accept a trial position when the worst angle rises by more than 10⁻⁹ degrees
Trial positions are the neighbor average, fractions toward it, and steps off the plane of a sliver's opposite face. A trial is kept only if the worst angle improves by more than 1e-9°, and a position that inverts any tet scores −1. Three sweeps at most.
tetmesh.js
lines 766 to 806
PartCraft 0.8.1
/**
* Interior vertices of tets whose smallest dihedral angle is under `minDih`: moved to the centre of their neighbours,
* or part way there, when that raises the worst angle of the tets around them (3 sweeps). → vertices moved
*/
export function improveInterior(nodes, tets, btri, minDih = 10) {
const m = tets.length / 4, nn = nodes.length / 3, onB = new Uint8Array(nn);
for (const v of btri) onB[v] = 1;
const vt = Array.from({ length: nn }, () => []);
for (let e = 0; e < m; e++) for (let k = 0; k < 4; k++) vt[tets[4 * e + k]].push(e);
const q = (e) => tetQuality(nodes, [tets[4 * e], tets[4 * e + 1], tets[4 * e + 2], tets[4 * e + 3]]);
const worst = (v) => { let w = 180; for (const e of vt[v]) { if (tetVolume(nodes, tets[4 * e], tets[4 * e + 1], tets[4 * e + 2], tets[4 * e + 3]) <= 0) return -1; w = Math.min(w, q(e).minDih); } return w; };
let moved = 0;
for (let sweep = 0; sweep < 3; sweep++) {
const cand = new Set();
for (let e = 0; e < m; e++) if (q(e).minDih < minDih) for (let k = 0; k < 4; k++) { const v = tets[4 * e + k]; if (!onB[v]) cand.add(v); }
if (!cand.size) break;
for (const v of cand) {
const nb = new Set(); for (const e of vt[v]) for (let k = 0; k < 4; k++) { const w = tets[4 * e + k]; if (w !== v) nb.add(w); }
const c = [0, 0, 0]; for (const w of nb) for (let a = 0; a < 3; a++) c[a] += nodes[3 * w + a] / nb.size;
const p0 = [nodes[3 * v], nodes[3 * v + 1], nodes[3 * v + 2]], w0 = worst(v);
let best = w0, bp = null, ls = 0;
for (const w of nb) ls += Math.hypot(nodes[3 * w] - p0[0], nodes[3 * w + 1] - p0[1], nodes[3 * w + 2] - p0[2]) / nb.size;
const tries = [1, 0.5, 0.25].map((f) => [0, 1, 2].map((a) => p0[a] + f * (c[a] - p0[a])));
// and off the plane of each bad tet's opposite face (a sliver's vertex pushed out of its plane)
for (const e of vt[v]) {
if (q(e).minDih >= minDih) continue;
const o = [tets[4 * e], tets[4 * e + 1], tets[4 * e + 2], tets[4 * e + 3]].filter((x) => x !== v), P = o.map((x) => [nodes[3 * x], nodes[3 * x + 1], nodes[3 * x + 2]]);
const u = [0, 1, 2].map((a) => P[1][a] - P[0][a]), w2 = [0, 1, 2].map((a) => P[2][a] - P[0][a]), n = [u[1] * w2[2] - u[2] * w2[1], u[2] * w2[0] - u[0] * w2[2], u[0] * w2[1] - u[1] * w2[0]], nl = Math.hypot(...n) || 1;
const sgn = Math.sign(n[0] * (p0[0] - P[0][0]) + n[1] * (p0[1] - P[0][1]) + n[2] * (p0[2] - P[0][2])) || 1;
for (const f of [0.05, 0.1, 0.2, 0.3]) tries.push([0, 1, 2].map((a) => p0[a] + sgn * f * ls * n[a] / nl));
}
for (const t of tries) {
for (let a = 0; a < 3; a++) nodes[3 * v + a] = t[a];
const w = worst(v); if (w > best + 1e-9) { best = w; bp = t; }
}
const p = bp || p0; for (let a = 0; a < 3; a++) nodes[3 * v + a] = p[a];
if (bp) moved++;
}
}
return moved;
}
Jacobian check of a quadratic element
64-bit
The Jacobian determinant at the four corners and four Gauss points of a 10-node tetrahedron.
Eight 3 × 3 determinants by cofactors from sums of coordinate times shape derivative. The Gauss-point coordinates are the same 16-digit constants as in the solver.
tetmesh.js
lines 812 to 823
PartCraft 0.8.1
/** Jacobian determinants of a T10 at its corners and the 4 Gauss points (all must be > 0). */
export function t10DetJ(X, el) {
const pts = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [0.1381966011250105, 0.1381966011250105, 0.1381966011250105], [0.5854101966249685, 0.1381966011250105, 0.1381966011250105], [0.1381966011250105, 0.5854101966249685, 0.1381966011250105], [0.1381966011250105, 0.1381966011250105, 0.5854101966249685]];
return pts.map(([x, y, z]) => {
const L = [1 - x - y - z, x, y, z], dL = [[-1, -1, -1], [1, 0, 0], [0, 1, 0], [0, 0, 1]];
const J = [0, 0, 0, 0, 0, 0, 0, 0, 0];
const add = (node, d) => { for (let r = 0; r < 3; r++) for (let c = 0; c < 3; c++) J[3 * r + c] += X[3 * node + r] * d[c]; };
for (let i = 0; i < 4; i++) add(el[i], dL[i].map((v) => (4 * L[i] - 1) * v));
T10_EDGES.forEach(([i, j], k) => add(el[4 + k], [0, 1, 2].map((c) => 4 * (L[i] * dL[j][c] + L[j] * dL[i][c]))));
return J[0] * (J[4] * J[8] - J[5] * J[7]) - J[1] * (J[3] * J[8] - J[5] * J[6]) + J[2] * (J[3] * J[7] - J[4] * J[6]);
});
}
Midside nodes onto the exact surface
64-bit
Places the middle node of each boundary edge on the true surface so curved faces stay curved.
start at (a + b) ÷ 2; move onto the faces when every Jacobian stays > 0.2 × its straight value
A move below 1e-12 × the box diagonal is skipped. After the move the smallest of the eight Jacobian determinants of every element around the node must exceed 0.2 × its value before; otherwise the node goes back to the straight midpoint and is counted.
tetmesh.js
lines 858 to 865
PartCraft 0.8.1
const id = emap.get(key), p = [X[3 * id], X[3 * id + 1], X[3 * id + 2]], q = projectToFaces(p, [...fs], faces, bvh, mesh.triFace);
const d = Math.hypot(q[0] - p[0], q[1] - p[1], q[2] - p[2]);
if (d < 1e-12 * bvh.diag) continue;
// a moved midside node must not fold its elements (detJ > 0.2 × the straight element's, everywhere checked)
const els = around.get(id) || [], before = els.map((e) => Math.min(...t10DetJ(X, tets10.subarray(10 * e, 10 * e + 10))));
X[3 * id] = q[0]; X[3 * id + 1] = q[1]; X[3 * id + 2] = q[2];
if (els.every((e, i) => Math.min(...t10DetJ(X, tets10.subarray(10 * e, 10 * e + 10))) > 0.2 * before[i])) curved++;
else { X[3 * id] = p[0]; X[3 * id + 1] = p[1]; X[3 * id + 2] = p[2]; kept++; }
In detail
Method
Isosurface stuffing on a graded body-centered cubic lattice; the file header cites Labelle and Shewchuk 2007 and says the paper's angle bounds are not claimed. Longest-edge bisection for grading. Closest point on a triangle after Ericson 5.1.5; ray against triangle by Möller–Trumbore; inside test by ray parity along three fixed skewed directions. Boundary vertices are projected onto planes, cylinders and recorded surfaces; midside nodes of the quadratic elements likewise.
Number types
64-bit numbers; Float64Array for nodes, boxes and the size grid; Int32Array for elements; Int8Array and Uint8Array for vertex classes. The tessellation's vertex array is read as it is handed in.
Tolerances
On the surface: within 1e-8 × the box diagonal. Ray parallel to a triangle: |det| ≤ 1e-14 × |e₁| |e₂| |d|. Barycentric slack 1e-12 (1e-7 for first-hit queries); a hit within 1e-9 of a triangle edge is treated as ambiguous. Boxes padded by 1e-9 × the diagonal. Warp fractions 0.24999, 0.41189 and 0.33 of the edge. Quad diagonals tie within 1e-9 relative. Longest-edge ties: squared lengths compared after rounding to 1e-7 of a reference. Slivers: minimum dihedral under 5° dropped. Edge snapping within 0.45 × the local size, keeping 12°; projection skipped under 1e-12 × the diagonal or over 0.25 × the size, keeping 10°. A curved midside node is kept when every Jacobian stays above 0.2 × its straight value. Alternating projections stop at 1e-12 × the diagonal or 12 rounds. Lattice stretched at most 8 % per axis.
Limits
Faces without an exact record follow the tessellation and its chord error. The inside test is a majority of up to three rays; when every ray grazes an edge the point is taken as outside. A cut point that the ray misses is found by 40 bisections of the inside test.
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.