diff --git a/Cargo.lock b/Cargo.lock index 703fb31e..1e74cafe 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3368,6 +3368,7 @@ dependencies = [ "serde", "serde_json", "solidcraft-geom", + "spade", "thiserror 2.0.21", "truck-meshalgo", "truck-modeling", diff --git a/Cargo.toml b/Cargo.toml index 9e90ab7b..bc4f8618 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -37,6 +37,9 @@ truck-shapeops = { version = "0.4", default-features = false } truck-meshalgo = { version = "0.4", default-features = false, features = ["tessellation", "analyzers"] } truck-stepio = { version = "0.3", default-features = false } truck-topology = "0.6" +# constrained Delaunay + Ruppert refinement for trimmed faces (MIT OR Apache-2.0); already in the +# tree through truck-meshalgo. +spade = { version = "2.15", default-features = false, features = ["std"] } # packages (3MF) zip = { version = "3", default-features = false, features = ["deflate"] } diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index f15ce9ec..4486ddbb 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -16,6 +16,7 @@ pub mod fill; mod mesh; mod plane; mod profile; +pub mod tess; pub mod tspline; mod vec; diff --git a/crates/geom/src/tess.rs b/crates/geom/src/tess.rs new file mode 100644 index 00000000..3bd91aae --- /dev/null +++ b/crates/geom/src/tess.rs @@ -0,0 +1,495 @@ +//! Display tessellation plans: the parameter-space triangulation of every face of a body, +//! with the surfaces to evaluate it on, kept apart from the evaluation itself. +//! +//! The kernel decides *where* the vertices go (in each face's (u, v) parameters, with the +//! shared edge nodes fixed in 3D so neighbouring faces meet exactly) and *how* they are +//! connected. Evaluating the surfaces at those parameters is the bulk of the work and has no +//! branching on topology, so it can run anywhere: [`TessPlan::evaluate`] does it on the CPU, +//! and the GPU evaluator runs the same maths in a compute shader. Both evaluate the +//! [`EvalSurface`] forms defined here, so the two cannot drift apart. +//! +//! References: Piegl & Tiller, *The NURBS Book* (2nd ed.), A2.1 (span), A2.3 (basis +//! functions and derivatives), §4.5 (derivatives of a rational surface). + +use crate::bspline::MAX_SPLINE_DEGREE; +use crate::{Mesh, Vec3}; + +/// Most knots in one evaluable knot vector (bounds the shader's span search). +pub const MAX_EVAL_KNOTS: usize = 1 << 12; +/// Most control points in one evaluable surface. +pub const MAX_EVAL_POINTS: usize = 1 << 18; +/// Basis slots for the highest evaluable degree. +const W: usize = MAX_SPLINE_DEGREE + 1; + +/// A rational B-spline curve: homogeneous control points `[x·w, y·w, z·w, w]`. Outside its +/// domain it continues as the polynomial of the end span (a line extends for ever). +#[derive(Clone, Debug, PartialEq)] +pub struct EvalCurve { + pub degree: usize, + pub knots: Vec, + pub ctrl: Vec<[f64; 4]>, +} + +/// A surface in a form both evaluators understand. +#[derive(Clone, Debug, PartialEq)] +pub enum EvalSurface { + /// `origin + u·du + v·dv`. + Plane { origin: Vec3, du: Vec3, dv: Vec3 }, + /// A rational tensor B-spline; control point `(i, j)` is `ctrl[j·count[0] + i]`. + Nurbs { degree: [usize; 2], knots: [Vec; 2], count: [usize; 2], ctrl: Vec<[f64; 4]> }, + /// `profile(u)` turned by the angle `v` (right-handed) about the axis through `origin` + /// along the unit `axis`, then mapped by the affine `transform` (row-major 3×4); with + /// `swap` the parameters are exchanged first (the surface is evaluated at `(v, u)`). + Revolved { profile: EvalCurve, origin: Vec3, axis: Vec3, transform: [[f64; 4]; 3], swap: bool }, +} + +/// One face's share of a [`TessPlan`]. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct FacePlan { + /// The B-rep face (index in the body's face order). + pub face: u32, + /// Index into [`TessPlan::surfaces`]; `None` when every position and normal is given. + pub surface: Option, + /// Parameters of every vertex. + pub uv: Vec<[f64; 2]>, + /// Positions of the first `fixed.len()` vertices: the shared edge nodes, used exactly as + /// they are so the neighbouring face gets bit-identical boundary points. + pub fixed: Vec, + /// Normals of every vertex, only for faces without an evaluable surface. + pub normals: Vec, + /// The face is reversed against its surface: normals point the other way. + pub flip: bool, + /// Counter-clockwise seen from outside (orientation already applied). + pub triangles: Vec<[u32; 3]>, +} + +/// A body's display tessellation, ready to evaluate. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct TessPlan { + pub surfaces: Vec, + pub faces: Vec, + /// B-rep edges as polylines (the shared edge nodes), in the body's edge order. + pub edges: Vec>, + /// Per edge: the faces it bounds. + pub edge_faces: Vec>, +} + +impl TessPlan { + /// Evaluate on the CPU into an indexed mesh (vertices are shared within a face, not + /// across faces, so creases keep their own normals). + pub fn evaluate(&self) -> Mesh { + let mut out = Mesh::default(); + for f in &self.faces { + let base = out.positions.len(); + let surface = f.surface.and_then(|i| self.surfaces.get(i as usize)); + for (k, uv) in f.uv.iter().enumerate() { + let (p, n) = match surface { + Some(s) => s.eval(uv[0], uv[1]), + None => (Vec3::ZERO, f.normals.get(k).copied().unwrap_or(Vec3::Z)), + }; + let p = f.fixed.get(k).copied().unwrap_or(p); + out.positions.push(p); + out.normals.push(if f.flip { -n } else { n }); + } + let n = f.uv.len(); + for t in &f.triangles { + if t.iter().all(|&i| (i as usize) < n) { + let at = |i: u32| u32::try_from(base + i as usize).unwrap_or(u32::MAX); + out.triangles.push([at(t[0]), at(t[1]), at(t[2])]); + out.tri_face.push(f.face); + } + } + } + out.edges = self.edges.clone(); + out.edge_faces = self.edge_faces.clone(); + out + } + + /// Vertices across all faces. + pub fn vertex_count(&self) -> usize { + self.faces.iter().map(|f| f.uv.len()).sum() + } +} + +impl EvalCurve { + /// Point and first derivative at `t`; `None` if the curve is malformed or its weight + /// vanishes there. + pub fn eval(&self, t: f64) -> Option<(Vec3, Vec3)> { + let p = self.degree; + let n = self.ctrl.len(); + if p == 0 || p > MAX_SPLINE_DEGREE || n <= p || self.knots.len() != n + p + 1 || !t.is_finite() { + return None; + } + let span = span_of(&self.knots, p, n, t); + let ders = basis_ders(&self.knots, span, p, t)?; + let (mut a, mut da) = ([0.0; 4], [0.0; 4]); + for (k, (n0, n1)) in ders[0].iter().zip(&ders[1]).take(p + 1).enumerate() { + let c = self.ctrl.get(span + k - p)?; + for d in 0..4 { + a[d] += n0 * c[d]; + da[d] += n1 * c[d]; + } + } + rational(a, da) + } +} + +impl EvalSurface { + /// Point and unit normal (`S_u × S_v`, oriented as the surface) at `(u, v)`. A degenerate + /// point (a cone's apex, a sphere's pole) takes its normal from a step inside; anything + /// that still fails gives the origin and +Z rather than a NaN. + pub fn eval(&self, u: f64, v: f64) -> (Vec3, Vec3) { + let Some((p, du, dv)) = self.derivs(u, v) else { return (Vec3::ZERO, Vec3::Z) }; + let n = self.oriented(du.cross(dv)); + if let Some(n) = n.normalized().filter(|_| n.len2() > 1e-20 * (du.len2() * dv.len2()).max(1e-300)) { + return (p, n); + } + // Step into the parameter domain around the degenerate point. + let h = 1e-6 * (1.0 + u.abs().max(v.abs())); + for (su, sv) in [(1.0, 1.0), (-1.0, 1.0), (1.0, -1.0), (-1.0, -1.0)] { + if let Some((_, a, b)) = self.derivs(u + su * h, v + sv * h) + && let Some(n) = self.oriented(a.cross(b)).normalized() + { + return (p, n); + } + } + (p, Vec3::Z) + } + + /// Mirror transforms reverse the cross product of the mapped derivatives. + fn oriented(&self, n: Vec3) -> Vec3 { + match self { + EvalSurface::Revolved { transform: m, .. } if det3(m) < 0.0 => -n, + _ => n, + } + } + + /// Point and first partial derivatives. + pub fn derivs(&self, u: f64, v: f64) -> Option<(Vec3, Vec3, Vec3)> { + if !(u.is_finite() && v.is_finite()) { + return None; + } + match self { + EvalSurface::Plane { origin, du, dv } => Some((*origin + *du * u + *dv * v, *du, *dv)), + EvalSurface::Nurbs { degree, knots, count, ctrl } => nurbs_derivs(*degree, knots, *count, ctrl, u, v), + EvalSurface::Revolved { profile, origin, axis, transform, swap } => { + let (s, a) = if *swap { (v, u) } else { (u, v) }; + let (c, dc) = profile.eval(s)?; + let (sin, cos) = a.sin_cos(); + let rot = |x: Vec3| x * cos + axis.cross(x) * sin + *axis * (axis.dot(x) * (1.0 - cos)); + let r = c - *origin; + let p = *origin + rot(r); + let ds = rot(dc); + // d/dθ of the rotated point: axis × (rotated radius). + let da = axis.cross(rot(r)); + let (p, ds, da) = (affine(transform, p), linear(transform, ds), linear(transform, da)); + Some(if *swap { (p, da, ds) } else { (p, ds, da) }) + } + } + } +} + +fn affine(m: &[[f64; 4]; 3], p: Vec3) -> Vec3 { + linear(m, p) + Vec3::new(m[0][3], m[1][3], m[2][3]) +} + +fn linear(m: &[[f64; 4]; 3], v: Vec3) -> Vec3 { + let row = |r: &[f64; 4]| r[0] * v.x + r[1] * v.y + r[2] * v.z; + Vec3::new(row(&m[0]), row(&m[1]), row(&m[2])) +} + +fn det3(m: &[[f64; 4]; 3]) -> f64 { + m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) +} + +/// The point and derivative of a rational function from its homogeneous value and derivative. +fn rational(a: [f64; 4], da: [f64; 4]) -> Option<(Vec3, Vec3)> { + let w = a[3]; + if !(w.abs() > 1e-300) { + return None; + } + let p = Vec3::new(a[0] / w, a[1] / w, a[2] / w); + let d = (Vec3::new(da[0], da[1], da[2]) - p * da[3]) * (1.0 / w); + (p.is_finite() && d.is_finite()).then_some((p, d)) +} + +fn nurbs_derivs(degree: [usize; 2], knots: &[Vec; 2], count: [usize; 2], ctrl: &[[f64; 4]], u: f64, v: f64) -> Option<(Vec3, Vec3, Vec3)> { + let [pu, pv] = degree; + let [nu, nv] = count; + let ok = |p: usize, n: usize, k: &[f64]| (1..=MAX_SPLINE_DEGREE).contains(&p) && n > p && k.len() == n + p + 1; + if !ok(pu, nu, &knots[0]) || !ok(pv, nv, &knots[1]) || nu.checked_mul(nv)? != ctrl.len() { + return None; + } + let (su, sv) = (span_of(&knots[0], pu, nu, u), span_of(&knots[1], pv, nv, v)); + let (bu, bv) = (basis_ders(&knots[0], su, pu, u)?, basis_ders(&knots[1], sv, pv, v)?); + let (mut s, mut s_u, mut s_v) = ([0.0; 4], [0.0; 4], [0.0; 4]); + for (b, (nv0, nv1)) in bv[0].iter().zip(&bv[1]).take(pv + 1).enumerate() { + for (a, (nu0, nu1)) in bu[0].iter().zip(&bu[1]).take(pu + 1).enumerate() { + let c = ctrl.get((sv + b - pv) * nu + su + a - pu)?; + for d in 0..4 { + s[d] += nu0 * nv0 * c[d]; + s_u[d] += nu1 * nv0 * c[d]; + s_v[d] += nu0 * nv1 * c[d]; + } + } + } + let (p, du) = rational(s, s_u)?; + let (_, dv) = rational(s, s_v)?; + Some((p, du, dv)) +} + +/// The knot span used at `t`: the last non-empty span starting at or before `t`, clamped to +/// the valid spans `p ≤ k < n` so parameters just outside the domain extend the end spans. +fn span_of(knots: &[f64], p: usize, n: usize, t: f64) -> usize { + let mut k = p; + for i in p..n { + if let (Some(&a), Some(&b)) = (knots.get(i), knots.get(i + 1)) + && a < b + { + if a <= t { + k = i; + } else { + break; + } + } + } + k +} + +/// Non-zero basis functions of degree `p` on span `span` and their first derivatives +/// (NURBS Book A2.3 with n = 1). `ders[0][j]` is `N_{span−p+j}`, `ders[1][j]` its derivative. +fn basis_ders(knots: &[f64], span: usize, p: usize, t: f64) -> Option<[[f64; W]; 2]> { + if p == 0 || p >= W { + return None; + } + let kn = |i: usize| knots.get(i).copied(); + // ndu[j][r]: basis values (upper triangle incl. diagonal), knot differences (lower). + let mut ndu = [[0.0f64; W]; W]; + let (mut left, mut right) = ([0.0f64; W], [0.0f64; W]); + ndu[0][0] = 1.0; + for j in 1..=p { + left[j] = t - kn((span + 1).checked_sub(j)?)?; + right[j] = kn(span + j)? - t; + let mut saved = 0.0; + for r in 0..j { + ndu[j][r] = right[r + 1] + left[j - r]; + let temp = if ndu[j][r] != 0.0 { ndu[r][j - 1] / ndu[j][r] } else { 0.0 }; + ndu[r][j] = saved + right[r + 1] * temp; + saved = left[j - r] * temp; + } + ndu[j][j] = saved; + } + let mut ders = [[0.0f64; W]; 2]; + for j in 0..=p { + ders[0][j] = ndu[j][p]; + } + // First derivative: p·(N_{r,p−1}/(u_{r+p}−u_r) − N_{r+1,p−1}/(u_{r+p+1}−u_{r+1})). + for r in 0..=p { + let mut d = 0.0; + if r >= 1 { + let den = ndu[p][r - 1]; + if den != 0.0 { + d += ndu[r - 1][p - 1] / den; + } + } + if r < p { + let den = ndu[p][r]; + if den != 0.0 { + d -= ndu[r][p - 1] / den; + } + } + ders[1][r] = d * p as f64; + } + Some(ders) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: Vec3, b: Vec3, tol: f64) -> bool { + a.dist(b) <= tol + } + + fn line(a: Vec3, b: Vec3) -> EvalCurve { + EvalCurve { degree: 1, knots: vec![0.0, 0.0, 1.0, 1.0], ctrl: vec![[a.x, a.y, a.z, 1.0], [b.x, b.y, b.z, 1.0]] } + } + + const IDENTITY: [[f64; 4]; 3] = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]; + + #[test] + fn plane_point_and_normal() { + let s = EvalSurface::Plane { origin: Vec3::new(1.0, 2.0, 3.0), du: Vec3::X, dv: Vec3::Y }; + let (p, n) = s.eval(2.0, 3.0); + assert!(close(p, Vec3::new(3.0, 5.0, 3.0), 1e-12)); + assert!(close(n, Vec3::Z, 1e-12)); + } + + #[test] + fn revolved_line_is_a_cylinder() { + // A vertical line at radius 2 turned about Z. + let s = EvalSurface::Revolved { + profile: line(Vec3::new(2.0, 0.0, 0.0), Vec3::new(2.0, 0.0, 5.0)), + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: IDENTITY, + swap: false, + }; + for &(u, v) in &[(0.0, 0.0), (0.5, 1.0), (1.0, 4.0), (0.25, 6.0)] { + let (p, n) = s.eval(u, v); + assert!(close(p, Vec3::new(2.0 * v.cos(), 2.0 * v.sin(), 5.0 * u), 1e-12), "{p:?}"); + // S_u = +Z, S_v = tangential: the normal points inward (−radial) for this profile. + assert!(close(n, Vec3::new(-v.cos(), -v.sin(), 0.0), 1e-12), "{n:?}"); + } + } + + #[test] + fn swapped_parameters_reverse_the_normal() { + let mk = |swap| EvalSurface::Revolved { + profile: line(Vec3::new(2.0, 0.0, 0.0), Vec3::new(2.0, 0.0, 5.0)), + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: IDENTITY, + swap, + }; + let (p0, n0) = mk(false).eval(0.3, 1.1); + let (p1, n1) = mk(true).eval(1.1, 0.3); + assert!(close(p0, p1, 1e-12)); + assert!(close(n0, -n1, 1e-12)); + } + + #[test] + fn mirror_transform_keeps_the_normal_outward() { + let mirror = [[-1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]; + let s = EvalSurface::Revolved { + profile: line(Vec3::new(2.0, 0.0, 0.0), Vec3::new(2.0, 0.0, 5.0)), + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: mirror, + swap: false, + }; + let (p, n) = s.eval(0.5, 0.0); + assert!(close(p, Vec3::new(-2.0, 0.0, 2.5), 1e-12)); + // The unmirrored normal at θ = 0 is −X; mirrored it is +X, still pointing to the axis. + assert!(close(n, Vec3::X, 1e-12), "{n:?}"); + } + + #[test] + fn rational_quarter_circle_profile_is_exact() { + // A quarter circle of radius 3 in the XZ plane, turned about Z: a sphere zone. + let w = std::f64::consts::FRAC_1_SQRT_2; + let profile = EvalCurve { + degree: 2, + knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + ctrl: vec![[3.0, 0.0, 0.0, 1.0], [3.0 * w, 0.0, 3.0 * w, w], [0.0, 0.0, 3.0, 1.0]], + }; + let s = EvalSurface::Revolved { profile, origin: Vec3::ZERO, axis: Vec3::Z, transform: IDENTITY, swap: false }; + for i in 0..=10 { + for j in 0..=6 { + let (u, v) = (i as f64 / 10.0, j as f64); + let (p, n) = s.eval(u, v); + assert!((p.len() - 3.0).abs() < 1e-12, "{p:?}"); + // At the pole (u = 1) the normal comes from a step inside: within ~1e-6 rad. + assert!(n.cross(p).len() < 3e-5, "normal {n:?} not radial at {p:?}"); + } + } + } + + #[test] + fn bilinear_nurbs_patch_and_derivatives() { + // Degree (1, 1) over [0, 2] × [0, 1] with corners lifted: a hyperbolic paraboloid. + let ctrl = vec![[0.0, 0.0, 0.0, 1.0], [2.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 1.0], [2.0, 1.0, 1.0, 1.0]]; + let s = EvalSurface::Nurbs { degree: [1, 1], knots: [vec![0.0, 0.0, 2.0, 2.0], vec![0.0, 0.0, 1.0, 1.0]], count: [2, 2], ctrl }; + let Some((p, du, dv)) = s.derivs(1.0, 0.5) else { panic!("evaluates") }; + assert!(close(p, Vec3::new(1.0, 0.5, 0.25), 1e-12), "{p:?}"); + assert!(close(du, Vec3::new(1.0, 0.0, 0.25), 1e-12), "{du:?}"); + assert!(close(dv, Vec3::new(0.0, 1.0, 0.5), 1e-12), "{dv:?}"); + } + + #[test] + fn cubic_nurbs_matches_the_geom_surface() { + // Cross-check the evaluator against `bspline::NurbsSurface::eval` on a weighted patch. + let knots = [vec![0.0, 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, 1.0, 1.0], vec![0.0, 0.0, 0.0, 1.0, 2.0, 2.0, 2.0]]; + let (nu, nv) = (5, 4); + let mut pts = Vec::new(); + let mut wts = Vec::new(); + for j in 0..nv { + for i in 0..nu { + pts.push(Vec3::new(i as f64, j as f64, ((i * 7 + j * 3) % 5) as f64 * 0.3)); + wts.push(1.0 + ((i + 2 * j) % 3) as f64 * 0.25); + } + } + let Ok(reference) = crate::bspline::NurbsSurface::new([3, 2], knots.clone(), &pts, &wts) else { panic!("valid surface") }; + let ctrl = reference.homogeneous().iter().map(|h| [h.x, h.y, h.z, h.w]).collect(); + let s = EvalSurface::Nurbs { degree: [3, 2], knots, count: [nu, nv], ctrl }; + for i in 0..=8 { + for j in 0..=8 { + let (u, v) = (i as f64 / 8.0, j as f64 / 4.0); + let Ok(want) = reference.eval(u, v) else { panic!("in domain") }; + let Some((p, du, dv)) = s.derivs(u, v) else { panic!("evaluates") }; + assert!(close(p, want, 1e-12), "({u}, {v}): {p:?} vs {want:?}"); + // Derivatives against central differences (one-sided at the domain ends). + let h = 1e-6; + let fd = |a: (f64, f64), b: (f64, f64)| { + let (Some(pa), Some(pb)) = (s.derivs(a.0, a.1), s.derivs(b.0, b.1)) else { panic!("evaluates") }; + (pb.0 - pa.0) * (1.0 / (b.0 - a.0 + b.1 - a.1)) + }; + let fu = fd((u.max(h) - h, v), ((u + h).min(1.0), v)); + let fv = fd((u, v.max(h) - h), (u, (v + h).min(2.0))); + assert!(close(du, fu, 1e-4 * (1.0 + du.len())), "du at ({u}, {v}): {du:?} vs {fu:?}"); + assert!(close(dv, fv, 1e-4 * (1.0 + dv.len())), "dv at ({u}, {v}): {dv:?} vs {fv:?}"); + } + } + } + + #[test] + fn hostile_surfaces_do_not_panic() { + let bad = [ + EvalSurface::Nurbs { degree: [9, 1], knots: [vec![0.0; 3], vec![0.0; 3]], count: [1, 1], ctrl: vec![] }, + EvalSurface::Nurbs { + degree: [1, 1], + knots: [vec![0.0, 0.0, 1.0, 1.0], vec![0.0, 0.0, 1.0, 1.0]], + count: [2, 2], + ctrl: vec![[0.0; 4]; 4], + }, + EvalSurface::Revolved { + profile: EvalCurve { degree: 0, knots: vec![], ctrl: vec![] }, + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: IDENTITY, + swap: false, + }, + ]; + for s in &bad { + for &(u, v) in &[(0.5, 0.5), (f64::NAN, 0.0), (1e300, -1e300)] { + let (p, n) = s.eval(u, v); + assert!(p.is_finite() && n.is_finite()); + } + } + } + + #[test] + fn evaluate_uses_fixed_points_and_drops_bad_indices() { + let plan = TessPlan { + surfaces: vec![EvalSurface::Plane { origin: Vec3::ZERO, du: Vec3::X, dv: Vec3::Y }], + faces: vec![FacePlan { + face: 3, + surface: Some(0), + uv: vec![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], + fixed: vec![Vec3::new(0.0, 0.0, 1e-9)], + normals: vec![], + flip: true, + triangles: vec![[0, 1, 2], [0, 1, 7]], + }], + edges: vec![], + edge_faces: vec![], + }; + let m = plan.evaluate(); + assert_eq!(m.positions.len(), 3); + assert_eq!(m.positions.first().copied(), Some(Vec3::new(0.0, 0.0, 1e-9))); + assert_eq!(m.triangles, vec![[0, 1, 2]]); + assert_eq!(m.tri_face, vec![3]); + assert!(m.normals.iter().all(|n| close(*n, -Vec3::Z, 1e-12))); + } +} diff --git a/crates/kernel/Cargo.toml b/crates/kernel/Cargo.toml index fe69d82c..18712825 100644 --- a/crates/kernel/Cargo.toml +++ b/crates/kernel/Cargo.toml @@ -19,3 +19,4 @@ truck-shapeops = { workspace = true } truck-meshalgo = { workspace = true } truck-stepio = { workspace = true } truck-topology = { workspace = true } +spade = { workspace = true } diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index 7e2a3ac7..31d8a0bd 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -1047,6 +1047,18 @@ impl Body { }) } + /// The display tessellation as a plan for the CPU or GPU evaluator: every face's + /// triangulation in its own parameters, with each edge sampled once (chord tolerance + /// `tol`, mm) so faces meet exactly. + pub fn display_plan(&self, tol: f64) -> Result { + if self.is_mesh() { + return Err(KernelError::Invalid("a mesh body has no surfaces to tessellate".into())); + } + let tol = if tol.is_finite() && tol > 0.0 { tol } else { 0.05 }; + let tol = tol.max(self.size() * 1e-6); + guard("tessellate", || crate::tess::plan(&self.solid, tol).map(|(plan, _)| plan)) + } + /// Faces with no triangles in a coarse tessellation (the mesher failed or gave up). pub fn unmeshed_faces(&self) -> usize { if self.is_mesh() { diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 98f48baf..4c6cccc5 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -54,6 +54,7 @@ mod native; mod numbered; mod offset; mod ops; +mod par; mod plane_fit; mod polybool; mod polyhedron; @@ -72,6 +73,7 @@ mod step_out; mod surface_kind; mod surfaces; mod sweep_path; +mod tess; #[cfg(test)] mod thread_tests; mod topo; diff --git a/crates/kernel/src/par.rs b/crates/kernel/src/par.rs new file mode 100644 index 00000000..930a9868 --- /dev/null +++ b/crates/kernel/src/par.rs @@ -0,0 +1,72 @@ +//! Ordered, bounded parallel map over independent items (desktop only; the browser build runs +//! serially). Each item runs under [`crate::guard`]; a worker that cannot be spawned runs its +//! chunk on the calling thread, and a worker that dies reports an internal error for its items. + +use crate::{KernelError, Result, guard}; + +/// Most worker threads. +const MAX_WORKERS: usize = 8; + +pub(crate) fn workers() -> usize { + #[cfg(not(target_arch = "wasm32"))] + { + std::thread::available_parallelism().map_or(1, |n| n.get().min(MAX_WORKERS)) + } + #[cfg(target_arch = "wasm32")] + { + 1 + } +} + +/// `f` over `items`, results in item order. Runs serially below `min_items` items. +pub(crate) fn map(what: &'static str, items: &[T], min_items: usize, f: impl Fn(&T) -> Result + Sync) -> Vec> { + let one = |item: &T| guard(what, || f(item)); + #[cfg(not(target_arch = "wasm32"))] + { + let workers = workers(); + if items.len() >= min_items.max(2) && workers > 1 { + return std::thread::scope(|scope| { + let size = items.len().div_ceil(workers); + let one = &one; + let pending: Vec<_> = items + .chunks(size) + .map(|chunk| { + let handle = + std::thread::Builder::new().name(what.into()).spawn_scoped(scope, move || chunk.iter().map(one).collect::>()).ok(); + (chunk, handle) + }) + .collect(); + let mut out = Vec::with_capacity(items.len()); + for (chunk, handle) in pending { + match handle { + Some(h) => match h.join() { + Ok(results) => out.extend(results), + Err(_) => out.extend(chunk.iter().map(|_| Err(KernelError::Internal(format!("{what}: worker stopped"))))), + }, + None => out.extend(chunk.iter().map(one)), + } + } + out + }); + } + } + items.iter().map(one).collect() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn keeps_order_and_errors() { + let items: Vec = (0..1000).collect(); + let out = map("test", &items, 1, |x| if *x == 500 { Err(KernelError::Invalid("five hundred".into())) } else { Ok(x * 2) }); + assert_eq!(out.len(), 1000); + for (i, r) in out.iter().enumerate() { + match r { + Ok(v) => assert_eq!(*v, 2 * i as u32), + Err(_) => assert_eq!(i, 500), + } + } + } +} diff --git a/crates/kernel/src/tess/cdt.rs b/crates/kernel/src/tess/cdt.rs new file mode 100644 index 00000000..1c317ced --- /dev/null +++ b/crates/kernel/src/tess/cdt.rs @@ -0,0 +1,276 @@ +//! Faces the structured grid cannot take (trimmed faces, holes): a constrained Delaunay +//! triangulation of the boundary loops in (u, v) scaled by the surface's mean metric, with +//! interior seed points, refinement wherever a triangle sags more than the tolerance from the +//! surface, then Ruppert refinement for shape. Constraint edges are never split, so the +//! boundary stays exactly the shared edge nodes. + +use std::collections::HashMap; + +use solidcraft_geom::Vec3; + +use spade::handles::FixedVertexHandle; +use spade::{AngleLimit, ConstrainedDelaunayTriangulation, Point2, RefinementParameters, Triangulation}; + +type Cdt = ConstrainedDelaunayTriangulation>; + +/// Smallest angle Ruppert refinement aims for. +const MIN_ANGLE_DEG: f64 = 20.0; + +/// Rounds of sag refinement. +const SAG_ROUNDS: usize = 16; + +/// The surface, for refining where the triangles leave it. +pub(crate) struct Sag<'a> { + pub eval: &'a dyn Fn(f64, f64) -> Option, + pub tol: f64, + /// Scale of u and v in the triangulation (the surface's mean |S_u|, |S_v|). + pub scale: [f64; 2], +} + +/// How to fill the inside of the loops. +pub(crate) struct Fill<'a> { + /// Interior points (already inside the loops). + pub seeds: &'a [[f64; 2]], + /// Refinement, adding at most this many points in all. + pub refine: Option<(Sag<'a>, usize)>, +} + +/// Triangulate the region bounded by `loops` (ranges of `boundary`, each closed). Boundary +/// vertices keep their indices; added points follow. Triangles are counter-clockwise in +/// (u, v). `None` when a loop cannot be kept as constraints (it crosses itself or another). +pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range], fill: Fill) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + if boundary.iter().flatten().any(|x| !x.is_finite()) { + return None; + } + let scale = match &fill.refine { + Some((sag, _)) if sag.scale.iter().all(|k| k.is_finite() && *k > 0.0) => sag.scale, + _ => [1.0, 1.0], + }; + let to = |p: &[f64; 2]| Point2::new(p[0] * scale[0], p[1] * scale[1]); + let from = |p: Point2| [p.x / scale[0], p.y / scale[1]]; + let mut cdt = Cdt::new(); + let mut index: HashMap = HashMap::new(); + let mut handles: Vec = Vec::with_capacity(boundary.len()); + for (i, p) in boundary.iter().enumerate() { + let h = cdt.insert(to(p)).ok()?; + index.entry(h).or_insert(u32::try_from(i).ok()?); + handles.push(h); + } + for r in loops { + let n = r.len(); + for k in 0..n { + let (a, b) = (*handles.get(r.start + k)?, *handles.get(r.start + (k + 1) % n)?); + if a == b || cdt.exists_constraint(a, b) { + continue; + } + if !cdt.can_add_constraint(a, b) { + return None; + } + cdt.add_constraint(a, b); + } + } + for s in fill.seeds { + if s.iter().all(|x| x.is_finite()) { + let _ = cdt.insert(to(s)); + } + } + if let Some((sag, budget)) = &fill.refine { + let constraints = cdt.num_constraints(); + let added = refine_sag(&mut cdt, sag, *budget, &from); + let params = RefinementParameters::::new() + .keep_constraint_edges() + .exclude_outer_faces(true) + .with_angle_limit(AngleLimit::from_deg(MIN_ANGLE_DEG)) + .with_max_additional_vertices(budget.saturating_sub(added)); + cdt.refine(params); + // A split constraint would put a point on the boundary the neighbour does not have. + if cdt.num_constraints() != constraints { + return None; + } + } + let mut uv: Vec<[f64; 2]> = boundary.to_vec(); + let mut tris = Vec::new(); + let region = region(&cdt); + for face in cdt.inner_faces() { + if !region.get(face.index()).copied().unwrap_or(false) { + continue; + } + let vs = face.vertices(); + let mut t = [0u32; 3]; + for (slot, v) in t.iter_mut().zip(vs) { + *slot = match index.get(&v.fix()) { + Some(&i) => i, + None => { + let i = u32::try_from(uv.len()).ok()?; + uv.push(from(v.position())); + index.insert(v.fix(), i); + i + } + }; + } + tris.push(t); + } + let extra = uv.split_off(boundary.len()); + Some((extra, tris)) +} + +/// Insert the parameters of the surface point under every triangle centre and every +/// non-constraint edge middle that lies more than the tolerance from the triangle, until none +/// does or `budget` points were added. Returns the points added. +fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2) -> [f64; 2]) -> usize { + let mut added = 0; + let mut cache: HashMap> = HashMap::new(); + for _ in 0..SAG_ROUNDS { + let mut todo: Vec> = Vec::new(); + let mut at = |h: FixedVertexHandle, p: Point2| { + *cache.entry(h).or_insert_with(|| { + let q = from(p); + (sag.eval)(q[0], q[1]) + }) + }; + let region = region(cdt); + for face in cdt.inner_faces() { + if !region.get(face.index()).copied().unwrap_or(false) { + continue; + } + let vs = face.vertices(); + let ps = vs.map(|v| v.position()); + let c = Point2::new((ps[0].x + ps[1].x + ps[2].x) / 3.0, (ps[0].y + ps[1].y + ps[2].y) / 3.0); + let (Some(a), Some(b), Some(d)) = (at(vs[0].fix(), ps[0]), at(vs[1].fix(), ps[1]), at(vs[2].fix(), ps[2])) else { continue }; + let off = |q: Point2, chord: Vec3| { + let uv = from(q); + (sag.eval)(uv[0], uv[1]).is_some_and(|s| s.dist(chord) > sag.tol) + }; + if off(c, (a + b + d) * (1.0 / 3.0)) { + todo.push(c); + continue; + } + for (e, (pa, pb)) in face.adjacent_edges().into_iter().zip([(a, b), (b, d), (d, a)]) { + if e.is_constraint_edge() { + continue; + } + let [p, q] = e.positions(); + let m = Point2::new(0.5 * (p.x + q.x), 0.5 * (p.y + q.y)); + if off(m, (pa + pb) * 0.5) { + todo.push(m); + } + } + } + if todo.is_empty() { + break; + } + todo.sort_by(|a, b| a.x.total_cmp(&b.x).then(a.y.total_cmp(&b.y))); + todo.dedup(); + for p in todo { + if added >= budget { + return added; + } + if cdt.insert(p).is_ok() { + added += 1; + } + } + } + added +} + +/// Per face (by index): inside the region the constraint loops bound. A flood fill from the +/// outer face that flips parity each time it crosses a constraint edge (holes are out), linear +/// in the number of triangles. +fn region(cdt: &Cdt) -> Vec { + use spade::handles::{FixedFaceHandle, InnerTag}; + let mut parity: Vec> = vec![None; cdt.num_all_faces()]; + let mut todo: Vec<(FixedFaceHandle, bool)> = Vec::new(); + for e in cdt.convex_hull() { + // Hull edges have the outer face on one side; the other side is an inner face. + for side in [e, e.rev()] { + if let Some(inner) = side.face().as_inner() { + todo.push((inner.fix(), e.is_constraint_edge())); + } + } + } + while let Some((h, odd)) = todo.pop() { + match parity.get_mut(h.index()) { + Some(slot @ None) => *slot = Some(odd), + _ => continue, + } + for e in cdt.face(h).adjacent_edges() { + if let Some(next) = e.rev().face().as_inner() { + todo.push((next.fix(), odd ^ e.is_constraint_edge())); + } + } + } + parity.into_iter().map(|p| p.unwrap_or(false)).collect() +} + +/// The grid points `us` × `vs` inside the region bounded by the loops (even–odd, so holes are +/// out): per grid line, the loops' crossings sorted once, then each point by binary search. +pub(crate) fn inside_grid(us: &[f64], vs: &[f64], polys: &[&[[f64; 2]]]) -> Vec<[f64; 2]> { + let mut out = Vec::new(); + let mut xs: Vec = Vec::new(); + for &v in vs { + xs.clear(); + for poly in polys { + let n = poly.len(); + for k in 0..n { + let (Some(a), Some(b)) = (poly.get(k), poly.get((k + 1) % n)) else { continue }; + if (a[1] > v) != (b[1] > v) { + xs.push(a[0] + (v - a[1]) * (b[0] - a[0]) / (b[1] - a[1])); + } + } + } + xs.sort_by(f64::total_cmp); + for &u in us { + if xs.partition_point(|x| *x < u) % 2 == 1 { + out.push([u, v]); + } + } + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn area(uv: &[[f64; 2]], tris: &[[u32; 3]]) -> f64 { + tris.iter() + .map(|t| { + let (a, b, c) = (uv[t[0] as usize], uv[t[1] as usize], uv[t[2] as usize]); + 0.5 * ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) + }) + .sum() + } + + #[test] + fn square_with_a_hole_keeps_its_boundary() { + let mut b = vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + b.extend([[0.4, 0.4], [0.4, 0.6], [0.6, 0.6], [0.6, 0.4]]); + let loops = [0..4, 4..8]; + let flat = |u: f64, v: f64| Some(Vec3::new(u, v, 0.0)); + for refine in [false, true] { + let refine = refine.then_some((Sag { eval: &flat, tol: 1e-3, scale: [2.0, 1.0] }, 200)); + let label = refine.is_some(); + let Some((extra, tris)) = triangulate(&b, &loops, Fill { seeds: &[[0.2, 0.2], [0.5, 0.5]], refine }) else { panic!("meshes") }; + let mut uv = b.clone(); + uv.extend(extra); + assert!((area(&uv, &tris) - 0.96).abs() < 1e-12, "refined: {label}"); + // No added point on the boundary: boundary segments stay whole. + for r in &loops { + for k in r.clone() { + let next = if k + 1 == r.end { r.start } else { k + 1 }; + let (p, q) = (k as u32, next as u32); + let uses = + tris.iter().filter(|t| (0..3).any(|e| (t[e] == p && t[(e + 1) % 3] == q) || (t[e] == q && t[(e + 1) % 3] == p))).count(); + assert_eq!(uses, 1, "segment {p}-{q}"); + } + } + } + } + + #[test] + fn crossing_loops_are_refused() { + let b = [[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]; + let whole = std::slice::from_ref(&(0..4)); + assert!(triangulate(&b, whole, Fill { seeds: &[], refine: None }).is_none()); + } +} diff --git a/crates/kernel/src/tess/coons.rs b/crates/kernel/src/tess/coons.rs new file mode 100644 index 00000000..a7056600 --- /dev/null +++ b/crates/kernel/src/tess/coons.rs @@ -0,0 +1,235 @@ +//! Structured meshing of a four-cornered face whose sides are not parameter lines (a fillet +//! strip with mitred ends, the wall of a cross hole): a transfinite (Coons) blend of the four +//! sides gives an interior grid in (u, v), joined to the shared edge nodes by zipped strips. +//! A blend that folds (a strongly non-convex boundary) is refused, and the face goes to the +//! constrained Delaunay path. + +use solidcraft_geom::Vec3; + +/// Smallest turn (radians) of the boundary at a corner, in the tolerance metric. +const CORNER_TURN: f64 = 0.6; + +/// One side: boundary vertex indices from one corner to the next, with each node's position +/// along the side (0 to 1, by length in 3D). +struct Side { + nodes: Vec, + at: Vec, +} + +/// The four corners of a single boundary loop, by how sharply the loop turns there in (u, v) +/// scaled by `scale` (one unit per tolerance cell); `None` unless there are exactly four. +fn corners(uv: &[[f64; 2]], scale: [f64; 2]) -> Option<[usize; 4]> { + let n = uv.len(); + if n < 4 { + return None; + } + let p = |i: usize| uv.get(i % n).map(|q| [q[0] * scale[0], q[1] * scale[1]]); + let mut found = Vec::new(); + for i in 0..n { + let (a, b, c) = (p(i + n - 1)?, p(i)?, p(i + 1)?); + let (d0, d1) = ([b[0] - a[0], b[1] - a[1]], [c[0] - b[0], c[1] - b[1]]); + let (l0, l1) = (d0[0].hypot(d0[1]), d1[0].hypot(d1[1])); + // Two copies of a pole meet with no length between them: that is a corner. + if l0 == 0.0 || l1 == 0.0 { + found.push(i); + continue; + } + let turn = (d0[0] * d1[1] - d0[1] * d1[0]).atan2(d0[0] * d1[0] + d0[1] * d1[1]).abs(); + if turn > CORNER_TURN { + found.push(i); + } + } + let [a, b, c, d] = found[..] else { return None }; + Some([a, b, c, d]) +} + +fn side(fixed: &[Vec3], n: usize, from: usize, to: usize) -> Option { + let len = (to + n - from) % n; + if len == 0 { + return None; + } + let nodes: Vec = (0..=len).map(|k| ((from + k) % n) as u32).collect(); + let mut at = vec![0.0]; + let mut total = 0.0; + for w in nodes.windows(2) { + total += fixed.get(w[0] as usize)?.dist(*fixed.get(w[1] as usize)?); + at.push(total); + } + if !(total > 0.0) { + // A pole side: evenly along it. + let m = (nodes.len() - 1).max(1) as f64; + return Some(Side { at: (0..nodes.len()).map(|k| k as f64 / m).collect(), nodes }); + } + Some(Side { nodes, at: at.into_iter().map(|x| x / total).collect() }) +} + +/// The point of a side at position `s` (0 to 1), linear between its nodes. +fn along(side: &Side, uv: &[[f64; 2]], s: f64) -> Option<[f64; 2]> { + let k = side.at.partition_point(|x| *x < s).clamp(1, side.at.len().saturating_sub(1).max(1)); + let (a0, a1) = (*side.at.get(k - 1)?, *side.at.get(k)?); + let (p, q) = (uv.get(*side.nodes.get(k - 1)? as usize)?, uv.get(*side.nodes.get(k)? as usize)?); + let f = if a1 > a0 { ((s - a0) / (a1 - a0)).clamp(0.0, 1.0) } else { 0.0 }; + Some([p[0] + (q[0] - p[0]) * f, p[1] + (q[1] - p[1]) * f]) +} + +/// Interior line positions for a pair of opposite sides: the sides' own when they have the +/// same nodes (so the grid meets both node for node), else evenly as many as the denser side. +fn lines(a: &Side, b: &Side) -> Vec { + let inner = |s: &Side| s.at.get(1..s.at.len().saturating_sub(1)).map(<[f64]>::to_vec).unwrap_or_default(); + if a.at.len() == b.at.len() { + let l: Vec = inner(a).iter().zip(inner(b)).map(|(x, y)| 0.5 * (x + y)).collect(); + if !l.is_empty() { + return l; + } + } + let n = a.at.len().max(b.at.len()).saturating_sub(1).max(2); + (1..n).map(|i| i as f64 / n as f64).collect() +} + +/// Triangulate the loop `uv` (with its fixed 3D nodes) as a Coons patch: boundary vertices +/// keep their indices, the interior grid follows. Triangles counter-clockwise in (u, v). +pub(crate) fn coons(uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points: usize) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + let n = uv.len(); + let [c0, c1, c2, c3] = corners(uv, scale)?; + // Sides in loop order; the bottom (0) and right (1) run with the loop, the top (2) and + // left (3) against it, so s and t both run from corner c0. + let bottom = side(fixed, n, c0, c1)?; + let right = side(fixed, n, c1, c2)?; + let top = reversed(side(fixed, n, c2, c3)?); + let left = reversed(side(fixed, n, c3, c0)?); + let (ss, ts) = (lines(&bottom, &top), lines(&left, &right)); + if ss.len().checked_mul(ts.len())? > max_points { + return None; + } + let corner = |i: usize| uv.get(i).copied(); + let (p00, p10, p11, p01) = (corner(c0)?, corner(c1)?, corner(c2)?, corner(c3)?); + let mut out: Vec<[f64; 2]> = uv.to_vec(); + let base = u32::try_from(n).ok()?; + for &t in &ts { + for &s in &ss { + let (b, tp, l, r) = (along(&bottom, uv, s)?, along(&top, uv, s)?, along(&left, uv, t)?, along(&right, uv, t)?); + let blend = |d: usize| { + (1.0 - t) * b[d] + t * tp[d] + (1.0 - s) * l[d] + s * r[d] + - ((1.0 - s) * (1.0 - t) * p00[d] + s * (1.0 - t) * p10[d] + (1.0 - s) * t * p01[d] + s * t * p11[d]) + }; + out.push([blend(0), blend(1)]); + } + } + let (ns, nt) = (ss.len(), ts.len()); + let g = |i: usize, j: usize| base + (j * ns + i) as u32; + let mut tris = Vec::new(); + for j in 0..nt.saturating_sub(1) { + for i in 0..ns.saturating_sub(1) { + tris.push([g(i, j), g(i + 1, j), g(i + 1, j + 1)]); + tris.push([g(i, j), g(i + 1, j + 1), g(i, j + 1)]); + } + } + // Strips: each side against its grid row or column, by position along the side. In (s, t) + // the bottom and right strips come out counter-clockwise, the top and left clockwise. + let row = |j: usize| (0..ns).map(|i| g(i, j)).collect::>(); + let col = |i: usize| (0..nt).map(|j| g(i, j)).collect::>(); + zip(&bottom, &row(0), &ss, false, &mut tris)?; + zip(&right, &col(ns - 1), &ts, false, &mut tris)?; + zip(&top, &row(nt - 1), &ss, true, &mut tris)?; + zip(&left, &col(0), &ts, true, &mut tris)?; + // The (s, t) square maps to (u, v) keeping or reversing orientation, with the loop: all + // triangles must share one sign, which is then made counter-clockwise. + let area = |t: &[u32; 3]| -> Option { + let (a, b, c) = (out.get(t[0] as usize)?, out.get(t[1] as usize)?, out.get(t[2] as usize)?); + Some((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) + }; + let first = area(tris.first()?)?; + let sign = if first > 0.0 { 1.0 } else { -1.0 }; + for t in &tris { + if !(area(t)? * sign > 0.0) { + return None; + } + } + if sign < 0.0 { + for t in &mut tris { + t.swap(1, 2); + } + } + let extra = out.split_off(n); + Some((extra, tris)) +} + +fn reversed(s: Side) -> Side { + Side { nodes: s.nodes.into_iter().rev().collect(), at: s.at.into_iter().rev().map(|x| 1.0 - x).collect() } +} + +/// Join a side (its nodes at positions `side.at`) to a grid row at positions `inner_at`. +fn zip(side: &Side, inner: &[u32], inner_at: &[f64], flip: bool, tris: &mut Vec<[u32; 3]>) -> Option<()> { + let mut push = |a: u32, b: u32, c: u32| tris.push(if flip { [a, c, b] } else { [a, b, c] }); + let (outer, oat) = (&side.nodes, &side.at); + let (mut a, mut b) = (0usize, 0usize); + while a + 1 < outer.len() || b + 1 < inner.len() { + let advance_outer = if a + 1 >= outer.len() { + false + } else if b + 1 >= inner.len() { + true + } else { + oat.get(a + 1)? <= inner_at.get(b + 1)? + }; + if advance_outer { + push(*outer.get(a)?, *outer.get(a + 1)?, *inner.get(b)?); + a += 1; + } else { + push(*outer.get(a)?, *inner.get(b + 1)?, *inner.get(b)?); + b += 1; + } + } + Some(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn area(uv: &[[f64; 2]], tris: &[[u32; 3]]) -> f64 { + tris.iter() + .map(|t| { + let (a, b, c) = (uv[t[0] as usize], uv[t[1] as usize], uv[t[2] as usize]); + 0.5 * ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) + }) + .sum() + } + + /// A trapezoid (a strip with mitred ends), counter-clockwise, `k` nodes per side. + fn trapezoid(k: usize) -> Vec<[f64; 2]> { + let corners = [[0.0, 0.0], [10.0, 0.0], [8.0, 2.0], [2.0, 2.0]]; + let mut b = Vec::new(); + for w in 0..4 { + let (p, q) = (corners[w], corners[(w + 1) % 4]); + for i in 0..k { + let f = i as f64 / k as f64; + b.push([p[0] + (q[0] - p[0]) * f, p[1] + (q[1] - p[1]) * f]); + } + } + b + } + + #[test] + fn trapezoid_is_tiled_exactly() { + for k in [2, 5] { + for reverse in [false, true] { + let mut b = trapezoid(k); + if reverse { + b.reverse(); + } + let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); + let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000) else { panic!("meshes (k {k})") }; + let mut uv = b.clone(); + uv.extend(extra); + assert!((area(&uv, &tris) - 16.0).abs() < 1e-9, "k {k} reverse {reverse}: {}", area(&uv, &tris)); + } + } + } + + #[test] + fn five_corners_is_not_a_coons_patch() { + let b = [[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [0.5, 1.0], [0.0, 1.0]]; + let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); + assert!(coons(&b, &fixed, [1.0, 1.0], 100).is_none()); + } +} diff --git a/crates/kernel/src/tess/edges.rs b/crates/kernel/src/tess/edges.rs new file mode 100644 index 00000000..6dc20ce5 --- /dev/null +++ b/crates/kernel/src/tess/edges.rs @@ -0,0 +1,97 @@ +//! Edge nodes: every edge of a body sampled once, so the faces on either side share the very +//! same boundary points (watertight by construction). + +use solidcraft_geom::Vec3; +use truck_meshalgo::prelude::*; +use truck_modeling as mt; + +use crate::body::from_p3; + +/// Most nodes on one edge. +pub(crate) const MAX_EDGE_NODES: usize = 20_000; +/// Deepest bisection of one edge segment. +const MAX_DEPTH: u32 = 24; + +/// Nodes of `edge` from its absolute front to its absolute back, every chord within `tol` of +/// the curve (checked at the quarter points, so an S-bend in one segment is not missed). The +/// end nodes are the vertices' own points, so edges meeting at a vertex meet exactly. +pub(crate) fn sample_edge(edge: &mt::Edge, tol: f64) -> Vec { + let front = from_p3(edge.absolute_front().point()); + let back = from_p3(edge.absolute_back().point()); + let curve = edge.curve(); + if matches!(curve, mt::Curve::Line(_)) { + return vec![front, back]; + } + let (t0, t1) = curve.range_tuple(); + if !(t0.is_finite() && t1.is_finite() && t1 > t0) { + return vec![front, back]; + } + // Our evaluator where the curve has a form it knows (no allocation per point); truck's + // for intersection curves. + let fast = super::surface::eval_curve(&curve); + let at = |t: f64| { + match &fast { + Some(c) => c.eval(t).map(|(p, _)| p), + None => Some(from_p3(curve.subs(t))), + } + .filter(|p| p.is_finite()) + }; + // A closed curve's chord between its ends is empty: start from a few segments. + let initial: u32 = 4; + let t_at = |i: u32| t0 + (t1 - t0) * f64::from(i) / f64::from(initial); + let mut out_t = vec![t0]; + // Depth-first bisection, segment by segment, keeping the output ordered. + let mut stack: Vec<(f64, f64, u32)> = (0..initial).rev().map(|i| (t_at(i), t_at(i + 1), 0)).collect(); + while let Some((a, b, depth)) = stack.pop() { + let flat = depth >= MAX_DEPTH || out_t.len() + stack.len() >= MAX_EDGE_NODES || { + let (Some(pa), Some(pb)) = (at(a), at(b)) else { return vec![front, back] }; + [0.25, 0.5, 0.75].iter().all(|f| at(a + (b - a) * f).is_none_or(|q| q.dist_to_segment(pa, pb) <= tol)) + }; + if flat { + out_t.push(b); + } else { + let m = 0.5 * (a + b); + stack.push((m, b, depth + 1)); + stack.push((a, m, depth + 1)); + } + } + let mut pts: Vec = out_t.iter().filter_map(|&t| at(t)).collect(); + if pts.len() < 2 { + return vec![front, back]; + } + if let Some(p) = pts.first_mut() { + *p = front; + } + if let Some(p) = pts.last_mut() { + *p = back; + } + pts +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn circle_edges_are_within_tolerance_and_end_on_vertices() { + let Ok(cyl) = crate::cylinder(Vec3::ZERO, Vec3::Z, 10.0, 5.0) else { panic!("cylinder") }; + for tol in [0.1, 0.01, 0.001] { + for e in cyl.solid.edge_iter() { + let pts = sample_edge(&e, tol); + assert!(pts.len() >= 2); + assert_eq!(pts.first().copied(), Some(from_p3(e.absolute_front().point()))); + assert_eq!(pts.last().copied(), Some(from_p3(e.absolute_back().point()))); + if matches!(e.curve(), mt::Curve::Line(_)) { + assert_eq!(pts.len(), 2); + continue; + } + // Circle of radius 10: every chord's middle sags at most `tol` (plus rounding). + for w in pts.windows(2) { + let m = (w[0] + w[1]) * 0.5; + let r = (m.x * m.x + m.y * m.y).sqrt(); + assert!(10.0 - r <= tol * 1.0001, "sag {} at tol {tol}", 10.0 - r); + } + } + } + } +} diff --git a/crates/kernel/src/tess/grid.rs b/crates/kernel/src/tess/grid.rs new file mode 100644 index 00000000..9502a1a9 --- /dev/null +++ b/crates/kernel/src/tess/grid.rs @@ -0,0 +1,322 @@ +//! Structured meshing of a face whose boundary runs along its parameter box: a grid of +//! interior lines (spaced by the surface's own curvature), joined to the shared edge nodes on +//! each side by a zipped strip (after Cadium's structured and zipped grids: the larger side +//! fans into the smaller one, so the edge nodes are used exactly as they are). + +/// Most grid points of one face. +pub(crate) const MAX_FACE_GRID: usize = 150_000; + +/// The four sides of a parameter box, in counter-clockwise order. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum Side { + VMin, + UMax, + VMax, + UMin, +} + +const SIDES: [Side; 4] = [Side::VMin, Side::UMax, Side::VMax, Side::UMin]; + +impl Side { + /// Bit in a node's side mask. + fn bit(self) -> u8 { + match self { + Side::VMin => 1, + Side::UMax => 2, + Side::VMax => 4, + Side::UMin => 8, + } + } + /// The parameter that runs along this side. + fn along(self) -> usize { + match self { + Side::VMin | Side::VMax => 0, + Side::UMax | Side::UMin => 1, + } + } +} + +/// A face's boundary split into its four sides: for each side (in [`SIDES`] order) the +/// boundary vertex indices along it, ordered by increasing `along` parameter, corners included. +pub(crate) struct Quad { + sides: [Vec; 4], + bbox: [(f64, f64); 2], +} + +/// The boundary as four sides on its parameter box, when it is one: every node on the box, +/// exactly four corner nodes, and the chains between consecutive corners each on one side. +pub(crate) fn quad(boundary: &[[f64; 2]]) -> Option { + let n = boundary.len(); + if n < 4 { + return None; + } + let mut bbox = [(f64::MAX, f64::MIN); 2]; + for p in boundary { + for d in 0..2 { + bbox[d] = (bbox[d].0.min(p[d]), bbox[d].1.max(p[d])); + } + } + let ext = [bbox[0].1 - bbox[0].0, bbox[1].1 - bbox[1].0]; + if !(ext[0] > 0.0 && ext[1] > 0.0) || !(ext[0].is_finite() && ext[1].is_finite()) { + return None; + } + let eps = [ext[0] * 1e-7, ext[1] * 1e-7]; + let mask = |p: &[f64; 2]| { + let mut m = 0u8; + if (p[1] - bbox[1].0).abs() <= eps[1] { + m |= Side::VMin.bit(); + } + if (p[0] - bbox[0].1).abs() <= eps[0] { + m |= Side::UMax.bit(); + } + if (p[1] - bbox[1].1).abs() <= eps[1] { + m |= Side::VMax.bit(); + } + if (p[0] - bbox[0].0).abs() <= eps[0] { + m |= Side::UMin.bit(); + } + m + }; + let masks: Vec = boundary.iter().map(mask).collect(); + if masks.contains(&0) { + return None; + } + let corners: Vec = (0..n).filter(|&i| masks.get(i).is_some_and(|m| m.count_ones() == 2)).collect(); + if corners.len() != 4 { + return None; + } + let mut sides: [Vec; 4] = Default::default(); + let mut seen = 0u8; + for (k, &c0) in corners.iter().enumerate() { + let c1 = *corners.get((k + 1) % 4)?; + let len = (c1 + n - c0) % n; + let chain: Vec = (0..=len).map(|i| (c0 + i) % n).collect(); + let common = chain.iter().fold(0xFu8, |m, &i| m & masks.get(i).copied().unwrap_or(0)); + let side = SIDES.into_iter().find(|s| common & s.bit() != 0)?; + if seen & side.bit() != 0 { + return None; + } + seen |= side.bit(); + let d = side.along(); + let mut chain: Vec = chain.into_iter().map(|i| i as u32).collect(); + let coord = |i: &u32| boundary.get(*i as usize).map_or(0.0, |p| p[d]); + if chain.first().map(coord) > chain.last().map(coord) { + chain.reverse(); + } + // The edge nodes must advance along the side. + if chain.windows(2).any(|w| coord(&w[0]) >= coord(&w[1])) && chain.len() > 2 { + return None; + } + let slot = SIDES.iter().position(|s| *s == side)?; + *sides.get_mut(slot)? = chain; + } + Some(Quad { sides, bbox }) +} + +/// Lines along direction `d` taken from the boundary, when the two sides running along `d` +/// carry the same parameters (as opposite iso-lines of one surface sampled alike do): the grid +/// then meets both sides node for node. Interior values only. +pub(crate) fn aligned_lines(boundary: &[[f64; 2]], q: &Quad, d: usize) -> Option> { + let (a, b) = if d == 0 { (q.sides.first()?, q.sides.get(2)?) } else { (q.sides.get(1)?, q.sides.get(3)?) }; + if a.len() != b.len() || a.len() < 3 { + return None; + } + let (lo, hi) = q.bbox.get(d).copied()?; + let eps = (hi - lo) * 1e-9; + let mut lines = Vec::with_capacity(a.len()); + for (i, j) in a.iter().zip(b) { + let (x, y) = (boundary.get(*i as usize)?[d], boundary.get(*j as usize)?[d]); + if (x - y).abs() > eps { + return None; + } + if x > lo + eps && x < hi - eps { + lines.push(x); + } + } + Some(lines) +} + +/// A face straight across direction `1 − d` (no interior lines there, as a cylinder wall along +/// its axis): the two sides running along `d` zipped straight to each other, when the other two +/// sides are single segments. Triangles counter-clockwise, or `None`. +pub(crate) fn strip(boundary: &[[f64; 2]], q: &Quad, d: usize) -> Option> { + // Sides along u are VMin (low) and VMax (high); along v, UMax (high) and UMin (low). + let (low, high, cross) = if d == 0 { (0, 2, [1, 3]) } else { (3, 1, [0, 2]) }; + if cross.iter().any(|&k| q.sides.get(k).is_none_or(|c| c.len() != 2)) { + return None; + } + let mut tris = Vec::new(); + zip(boundary, q.sides.get(low)?, q.sides.get(high)?, d, d == 1, &mut tris)?; + for t in &tris { + if !(signed_area(boundary, *t)? > 0.0) { + return None; + } + } + Some(tris) +} + +/// The interior lines of a structured face: `us`, `vs` strictly inside the box. Each +/// direction gets at least one line so every side has an inner chain to zip to. +pub(crate) fn interior(lines: Vec, (a, b): (f64, f64)) -> Vec { + let inner: Vec = lines.into_iter().filter(|x| *x > a && *x < b).collect(); + if inner.is_empty() { vec![0.5 * (a + b)] } else { inner } +} + +/// Triangulate a [`Quad`] face: boundary vertices are `0..boundary.len()`; the interior grid +/// (`us` × `vs`, strictly inside the box) follows, row by row. Returns the interior +/// parameters and the triangles (counter-clockwise in (u, v)), or `None` if a strip would +/// fold over. +pub(crate) fn structured(boundary: &[[f64; 2]], q: &Quad, us: &[f64], vs: &[f64]) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + let (nu, nv) = (us.len(), vs.len()); + if nu == 0 || nv == 0 || nu.checked_mul(nv)? > MAX_FACE_GRID { + return None; + } + let n = u32::try_from(boundary.len()).ok()?; + let g = |i: usize, j: usize| n + (j * nu + i) as u32; + let mut uv: Vec<[f64; 2]> = boundary.to_vec(); + for &v in vs { + for &u in us { + uv.push([u, v]); + } + } + let mut tris: Vec<[u32; 3]> = Vec::with_capacity(2 * (nu + 1) * (nv + 1) + 2 * boundary.len()); + for j in 0..nv.saturating_sub(1) { + for i in 0..nu.saturating_sub(1) { + let (a, b, c, d) = (g(i, j), g(i + 1, j), g(i + 1, j + 1), g(i, j + 1)); + tris.push([a, b, c]); + tris.push([a, c, d]); + } + } + for (slot, side) in SIDES.iter().enumerate() { + let outer = q.sides.get(slot)?; + let inner: Vec = match side { + Side::VMin => (0..nu).map(|i| g(i, 0)).collect(), + Side::VMax => (0..nu).map(|i| g(i, nv - 1)).collect(), + Side::UMin => (0..nv).map(|j| g(0, j)).collect(), + Side::UMax => (0..nv).map(|j| g(nu - 1, j)).collect(), + }; + // The strip lies inside the box from its side: below the top and left of the right + // side its triangles come out clockwise. + let flip = matches!(side, Side::VMax | Side::UMin); + zip(&uv, outer, &inner, side.along(), flip, &mut tris)?; + } + // Every triangle counter-clockwise: one that is not has folded over. + for t in &tris { + if !(signed_area(&uv, *t)? > 0.0) { + return None; + } + } + let extra = uv.split_off(boundary.len()); + Some((extra, tris)) +} + +fn signed_area(uv: &[[f64; 2]], t: [u32; 3]) -> Option { + let (a, b, c) = (uv.get(t[0] as usize)?, uv.get(t[1] as usize)?, uv.get(t[2] as usize)?); + Some((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) +} + +/// Join a side's outer chain (edge nodes) to the inner chain (grid line), both ordered along +/// direction `d`: advance whichever chain's next node comes first along the side, so nodes at +/// the same parameter on both chains are joined straight across. +fn zip(uv: &[[f64; 2]], outer: &[u32], inner: &[u32], d: usize, flip: bool, tris: &mut Vec<[u32; 3]>) -> Option<()> { + let mut push = |a: u32, b: u32, c: u32| tris.push(if flip { [a, c, b] } else { [a, b, c] }); + let coord = |chain: &[u32], k: usize| chain.get(k).and_then(|i| uv.get(*i as usize)).map(|p| p[d]); + let (mut a, mut b) = (0usize, 0usize); + while a + 1 < outer.len() || b + 1 < inner.len() { + let advance_outer = if a + 1 >= outer.len() { + false + } else if b + 1 >= inner.len() { + true + } else { + coord(outer, a + 1)? <= coord(inner, b + 1)? + }; + if advance_outer { + push(*outer.get(a)?, *outer.get(a + 1)?, *inner.get(b)?); + a += 1; + } else { + push(*outer.get(a)?, *inner.get(b + 1)?, *inner.get(b)?); + b += 1; + } + } + Some(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A unit square boundary with `k` nodes per side, counter-clockwise from (0, 0). + fn square(k: usize) -> Vec<[f64; 2]> { + let mut b = Vec::new(); + for i in 0..k { + b.push([i as f64 / k as f64, 0.0]); + } + for i in 0..k { + b.push([1.0, i as f64 / k as f64]); + } + for i in 0..k { + b.push([1.0 - i as f64 / k as f64, 1.0]); + } + for i in 0..k { + b.push([0.0, 1.0 - i as f64 / k as f64]); + } + b + } + + fn area(uv: &[[f64; 2]], tris: &[[u32; 3]]) -> f64 { + tris.iter().map(|t| signed_area(uv, *t).unwrap_or(f64::NAN) * 0.5).sum() + } + + #[test] + fn grid_and_strips_tile_the_box_exactly() { + for (k, us, vs) in [(1, vec![0.5], vec![0.5]), (5, vec![0.3, 0.6], vec![0.2, 0.4, 0.8]), (13, vec![0.5], vec![0.1, 0.9])] { + let b = square(k); + let Some(q) = quad(&b) else { panic!("a quad") }; + let Some((extra, tris)) = structured(&b, &q, &us, &vs) else { panic!("meshes") }; + let mut uv = b.clone(); + uv.extend(extra); + assert!((area(&uv, &tris) - 1.0).abs() < 1e-12, "k = {k}"); + assert!(tris.iter().all(|t| signed_area(&uv, *t).is_some_and(|a| a > 0.0))); + // Every boundary segment is used by exactly one triangle (no node skipped). + let n = b.len() as u32; + for i in 0..n { + let (p, q) = (i, (i + 1) % n); + let uses = tris.iter().filter(|t| (0..3).any(|e| (t[e] == p && t[(e + 1) % 3] == q) || (t[e] == q && t[(e + 1) % 3] == p))).count(); + assert_eq!(uses, 1, "boundary segment {p}-{q}"); + } + } + } + + #[test] + fn a_pole_side_is_a_two_node_chain() { + // A triangle-like patch whose top side is a pole (two copies of one point in 3D). + let b = vec![[0.0, 0.0], [0.5, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + let Some(q) = quad(&b) else { panic!("a quad") }; + let Some((_, tris)) = structured(&b, &q, &[0.5], &[0.5]) else { panic!("meshes") }; + assert!(!tris.is_empty()); + } + + #[test] + fn straight_faces_are_one_strip() { + // Five nodes along the bottom, three along the top, single-segment ends. + let b = vec![[0.0, 0.0], [0.25, 0.0], [0.5, 0.0], [0.75, 0.0], [1.0, 0.0], [1.0, 1.0], [0.5, 1.0], [0.0, 1.0]]; + let Some(q) = quad(&b) else { panic!("a quad") }; + let Some(tris) = strip(&b, &q, 0) else { panic!("a strip") }; + assert_eq!(tris.len(), 6); + assert!((area(&b, &tris) - 1.0).abs() < 1e-12); + // Transposed: the long sides run along v. + let t: Vec<[f64; 2]> = b.iter().map(|p| [p[1], p[0]]).rev().collect(); + let Some(q) = quad(&t) else { panic!("a quad") }; + let Some(tris) = strip(&t, &q, 1) else { panic!("a strip") }; + assert!((area(&t, &tris) - 1.0).abs() < 1e-12); + assert!(strip(&b, &q, 0).is_none() || strip(&t, &q, 0).is_none()); + } + + #[test] + fn not_a_quad() { + // A hole-like loop inside the box, and a five-cornered loop. + assert!(quad(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]).is_none()); + let pentagon = [[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [0.5, 1.0], [0.0, 1.0]]; + assert!(quad(&pentagon).is_none()); + } +} diff --git a/crates/kernel/src/tess/mod.rs b/crates/kernel/src/tess/mod.rs new file mode 100644 index 00000000..06509a1f --- /dev/null +++ b/crates/kernel/src/tess/mod.rs @@ -0,0 +1,307 @@ +//! Display tessellation: each face's triangulation in its own parameters, with every edge +//! sampled once so neighbouring faces share their boundary points exactly, handed over as a +//! [`TessPlan`] for the CPU or GPU evaluator. +//! +//! Per face, cheapest first: +//! - planar faces: a constrained Delaunay triangulation of the edge nodes alone (a plane needs +//! no interior points); +//! - faces whose boundary runs along their parameter box (cylinder halves, sphere and torus +//! patches, fillets): a structured grid zipped to the edge nodes; +//! - four-cornered faces whose sides are not parameter lines (mitred fillet strips, cross +//! hole walls): a Coons blend of the four sides, zipped the same way; +//! - anything else (trimmed faces, holes): constrained Delaunay with interior seeds and +//! Ruppert refinement. +//! +//! A face the planner cannot handle (its surface has no evaluable form, a boundary that does +//! not project or crosses itself) fails the whole plan, and the caller meshes the body the +//! old way. + +pub(crate) mod cdt; +pub(crate) mod coons; +pub(crate) mod edges; +pub(crate) mod grid; +pub(crate) mod surface; +pub(crate) mod uv; + +use std::collections::HashMap; + +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::{EvalSurface, FacePlan, TessPlan}; +use truck_modeling as mt; + +use crate::body::{Solid, edge_order}; +use crate::{KernelError, Result}; + +/// How the faces of a plan were meshed. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct PlanStats { + pub planar: usize, + pub structured: usize, + pub delaunay: usize, +} + +/// Most lines a curved face's grid gets in one direction. +const MAX_LINES: usize = 1024; +/// Most points Ruppert refinement may add to one face, beyond its seeds. +const MAX_REFINE: usize = 20_000; + +fn unsupported(what: impl std::fmt::Display) -> KernelError { + KernelError::Failed(format!("not supported yet: tessellation: {what}")) +} + +/// Plan the display tessellation of `solid` with chord tolerance `tol`. +pub(crate) fn plan(solid: &Solid, tol: f64) -> Result<(TessPlan, PlanStats)> { + let faces: Vec<&mt::Face> = solid.face_iter().collect(); + let (edges, edge_faces) = edge_order(faces.iter().copied()); + // Every edge once, in parallel; faces then read their nodes. + let sampled = crate::par::map("edge nodes", &edges, 64, |e| Ok(edges::sample_edge(e, tol))); + let mut nodes: HashMap<_, Vec> = HashMap::with_capacity(edges.len()); + for (e, pts) in edges.iter().zip(sampled) { + nodes.insert(e.id(), pts?); + } + let oriented = |e: &mt::Edge| -> Option> { + let mut pts = nodes.get(&e.id())?.clone(); + if !e.orientation() { + pts.reverse(); + } + Some(pts) + }; + let indexed: Vec<(usize, &mt::Face)> = faces.iter().copied().enumerate().collect(); + let planned = crate::par::map("face tessellation", &indexed, 8, |(fi, face)| { + let fi = *fi; + let surface = face.surface(); + let eval = surface::eval_surface(&surface).ok_or_else(|| unsupported(format!("face {fi}: surface not evaluable")))?; + // Boundary nodes, loop by loop (each edge's last node is the next edge's first). + let mut loops3: Vec> = Vec::new(); + for wire in face.absolute_boundaries() { + let mut pts = Vec::new(); + for e in wire.edge_iter() { + let mut ep = oriented(e).ok_or_else(|| unsupported("edge without nodes"))?; + ep.pop(); + pts.extend(ep); + } + if pts.len() >= 2 { + loops3.push(pts); + } + } + let mut stats = PlanStats::default(); + let fp = plan_face(fi, &eval, &loops3, tol, &mut stats)?; + let fp = FacePlan { + flip: !face.orientation(), + triangles: if face.orientation() { fp.triangles } else { fp.triangles.into_iter().map(|[a, b, c]| [a, c, b]).collect() }, + ..fp + }; + Ok((fp, eval, stats)) + }); + let mut plan = TessPlan::default(); + let mut stats = PlanStats::default(); + for r in planned { + let (fp, eval, s) = r?; + let surface = u32::try_from(plan.surfaces.len()).map_err(|_| unsupported("too many faces"))?; + plan.surfaces.push(eval); + plan.faces.push(FacePlan { surface: Some(surface), ..fp }); + stats.planar += s.planar; + stats.structured += s.structured; + stats.delaunay += s.delaunay; + } + plan.edges = edges.iter().filter_map(oriented).collect(); + plan.edge_faces = edge_faces.into_iter().map(|fs| fs.into_iter().map(|f| u32::try_from(f).unwrap_or(u32::MAX)).collect()).collect(); + Ok((plan, stats)) +} + +/// One face: boundary loops in 3D (shared edge nodes) to parameters and triangles. +fn plan_face(fi: usize, s: &EvalSurface, loops3: &[Vec], tol: f64, stats: &mut PlanStats) -> Result { + let dom = uv::domain(s); + let accept = (tol * 10.0).max(1e-6); + let fail = |why: &str| unsupported(format!("face {fi}: {why}")); + // Project every loop, unwrapping periodic parameters along it and splitting poles. + let mut fixed: Vec = Vec::new(); + let mut uvs: Vec<[f64; 2]> = Vec::new(); + let mut ranges: Vec> = Vec::new(); + for pts in loops3 { + let mut luv: Vec<[f64; 2]> = Vec::with_capacity(pts.len()); + let mut prev: Option<[f64; 2]> = None; + for p in pts { + let mut q = uv::project(s, &dom, *p, prev, accept).ok_or_else(|| fail("boundary does not lie on the surface"))?; + if let Some(pr) = prev { + for d in 0..2 { + q[d] = uv::unwrap(q[d], pr[d], dom.period[d]); + } + } + luv.push(q); + prev = Some(q); + } + // A loop that winds once around a periodic direction does not close in (u, v). + if let (Some(first), Some(last)) = (luv.first(), luv.last()) { + for d in 0..2 { + if let Some(per) = dom.period[d] + && (uv::unwrap(first[d], last[d], Some(per)) - first[d]).abs() > per * 0.25 + { + return Err(fail("boundary winds around the surface")); + } + } + } + let start = uvs.len(); + let n = luv.len(); + for (k, (q, p)) in luv.iter().zip(pts).enumerate() { + // At a pole the parameter along the degenerate direction is free: one copy with + // the previous node's value, one with the next node's. + let (prev, next) = (luv.get((k + n - 1) % n).copied(), luv.get((k + 1) % n).copied()); + let pole = (0..2).find(|&d| uv::degenerate(s, *q, d)); + match (pole, prev, next) { + (Some(d), Some(a), Some(b)) => { + let mut qa = *q; + let mut qb = *q; + qa[d] = a[d]; + qb[d] = b[d]; + uvs.push(qa); + fixed.push(*p); + if (qa[d] - qb[d]).abs() > 1e-12 { + uvs.push(qb); + fixed.push(*p); + } + } + _ => { + uvs.push(*q); + fixed.push(*p); + } + } + } + ranges.push(start..uvs.len()); + } + if uvs.len() < 3 { + return Err(fail("degenerate boundary")); + } + let planar = matches!(s, EvalSurface::Plane { .. }); + let (extra, triangles) = if planar { + stats.planar += 1; + cdt::triangulate(&uvs, &ranges, cdt::Fill { seeds: &[], refine: None }).ok_or_else(|| fail("boundary crosses itself"))? + } else { + let mut bbox = [(f64::MAX, f64::MIN); 2]; + for p in &uvs { + for d in 0..2 { + bbox[d] = (bbox[d].0.min(p[d]), bbox[d].1.max(p[d])); + } + } + let mut us = uv::divide(s, 0, bbox[0], bbox[1], tol, MAX_LINES); + let mut vs = uv::divide(s, 1, bbox[1], bbox[0], tol, MAX_LINES); + thin(&mut us, &mut vs, grid::MAX_FACE_GRID); + let structured = match (ranges.len(), grid::quad(&uvs)) { + (1, Some(q)) => { + // Interior lines: the boundary's own where opposite sides agree, else the + // tolerance division. + let lines = |d: usize, div: &Vec| { + grid::aligned_lines(&uvs, &q, d).unwrap_or_else(|| div.iter().copied().filter(|x| *x > bbox[d].0 && *x < bbox[d].1).collect()) + }; + let (lu, lv) = (lines(0, &us), lines(1, &vs)); + let strip = match (lu.is_empty(), lv.is_empty()) { + (false, true) => grid::strip(&uvs, &q, 0), + (true, false) => grid::strip(&uvs, &q, 1), + _ => None, + }; + strip.map(|t| (Vec::new(), t)).or_else(|| grid::structured(&uvs, &q, &grid::interior(lu, bbox[0]), &grid::interior(lv, bbox[1]))) + } + _ => None, + }; + let cells = |lines: &Vec, (a, b): (f64, f64)| (lines.len().saturating_sub(1).max(1) as f64) / (b - a).max(f64::MIN_POSITIVE); + let scale = [cells(&us, bbox[0]), cells(&vs, bbox[1])]; + let structured = match structured { + None if ranges.len() == 1 => coons::coons(&uvs, &fixed, scale, grid::MAX_FACE_GRID), + other => other, + }; + match structured { + Some(r) => { + stats.structured += 1; + r + } + None => { + stats.delaunay += 1; + let polys: Vec<&[[f64; 2]]> = ranges.iter().filter_map(|r| uvs.get(r.clone())).collect(); + let inner = |lines: &Vec, (a, b): (f64, f64)| lines.iter().copied().filter(|x| *x > a && *x < b).collect::>(); + let seeds = cdt::inside_grid(&inner(&us, bbox[0]), &inner(&vs, bbox[1]), &polys); + let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, _, _)| p); + // One unit per tolerance cell in each direction: a straight direction (one + // cell) keeps long triangles, Ruppert shapes them in that metric. + let sag = cdt::Sag { eval: &eval, tol, scale }; + let fill = cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)) }; + cdt::triangulate(&uvs, &ranges, fill).ok_or_else(|| fail("boundary crosses itself"))? + } + } + }; + // At a pole the boundary has two copies of one node: a triangle on both copies has two + // corners at one point and no area (its other two sides are one segment both ways). + let pole_pair = |i: u32, j: u32| i != j && fixed.get(i as usize).zip(fixed.get(j as usize)).is_some_and(|(p, q)| p == q); + let mut triangles = triangles; + triangles.retain(|t| !(pole_pair(t[0], t[1]) || pole_pair(t[1], t[2]) || pole_pair(t[2], t[0]))); + if !closed_boundary(&ranges, &fixed, &triangles) { + return Err(fail("triangulation does not follow the boundary")); + } + let mut all = uvs; + all.extend(extra); + Ok(FacePlan { face: u32::try_from(fi).unwrap_or(u32::MAX), surface: None, uv: all, fixed, normals: Vec::new(), flip: false, triangles }) +} + +/// Every boundary segment (consecutive nodes of a loop) is the side of exactly one triangle: +/// the face covers its boundary with nothing missing or doubled, so it meets its neighbours +/// exactly. Pole copies of one node (consecutive in a loop, at one point) are one vertex: the +/// segment between them has no length, and a side ending on either copy is the same side. +fn closed_boundary(ranges: &[std::ops::Range], fixed: &[Vec3], triangles: &[[u32; 3]]) -> bool { + let mut one: Vec = (0..fixed.len() as u32).collect(); + let root = |one: &[u32], mut i: u32| { + while let Some(&j) = one.get(i as usize).filter(|j| **j != i) { + i = j; + } + i + }; + let segment = |r: &std::ops::Range, k: usize| ((r.start + k) as u32, (r.start + (k + 1) % r.len()) as u32); + for r in ranges { + for k in 0..r.len() { + let (a, b) = segment(r, k); + if fixed.get(a as usize).zip(fixed.get(b as usize)).is_some_and(|(p, q)| p == q) { + let (ra, rb) = (root(&one, a), root(&one, b)); + if let Some(x) = one.get_mut(rb.max(ra) as usize) { + *x = ra.min(rb); + } + } + } + } + let side = |a: u32, b: u32| { + let (a, b) = (root(&one, a), root(&one, b)); + (a != b).then_some((a.min(b), a.max(b))) + }; + let mut sides: HashMap<(u32, u32), u32> = HashMap::with_capacity(triangles.len() * 3); + for t in triangles { + for e in 0..3 { + if let Some(k) = side(t[e], t[(e + 1) % 3]) { + *sides.entry(k).or_default() += 1; + } + } + } + ranges.iter().all(|r| { + (0..r.len()).all(|k| { + let (a, b) = segment(r, k); + side(a, b).is_none_or(|k| sides.get(&k).copied() == Some(1)) + }) + }) +} + +/// Keep at most `cap` grid points by taking every k-th line in both directions (ends kept). +fn thin(us: &mut Vec, vs: &mut Vec, cap: usize) { + let total = us.len().saturating_mul(vs.len()); + if total <= cap { + return; + } + let f = (total as f64 / cap.max(4) as f64).sqrt(); + let pick = |d: &mut Vec| { + let n = ((d.len() as f64 / f).ceil() as usize).max(2); + if d.len() > n { + let len = d.len(); + *d = (0..n).filter_map(|i| d.get(i * (len - 1) / (n - 1)).copied()).collect(); + } + }; + pick(us); + pick(vs); +} + +#[cfg(test)] +mod tests; diff --git a/crates/kernel/src/tess/surface.rs b/crates/kernel/src/tess/surface.rs new file mode 100644 index 00000000..f263e343 --- /dev/null +++ b/crates/kernel/src/tess/surface.rs @@ -0,0 +1,156 @@ +//! Truck surfaces in the evaluator's form ([`EvalSurface`]), which the CPU and GPU evaluators +//! share. Surfaces the evaluators cannot express (a revolved intersection curve) or that are +//! too large for them give `None`; their faces are evaluated here instead. + +use solidcraft_geom::Vec3; +use solidcraft_geom::bspline::MAX_SPLINE_DEGREE; +use solidcraft_geom::tess::{EvalCurve, EvalSurface, MAX_EVAL_KNOTS, MAX_EVAL_POINTS}; +use truck_meshalgo::prelude::*; +use truck_modeling as mt; + +use crate::body::from_p3; + +/// `surface` in the evaluator's form, parameterised exactly as truck parameterises it. +pub(crate) fn eval_surface(surface: &mt::Surface) -> Option { + let s = match surface { + mt::Surface::Plane(p) => { + EvalSurface::Plane { origin: from_p3(p.subs(0.0, 0.0)), du: from_v3(p.uder(0.0, 0.0)), dv: from_v3(p.vder(0.0, 0.0)) } + } + mt::Surface::BSplineSurface(b) => { + let ctrl = b.control_points().iter().map(|row| row.iter().map(|p| [p.x, p.y, p.z, 1.0]).collect()).collect(); + nurbs(b.degrees(), b.knot_vecs(), ctrl)? + } + mt::Surface::NurbsSurface(n) => { + let ctrl = n.control_points().iter().map(|row| row.iter().map(|h| [h.x, h.y, h.z, h.w]).collect()).collect(); + nurbs((n.udegree(), n.vdegree()), n.knot_vecs(), ctrl)? + } + mt::Surface::RevolutedCurve(proc) => { + let rev = proc.entity(); + let profile = eval_curve(rev.entity_curve())?; + let axis = from_v3(rev.axis()).normalized()?; + let m = proc.transform(); + // cgmath is column-major: m[column][row]. + let transform = [0, 1, 2].map(|r| [m[0][r], m[1][r], m[2][r], m[3][r]]); + if transform.iter().flatten().any(|x| !x.is_finite()) { + return None; + } + EvalSurface::Revolved { profile, origin: from_p3(rev.origin()), axis, transform, swap: !proc.orientation() } + } + }; + Some(s) +} + +fn from_v3(v: mt::Vector3) -> Vec3 { + Vec3::new(v.x, v.y, v.z) +} + +fn knots(k: &mt::KnotVec) -> Option> { + let v: Vec = k.iter().copied().collect(); + (v.len() <= MAX_EVAL_KNOTS && v.iter().all(|x| x.is_finite())).then_some(v) +} + +/// Truck's `control_points[i][j]` (i along u) as the evaluator's u-fastest array. +fn nurbs((pu, pv): (usize, usize), kv: &(mt::KnotVec, mt::KnotVec), rows: Vec>) -> Option { + let nu = rows.len(); + let nv = rows.first()?.len(); + if pu == 0 || pv == 0 || pu > MAX_SPLINE_DEGREE || pv > MAX_SPLINE_DEGREE || nu.checked_mul(nv)? > MAX_EVAL_POINTS { + return None; + } + let mut ctrl = Vec::with_capacity(nu * nv); + for j in 0..nv { + for row in &rows { + ctrl.push(<[[f64; 4]]>::get(row, j).copied()?); + } + } + if ctrl.iter().flatten().any(|x| !x.is_finite()) { + return None; + } + let knots = [knots(&kv.0)?, knots(&kv.1)?]; + (knots[0].len() == nu + pu + 1 && knots[1].len() == nv + pv + 1).then_some(EvalSurface::Nurbs { degree: [pu, pv], knots, count: [nu, nv], ctrl }) +} + +/// A profile curve in the evaluator's form (a line extends for ever, as truck's does). +pub(crate) fn eval_curve(c: &mt::Curve) -> Option { + let (degree, kv, ctrl): (usize, Vec, Vec<[f64; 4]>) = match c { + mt::Curve::Line(l) => (1, vec![0.0, 0.0, 1.0, 1.0], vec![[l.0.x, l.0.y, l.0.z, 1.0], [l.1.x, l.1.y, l.1.z, 1.0]]), + mt::Curve::BSplineCurve(b) => (b.degree(), knots(b.knot_vec())?, b.control_points().iter().map(|p| [p.x, p.y, p.z, 1.0]).collect()), + mt::Curve::NurbsCurve(n) => (n.degree(), knots(n.knot_vec())?, n.control_points().iter().map(|h| [h.x, h.y, h.z, h.w]).collect()), + mt::Curve::IntersectionCurve(_) => return None, + }; + let ok = (1..=MAX_SPLINE_DEGREE).contains(°ree) + && ctrl.len() <= MAX_EVAL_POINTS + && kv.len() == ctrl.len() + degree + 1 + && ctrl.iter().flatten().all(|x| x.is_finite()); + ok.then_some(EvalCurve { degree, knots: kv, ctrl }) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::body::Solid; + + /// Every face of `solid` evaluates as truck evaluates it, over its own parameter box. + fn assert_matches_truck(solid: &Solid) -> [usize; 4] { + let mut kinds = [0; 4]; + for face in solid.face_iter() { + let surface = face.surface(); + let kind = match surface { + mt::Surface::Plane(_) => 0, + mt::Surface::BSplineSurface(_) => 1, + mt::Surface::NurbsSurface(_) => 2, + mt::Surface::RevolutedCurve(_) => 3, + }; + let Some(eval) = eval_surface(&surface) else { panic!("not evaluable: {surface:?}") }; + // Sample the face's parameter box from its boundary vertices. + let pts: Vec<(f64, f64)> = face.vertex_iter().filter_map(|v| surface.search_parameter(v.point(), None, 100)).collect(); + if pts.is_empty() { + continue; + } + let (u0, u1) = pts.iter().fold((f64::MAX, f64::MIN), |(a, b), p| (a.min(p.0), b.max(p.0))); + let (v0, v1) = pts.iter().fold((f64::MAX, f64::MIN), |(a, b), p| (a.min(p.1), b.max(p.1))); + for i in 0..=6 { + for j in 0..=6 { + let (u, v) = (u0 + (u1 - u0) * i as f64 / 6.0, v0 + (v1 - v0) * j as f64 / 6.0); + let (p, n) = eval.eval(u, v); + let want = from_p3(surface.subs(u, v)); + assert!(p.dist(want) < 1e-9 * (1.0 + want.len()), "point at ({u}, {v}): {p:?} vs {want:?}"); + let wn = from_v3(surface.normal(u, v)); + if wn.len() > 0.5 { + assert!(n.dist(wn) < 1e-6, "normal at ({u}, {v}): {n:?} vs {wn:?} ({surface:?})"); + } + kinds[kind] += 1; + } + } + } + kinds + } + + fn solid(b: crate::Body) -> Solid { + (*b.solid).clone() + } + + #[test] + fn primitives_evaluate_like_truck() { + let Ok(bx) = crate::box_solid(Vec3::new(-1.0, -2.0, 0.0), Vec3::new(3.0, 4.0, 5.0)) else { panic!("box") }; + let Ok(cyl) = crate::cylinder(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 1.0, 1.0), 2.5, 7.0) else { panic!("cylinder") }; + let Ok(sph) = crate::sphere(Vec3::new(0.0, 0.0, 1.0), 4.0) else { panic!("sphere") }; + let Ok(tor) = crate::torus(Vec3::ZERO, 6.0, 1.5) else { panic!("torus") }; + let mut kinds = [0; 4]; + for b in [bx, cyl, sph, tor] { + for (k, n) in kinds.iter_mut().zip(assert_matches_truck(&solid(b))) { + *k += n; + } + } + // Planes, rational patches (sphere) and revolved profiles (cylinder, torus) all covered. + assert!(kinds[0] > 0 && kinds[3] > 0 && kinds[1] + kinds[2] > 0, "{kinds:?}"); + } + + #[test] + fn transformed_and_mirrored_bodies_evaluate_like_truck() { + let Ok(cyl) = crate::cylinder(Vec3::ZERO, Vec3::Z, 2.0, 5.0) else { panic!("cylinder") }; + // Column-major: x → 3 − x, z → z + 1. + let mirror = [[-1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [3.0, 0.0, 1.0, 1.0]]; + let Ok(m) = crate::transform_matrix(&cyl, mirror) else { panic!("transform") }; + assert!(assert_matches_truck(&solid(m)).iter().sum::() > 0); + } +} diff --git a/crates/kernel/src/tess/tests.rs b/crates/kernel/src/tess/tests.rs new file mode 100644 index 00000000..8b756cfb --- /dev/null +++ b/crates/kernel/src/tess/tests.rs @@ -0,0 +1,245 @@ +use std::collections::HashMap; + +use solidcraft_geom::{Mesh, Vec3}; + +use super::*; +use crate::{Body, BoolOp}; + +/// Every triangle edge, with vertices merged by exact position, is shared by exactly two +/// triangles: the mesh is closed and its faces meet exactly (no cracks, no T-junctions). +fn assert_watertight(m: &Mesh, name: &str) { + // By value: −0.0 and 0.0 are one position (adding 0.0 makes −0.0 into 0.0). + let key = |i: u32| m.positions.get(i as usize).map(|p| [(p.x + 0.0).to_bits(), (p.y + 0.0).to_bits(), (p.z + 0.0).to_bits()]); + let mut uses: HashMap<([u64; 3], [u64; 3]), i32> = HashMap::new(); + for t in &m.triangles { + for e in 0..3 { + let (Some(a), Some(b)) = (key(t[e]), key(t[(e + 1) % 3])) else { panic!("{name}: bad index") }; + if a == b { + continue; + } + // Directed: a closed, consistently oriented mesh uses each edge once each way. + *uses.entry((a, b)).or_default() += 1; + *uses.entry((b, a)).or_default() -= 1; + } + } + let open: Vec<_> = uses.iter().filter(|(_, n)| **n != 0).collect(); + assert!(open.is_empty(), "{name}: {} unmatched directed edges, e.g. {:?}", open.len(), open.first()); +} + +fn plan_body(b: &Body, tol: f64) -> (Mesh, PlanStats) { + match plan(&b.solid, tol) { + Ok((p, stats)) => (p.evaluate(), stats), + Err(e) => panic!("plan: {e}"), + } +} + +fn check(name: &str, b: &Body, tol: f64, volume: f64) -> PlanStats { + let (m, stats) = plan_body(b, tol); + assert_watertight(&m, name); + let got = m.measure().volume; + // Chord error tol over the surface: the volume is short by at most about area · tol. + let area = m.measure().area; + assert!((got - volume).abs() <= area * tol * 1.5 + 1e-9, "{name}: volume {got} vs {volume} (area {area}, tol {tol})"); + // Outward normals: the divergence theorem gives a positive volume only if they are. + assert!(got > 0.0, "{name}: inside out"); + let faces = b.solid.face_iter().count(); + assert_eq!(stats.planar + stats.structured + stats.delaunay, faces, "{name}: every face planned"); + stats +} + +fn ok(r: crate::Result) -> T { + match r { + Ok(v) => v, + Err(e) => panic!("{e}"), + } +} + +fn cut(a: &Body, b: &Body) -> Body { + match crate::boolean(a, b, BoolOp::Cut) { + Ok(Some(r)) => r, + other => panic!("cut: {other:?}"), + } +} + +const PI: f64 = std::f64::consts::PI; + +#[test] +fn primitives_are_watertight_with_the_right_volume() { + for tol in [0.05, 0.005] { + let s = check("box", &ok(crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 6.0, 4.0))), tol, 240.0); + assert_eq!(s, PlanStats { planar: 6, structured: 0, delaunay: 0 }); + let s = check("cylinder", &ok(crate::cylinder(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 1.0, 1.0), 2.5, 7.0)), tol, PI * 2.5 * 2.5 * 7.0); + assert_eq!(s.delaunay, 0, "cylinder walls are structured"); + let s = check("torus", &ok(crate::torus(Vec3::ZERO, 6.0, 1.5)), tol, 2.0 * PI * PI * 6.0 * 1.5 * 1.5); + assert_eq!(s.delaunay, 0, "torus patches are structured"); + } +} + +#[test] +fn holes_and_trimmed_faces_are_watertight() { + let tol = 0.01; + let bx = ok(crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 4.0))); + let pin = ok(crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0)); + check("box with a hole", &cut(&bx, &pin), tol, 400.0 - PI * 4.0 * 4.0); + // A cylinder with a cross hole: trimmed cylinder halves (Delaunay) and intersection-curve walls. + let cyl = ok(crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0)); + let cross = ok(crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0)); + let body = cut(&cyl, &cross); + let (m, stats) = plan_body(&body, tol); + assert_watertight(&m, "cross hole"); + assert!(stats.delaunay > 0, "{stats:?}"); + // Against the reference mesher's volume (no closed form for the two-cylinder intersection). + let Ok(reference) = body.tessellate(tol * 0.2) else { panic!("reference mesh") }; + let (a, b) = (m.measure().volume, reference.measure().volume); + assert!((a - b).abs() < 0.01 * b, "cross hole volume {a} vs {b}"); +} + +#[test] +fn faces_and_edges_keep_the_body_order() { + let body = ok(crate::cylinder(Vec3::ZERO, Vec3::Z, 2.0, 5.0)); + let (m, _) = plan_body(&body, 0.01); + let Ok(reference) = body.tessellate(0.01) else { panic!("reference mesh") }; + assert_eq!(m.edges.len(), reference.edges.len()); + assert_eq!(m.edge_faces, reference.edge_faces); + // Same faces in the same order: each face's area agrees with the reference mesher's. + let face_area = |m: &Mesh, f: u32| -> f64 { + m.triangles + .iter() + .zip(&m.tri_face) + .filter(|(_, tf)| **tf == f) + .map(|(t, _)| { + let p = |i: u32| m.positions.get(i as usize).copied().unwrap_or_default(); + (p(t[1]) - p(t[0])).cross(p(t[2]) - p(t[0])).len() * 0.5 + }) + .sum() + }; + for f in 0..body.solid.face_iter().count() as u32 { + let (a, b) = (face_area(&m, f), face_area(&reference, f)); + assert!((a - b).abs() < 0.01 * b.max(1e-9), "face {f}: {a} vs {b}"); + } + // Each edge polyline ends on the same points as the reference's. + for (a, b) in m.edges.iter().zip(&reference.edges) { + assert!(a.first().zip(b.first()).is_some_and(|(p, q)| p.dist(*q) < 1e-9)); + assert!(a.last().zip(b.last()).is_some_and(|(p, q)| p.dist(*q) < 1e-9)); + } +} + +/// Side-by-side wireframes, reference mesher (left) and plan (right): +/// `SOLIDCRAFT_TESS_SVG= cargo test -p solidcraft-kernel render_wireframes -- --ignored`. +#[test] +#[ignore] +fn render_wireframes() { + let Ok(dir) = std::env::var("SOLIDCRAFT_TESS_SVG") else { return }; + let tol = 0.02; + let bx = ok(crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 4.0))); + let pin = ok(crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0)); + let cross = ok(crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0)); + let filleted = + ok(crate::fillet(&bx, &[Vec3::new(5.0, 0.0, 4.0), Vec3::new(10.0, 5.0, 4.0), Vec3::new(5.0, 10.0, 4.0), Vec3::new(0.0, 5.0, 4.0)], 1.5)); + let shapes = [ + ("box_hole", cut(&bx, &pin)), + ("cylinder", pin.clone()), + ("cross_hole", cut(&pin, &cross)), + ("torus", ok(crate::torus(Vec3::ZERO, 6.0, 1.5))), + ("fillet", filleted), + ]; + for (name, body) in shapes { + let t0 = std::time::Instant::now(); + let Ok(old) = body.tessellate(tol) else { panic!("reference") }; + let t_old = t0.elapsed().as_secs_f64() * 1e3; + let t0 = std::time::Instant::now(); + let new = match plan(&body.solid, tol) { + Ok((p, st)) => Some((p.evaluate(), st)), + Err(e) => { + println!("{name}: no plan ({e})"); + None + } + }; + let t_new = t0.elapsed().as_secs_f64() * 1e3; + let mut svg = String::from(""); + let label = |x: f64, text: String| format!("{text}"); + svg += &label(20.0, format!("meshalgo (Delaunay): {} tris, {t_old:.1} ms", old.triangles.len())); + svg += &draw(&old, 0.0); + if let Some((m, st)) = &new { + svg += &label( + 820.0, + format!("plan: {} tris, {t_new:.1} ms [{} planar, {} grid, {} cdt]", m.triangles.len(), st.planar, st.structured, st.delaunay), + ); + svg += &draw(m, 800.0); + } + svg += ""; + let _ = std::fs::write(format!("{dir}/{name}.svg"), svg); + } +} + +/// Isometric, painter-sorted, flat-shaded triangles with their edges, in an 800×800 panel. +fn draw(m: &Mesh, x0: f64) -> String { + let (c30, s30) = (30f64.to_radians().cos(), 30f64.to_radians().sin()); + let proj = |p: Vec3| -> (f64, f64, f64) { + // Turn 45° about Z, tilt 30° about X; depth along the view. + let (x, y) = ((p.x - p.y) * std::f64::consts::FRAC_1_SQRT_2, (p.x + p.y) * std::f64::consts::FRAC_1_SQRT_2); + (x, p.z * c30 - y * s30, y * c30 + p.z * s30) + }; + let pts: Vec<(f64, f64, f64)> = m.positions.iter().map(|p| proj(*p)).collect(); + let (mut lo, mut hi) = ((f64::MAX, f64::MAX), (f64::MIN, f64::MIN)); + for p in &pts { + lo = (lo.0.min(p.0), lo.1.min(p.1)); + hi = (hi.0.max(p.0), hi.1.max(p.1)); + } + let k = 700.0 / (hi.0 - lo.0).max(hi.1 - lo.1).max(1e-9); + let sx = |p: (f64, f64, f64)| (x0 + 50.0 + (p.0 - lo.0) * k, 800.0 - (p.1 - lo.1) * k); + let mut tris: Vec<(f64, [u32; 3], u32)> = + m.triangles.iter().zip(&m.tri_face).map(|(t, f)| (t.iter().map(|&i| pts[i as usize].2).sum::(), *t, *f)).collect(); + tris.sort_by(|a, b| b.0.total_cmp(&a.0)); + let light = Vec3::new(0.3, -0.5, 0.8).normalized().unwrap_or(Vec3::new(0.0, 0.0, 1.0)); + let mut out = String::new(); + let tint = [(1.0, 1.0, 1.15), (1.15, 1.0, 1.0), (1.0, 1.15, 1.0), (1.1, 1.1, 0.9), (0.9, 1.1, 1.1), (1.1, 0.9, 1.1)]; + for (_, t, f) in tris { + let (tr, tg, tb) = tint[f as usize % tint.len()]; + let [a, b, c] = t.map(|i| m.positions[i as usize]); + let n = (b - a).cross(c - a).normalized().unwrap_or(Vec3::new(0.0, 0.0, 1.0)); + let shade = (90.0 + 150.0 * n.dot(light).max(0.0)) as u8; + let q: Vec<(f64, f64)> = t.iter().map(|&i| sx(pts[i as usize])).collect(); + // Back faces away (screen y points down, so outward triangles wind clockwise here). + if (q[1].0 - q[0].0) * (q[2].1 - q[0].1) - (q[2].0 - q[0].0) * (q[1].1 - q[0].1) > 0.0 { + continue; + } + let pp: Vec = q.iter().map(|(x, y)| format!("{x:.1},{y:.1}")).collect(); + out += &format!( + "", + pp.join(" "), + (f64::from(shade) * tr).min(255.0) as u8, + (f64::from(shade) * tg).min(255.0) as u8, + (f64::from(shade) * tb).min(255.0) as u8, + ); + } + out +} + +/// Blind drilled holes (a box cut by a revolved drill with a 118° point, as the hole feature +/// builds it): the point's apex is a pole of the cone faces, two copies of one boundary node. +/// The plan is closed and has no triangle on both copies (no zero-area triangles), also where +/// such a triangle carried a boundary side of the cone. +#[test] +fn a_drill_point_apex_leaves_no_zero_area_triangles() { + use solidcraft_geom::{Loop2, Plane, Region2, Vec2}; + // (box size, hole position, radius, depth) + let cases = + [(Vec3::new(40.0, 30.0, 15.0), Vec3::new(20.0, 15.0, 15.0), 3.4, 8.0), (Vec3::new(60.0, 40.0, 15.0), Vec3::new(50.0, 20.0, 15.0), 3.0, 10.0)]; + for (size, at, r, depth) in cases { + let bx = ok(crate::box_solid(Vec3::ZERO, size)); + let tip = r / (118f64 / 2.0).to_radians().tan(); + let Some(rp) = Plane::new(at, Vec3::X, Vec3::Z) else { panic!("plane") }; + let pts = [Vec2::new(0.0, 0.6), Vec2::new(r, 0.6), Vec2::new(r, -depth), Vec2::new(0.0, -depth - tip)]; + let region = Region2 { outer: Loop2::polygon(&pts), holes: vec![] }; + let Some(drill) = ok(crate::revolve(&rp, &[region], Vec2::ZERO, Vec2::Y, std::f64::consts::TAU)).pop() else { panic!("drill") }; + let body = cut(&bx, &drill); + let name = format!("hole r {r} depth {depth}"); + let (m, _) = plan_body(&body, 0.02); + assert_watertight(&m, &name); + for t in &m.triangles { + let [a, b, c] = t.map(|i| m.positions[i as usize]); + assert!((b - a).cross(c - a).len() > 1e-12, "{name}: zero-area triangle {a:?} {b:?} {c:?}"); + } + } +} diff --git a/crates/kernel/src/tess/uv.rs b/crates/kernel/src/tess/uv.rs new file mode 100644 index 00000000..2b2457fb --- /dev/null +++ b/crates/kernel/src/tess/uv.rs @@ -0,0 +1,246 @@ +//! Parameters of points on a surface, on the evaluator's own surface form (so the parameters +//! the plan stores are the ones the evaluators use), and the tolerance-driven divisions of a +//! face's parameter box. + +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::EvalSurface; + +/// The parameter box of a surface and its periods. +#[derive(Clone, Copy, Debug)] +pub(crate) struct Domain { + /// Bounded ranges (`None`: unbounded, as a plane's). + pub range: [Option<(f64, f64)>; 2], + pub period: [Option; 2], +} + +pub(crate) fn domain(s: &EvalSurface) -> Domain { + let tau = std::f64::consts::TAU; + match s { + EvalSurface::Plane { .. } => Domain { range: [None, None], period: [None, None] }, + EvalSurface::Nurbs { degree, knots, count, .. } => { + let r = |d: usize| Some((*knots[d].get(degree[d])?, *knots[d].get(count[d])?)); + Domain { range: [r(0), r(1)], period: [None, None] } + } + EvalSurface::Revolved { profile, swap, .. } => { + let r = (|| Some((*profile.knots.get(profile.degree)?, *profile.knots.get(profile.ctrl.len())?)))(); + let (a, b) = ((r, None), (Some((0.0, tau)), Some(tau))); + let (u, v) = if *swap { (b, a) } else { (a, b) }; + Domain { range: [u.0, v.0], period: [u.1, v.1] } + } + } +} + +impl Domain { + fn clamp(&self, uv: [f64; 2]) -> [f64; 2] { + let mut out = uv; + for (d, x) in out.iter_mut().enumerate() { + if let (Some((a, b)), None) = (self.range[d], self.period[d]) { + *x = x.clamp(a.min(b), a.max(b)); + } + } + out + } +} + +/// Damped Gauss–Newton (Levenberg–Marquardt) for the parameters nearest to `p`, from `start`. +/// Returns the parameters and their distance to `p`. +fn newton(s: &EvalSurface, dom: &Domain, p: Vec3, start: [f64; 2]) -> Option<([f64; 2], f64)> { + let mut uv = dom.clamp(start); + let (mut q, mut su, mut sv) = s.derivs(uv[0], uv[1])?; + let mut dist = q.dist(p); + let mut lambda = 1e-6; + for _ in 0..60 { + let r = p - q; + let (a, b, c) = (su.len2(), su.dot(sv), sv.len2()); + let (g0, g1) = (su.dot(r), sv.dot(r)); + let scale = a + c; + if !(scale > 0.0) || !scale.is_finite() { + break; + } + let (a2, c2) = (a + lambda * scale, c + lambda * scale); + let det = a2 * c2 - b * b; + if !(det.abs() > 0.0) { + break; + } + let step = [(c2 * g0 - b * g1) / det, (a2 * g1 - b * g0) / det]; + let next = dom.clamp([uv[0] + step[0], uv[1] + step[1]]); + let Some((nq, nsu, nsv)) = s.derivs(next[0], next[1]) else { break }; + let nd = nq.dist(p); + if nd <= dist { + let moved = (next[0] - uv[0]).abs() + (next[1] - uv[1]).abs(); + (uv, q, su, sv, dist) = (next, nq, nsu, nsv, nd); + lambda = (lambda * 0.1).max(1e-12); + if moved <= 1e-15 * (1.0 + uv[0].abs() + uv[1].abs()) { + break; + } + } else { + lambda *= 10.0; + if lambda > 1e8 { + break; + } + } + } + Some((uv, dist)) +} + +/// Parameters of `p` (a point on the surface within `accept`), from `hint` when that +/// converges, else from the best starts of a grid over the surface's domain. +pub(crate) fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> { + if let EvalSurface::Plane { origin, du, dv } = s { + let r = p - *origin; + let (a, b, c) = (du.len2(), du.dot(*dv), dv.len2()); + let det = a * c - b * b; + if !(det.abs() > 0.0) { + return None; + } + let (g0, g1) = (du.dot(r), dv.dot(r)); + let uv = [(c * g0 - b * g1) / det, (a * g1 - b * g0) / det]; + return (uv[0].is_finite() && uv[1].is_finite()).then_some(uv); + } + let mut best: Option<([f64; 2], f64)> = None; + if let Some(h) = hint + && let Some((uv, d)) = newton(s, dom, p, h) + { + if d <= accept { + return Some(uv); + } + best = Some((uv, d)); + } + let (ur, vr) = (dom.range[0]?, dom.range[1]?); + const N: usize = 16; + let mut starts: Vec<(f64, [f64; 2])> = Vec::with_capacity((N + 1) * (N + 1)); + for i in 0..=N { + for j in 0..=N { + let uv = [ur.0 + (ur.1 - ur.0) * i as f64 / N as f64, vr.0 + (vr.1 - vr.0) * j as f64 / N as f64]; + if let Some((q, _, _)) = s.derivs(uv[0], uv[1]) { + starts.push(((q - p).len2(), uv)); + } + } + } + starts.sort_by(|a, b| a.0.total_cmp(&b.0)); + for (_, start) in starts.iter().take(4) { + if let Some((uv, d)) = newton(s, dom, p, *start) + && best.is_none_or(|(_, bd)| d < bd) + { + best = Some((uv, d)); + } + } + best.filter(|(_, d)| *d <= accept).map(|(uv, _)| uv) +} + +/// `x` moved by whole periods to lie within half a period of `near`. +pub(crate) fn unwrap(x: f64, near: f64, period: Option) -> f64 { + match period { + Some(p) if p > 0.0 => x - ((x - near) / p).round() * p, + _ => x, + } +} + +/// Is the surface degenerate along direction `d` at `uv` (a pole: moving along `d` does not +/// move the point)? +pub(crate) fn degenerate(s: &EvalSurface, uv: [f64; 2], d: usize) -> bool { + let Some((_, su, sv)) = s.derivs(uv[0], uv[1]) else { return false }; + let (this, other) = if d == 0 { (su, sv) } else { (sv, su) }; + this.len() <= 1e-9 * other.len().max(1e-12) +} + +/// Division of `range` along direction `d` so that the surface between neighbouring lines +/// sags at most `tol` from its chord, checked along iso-lines across `across`; at most +/// `max_lines` lines, ends included. +pub(crate) fn divide(s: &EvalSurface, d: usize, range: (f64, f64), across: (f64, f64), tol: f64, max_lines: usize) -> Vec { + let (a0, a1) = range; + if matches!(s, EvalSurface::Plane { .. }) || !(a1 > a0) { + return vec![a0, a1]; + } + let at = |x: f64, y: f64| { + let (u, v) = if d == 0 { (x, y) } else { (y, x) }; + s.derivs(u, v).map(|(p, _, _)| p) + }; + const ACROSS: usize = 5; + let ys: Vec = (0..ACROSS).map(|k| across.0 + (across.1 - across.0) * k as f64 / (ACROSS - 1) as f64).collect(); + let flat = |a: f64, b: f64| { + ys.iter().all(|&y| { + let (Some(pa), Some(pb)) = (at(a, y), at(b, y)) else { return true }; + [0.25, 0.5, 0.75].iter().all(|f| at(a + (b - a) * f, y).is_none_or(|q| q.dist_to_segment(pa, pb) <= tol)) + }) + }; + // Start from quarter turns at most, so a closed iso-line is not taken for a flat one. + let initial = ((a1 - a0) / std::f64::consts::FRAC_PI_2).ceil().clamp(1.0, 64.0) as u32; + // Both ends exactly (a0 + (a1 − a0) · 1 may round away from a1): the end lines are the + // face's own parameter box, which the boundary nodes sit on. + let x_at = |i: u32| { + let t = f64::from(i) / f64::from(initial); + a0 * (1.0 - t) + a1 * t + }; + let mut out = vec![a0]; + let mut stack: Vec<(f64, f64, u32)> = (0..initial).rev().map(|i| (x_at(i), x_at(i + 1), 0)).collect(); + while let Some((a, b, depth)) = stack.pop() { + if depth >= 20 || out.len() + stack.len() >= max_lines.max(2) || flat(a, b) { + out.push(b); + } else { + let m = 0.5 * (a + b); + stack.push((m, b, depth + 1)); + stack.push((a, m, depth + 1)); + } + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + use solidcraft_geom::tess::EvalCurve; + + fn cylinder() -> EvalSurface { + EvalSurface::Revolved { + profile: EvalCurve { degree: 1, knots: vec![0.0, 0.0, 1.0, 1.0], ctrl: vec![[3.0, 0.0, 0.0, 1.0], [3.0, 0.0, 4.0, 1.0]] }, + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]], + swap: false, + } + } + + #[test] + fn projects_points_back_to_their_parameters() { + let s = cylinder(); + let dom = domain(&s); + for &(u, v) in &[(0.1, 0.2), (0.5, 3.0), (0.9, 6.0), (0.0, 1.0)] { + let (p, _) = s.eval(u, v); + let Some(uv) = project(&s, &dom, p, None, 1e-9) else { panic!("projects") }; + let (q, _) = s.eval(uv[0], uv[1]); + assert!(q.dist(p) < 1e-9); + assert!((uv[0] - u).abs() < 1e-9); + } + // Off the surface by more than `accept`: no parameters. + assert!(project(&s, &dom, Vec3::new(5.0, 0.0, 1.0), None, 1e-3).is_none()); + } + + #[test] + fn divides_the_turn_but_not_the_straight_direction() { + let s = cylinder(); + let turn = divide(&s, 1, (0.0, std::f64::consts::TAU), (0.0, 1.0), 0.01, 10_000); + let along = divide(&s, 0, (0.0, 1.0), (0.0, std::f64::consts::TAU), 0.01, 10_000); + assert_eq!(along, vec![0.0, 1.0]); + // Chord sag r(1 − cos(θ/2)) ≤ 0.01 with r = 3 needs θ ≤ 0.163: at least 39 segments. + assert!(turn.len() >= 40 && turn.len() < 80, "{}", turn.len()); + } + + /// The first and last lines are the range's ends to the last bit: the boundary nodes sit on + /// them, and a line a rounding inside would be taken for an interior one (a drill point's + /// cone, where a0 + (a1 − a0) is not a1). + #[test] + fn divides_with_the_ends_exactly() { + let (a0, a1) = (-1.010485317689877, 1.083909784703316); + let lines = divide(&cylinder(), 1, (a0, a1), (0.0, 1.0), 0.01, 10_000); + assert_eq!(lines.first().copied(), Some(a0)); + assert_eq!(lines.last().copied(), Some(a1)); + } + + #[test] + fn unwraps_across_the_period() { + let tau = std::f64::consts::TAU; + assert!((unwrap(0.1, tau - 0.1, Some(tau)) - (tau + 0.1)).abs() < 1e-12); + assert_eq!(unwrap(0.1, 6.0, None), 0.1); + } +}