Display tessellation plans (4/4): one GPU dispatch for the scene, intersection edges by Newton (#806)

Packs the scene's planned bodies into one GPU upload and dispatch, and evaluates intersection edges by seeded Newton refinement from a fitted B-spline through exactly solved leader points, falling back to the previous evaluation near tangency. Intersection nodes carry their faces' parameters. This changes speed, not what the app draws.

Part 4 of 4 of #303.

Co-authored-by: jyviko <[email protected]>

crates/kernel/src/tess/edges.rs

@@ -5,6 +5,8 @@ use solidcraft_geom::Vec3;
 use truck_meshalgo::prelude::*;
 use truck_modeling as mt;
 
+use solidcraft_geom::tess::EvalSurface;
+
 use super::size::Rules;
 use crate::body::from_p3;
 
@@ -13,11 +15,33 @@ pub(crate) const MAX_EDGE_NODES: usize = 20_000;
 /// Deepest bisection of one edge segment.
 const MAX_DEPTH: u32 = 24;
 
+/// An edge's nodes, from its absolute front to its absolute back.
+pub(crate) struct Nodes {
+    pub points: Vec<Vec3>,
+    /// For an intersection edge: its two surfaces and, per node, both surfaces' parameters
+    /// (u1, v1, u2, v2) where Newton placed the node (`None` where truck's evaluation did).
+    pub on: Option<([EvalSurface; 2], Vec<Option<[f64; 4]>>)>,
+}
+
+impl Nodes {
+    fn plain(points: Vec<Vec3>) -> Nodes {
+        Nodes { points, on: None }
+    }
+
+    /// The parameters of each node on `surface`, where the edge knows them: `None` unless
+    /// `surface` is one of the edge's two surfaces.
+    pub(crate) fn params_on(&self, surface: &EvalSurface) -> Option<Vec<Option<[f64; 2]>>> {
+        let (surfaces, params) = self.on.as_ref()?;
+        let side = surfaces.iter().position(|s| s == surface)?;
+        Some(params.iter().map(|q| q.map(|q| [q[2 * side], q[2 * side + 1]])).collect())
+    }
+}
+
 /// Nodes of `edge` from its absolute front to its absolute back: every chord within the
 /// tolerance of the curve (checked at the quarter points, so an S-bend in one segment is not
 /// missed), at most `h` long, and within the rules' facet angle. The end nodes are the
 /// vertices' own points, so edges meeting at a vertex meet exactly.
-pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec<Vec3> {
+pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Nodes {
     let tol = rules.tol;
     let front = from_p3(edge.absolute_front().point());
     let back = from_p3(edge.absolute_back().point());
@@ -32,29 +56,35 @@ pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec<Vec3> {
         if let Some(p) = pts.last_mut() {
             *p = back;
         }
-        return pts;
+        return Nodes::plain(pts);
     }
     let (t0, t1) = curve.range_tuple();
     if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
-        return vec![front, back];
+        return Nodes::plain(vec![front, back]);
     }
-    // Our evaluator where the curve has a form it knows (no allocation per point); truck's
-    // for intersection curves.
+    // Our evaluators: the curve's own form where it has one (no allocation per point); for an
+    // intersection curve, Newton on both surfaces from seeded parameters (see `ssi`), which
+    // also gives each node's parameters on both faces; truck's evaluation only where that does
+    // not converge.
     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))),
+    let ssi = if fast.is_none() { super::ssi::Intersection::new(&curve) } else { None };
+    type Pt = (Vec3, Option<[f64; 4]>);
+    let at = |t: f64| -> Option<Pt> {
+        match (&fast, &ssi) {
+            (Some(c), _) => c.eval(t).map(|(p, _)| (p, None)),
+            (None, Some(ix)) => ix.eval(t).map(|(p, q)| (p, Some(q))).or_else(|| Some((from_p3(curve.subs(t)), None))),
+            (None, None) => Some((from_p3(curve.subs(t)), None)),
         }
-        .filter(|p| p.is_finite())
+        .filter(|(p, _)| p.is_finite())
     };
+    let fallback = || Nodes::plain(vec![front, back]);
     // A closed curve's chord between its ends is empty: start it from a few segments. An open
     // one starts whole, so an arc the tolerance leaves flat stays one segment.
     let initial: u32 = if edge.absolute_front().id() == edge.absolute_back().id() { 4 } else { 1 };
     let t_at = |i: u32| t0 + (t1 - t0) * f64::from(i) / f64::from(initial);
-    let mut ends: Vec<Vec3> = Vec::with_capacity(initial as usize + 1);
+    let mut ends: Vec<Pt> = Vec::with_capacity(initial as usize + 1);
     for i in 0..=initial {
-        let Some(p) = at(t_at(i)) else { return vec![front, back] };
+        let Some(p) = at(t_at(i)) else { return fallback() };
         ends.push(p);
     }
     // Depth-first bisection, segment by segment, keeping the output ordered. Each segment
@@ -63,16 +93,17 @@ pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec<Vec3> {
     struct Seg {
         a: f64,
         b: f64,
-        pa: Vec3,
-        pb: Vec3,
-        mid: Option<Vec3>,
+        pa: Pt,
+        pb: Pt,
+        mid: Option<Pt>,
         depth: u32,
     }
-    let mut out: Vec<Vec3> = vec![front];
+    let Some(&first) = ends.first() else { return fallback() };
+    let mut out: Vec<Pt> = vec![(front, first.1)];
     let mut stack: Vec<Seg> = Vec::with_capacity(64);
     for i in (0..initial).rev() {
-        let (Some(pa), Some(pb)) = (<[Vec3]>::get(&ends, i as usize).copied(), <[Vec3]>::get(&ends, i as usize + 1).copied()) else {
-            return vec![front, back];
+        let (Some(pa), Some(pb)) = (<[Pt]>::get(&ends, i as usize).copied(), <[Pt]>::get(&ends, i as usize + 1).copied()) else {
+            return fallback();
         };
         stack.push(Seg { a: t_at(i), b: t_at(i + 1), pa, pb, mid: None, depth: 0 });
     }
@@ -86,11 +117,11 @@ pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec<Vec3> {
         let (q1, q3) = (at(a + 0.25 * (b - a)), at(a + 0.75 * (b - a)));
         let done = sg.depth >= MAX_DEPTH || out.len() + stack.len() >= MAX_EDGE_NODES;
         let fine = || {
-            let sag = [q1, pm, q3].iter().flatten().map(|q| q.dist_to_segment(pa, pb)).fold(0.0, f64::max);
-            let short = !(pa.dist(pb) > h * (1.0 + 1e-9));
+            let sag = [q1, pm, q3].iter().flatten().map(|q| q.0.dist_to_segment(pa.0, pb.0)).fold(0.0, f64::max);
+            let short = !(pa.0.dist(pb.0) > h * (1.0 + 1e-9));
             // The two half chords of an arc turning θ differ by θ/2.
             let straight = match (rules.turn, pm) {
-                (Some(t), Some(q)) => !t.exceeds(2.0 * turn(q - pa, pb - q), sag),
+                (Some(t), Some(q)) => !t.exceeds(2.0 * turn(q.0 - pa.0, pb.0 - q.0), sag),
                 _ => true,
             };
             short && straight && sag <= tol
@@ -104,12 +135,13 @@ pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec<Vec3> {
         }
     }
     if out.len() < 2 {
-        return vec![front, back];
+        return fallback();
     }
     if let Some(p) = out.last_mut() {
-        *p = back;
+        p.0 = back;
     }
-    out
+    let on = ssi.map(|ix| (ix.surfaces().clone(), out.iter().map(|p| p.1).collect()));
+    Nodes { points: out.into_iter().map(|p| p.0).collect(), on }
 }
 
 /// Angle (radians) between two directions; zero when either is empty.
@@ -130,7 +162,7 @@ mod tests {
         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, &Rules { tol, max_len: f64::INFINITY, turn: None }, f64::INFINITY);
+                let pts = sample_edge(&e, &Rules { tol, max_len: f64::INFINITY, turn: None }, f64::INFINITY).points;
                 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())));
@@ -152,7 +184,7 @@ mod tests {
     fn edges_are_split_to_the_side_length_and_turn() {
         let Ok(cyl) = crate::cylinder(Vec3::ZERO, Vec3::Z, 10.0, 5.0) else { panic!("cylinder") };
         for e in cyl.solid.edge_iter() {
-            let pts = sample_edge(&e, &Rules { tol: 1.0, max_len: f64::INFINITY, turn: None }, 1.2);
+            let pts = sample_edge(&e, &Rules { tol: 1.0, max_len: f64::INFINITY, turn: None }, 1.2).points;
             for w in pts.windows(2) {
                 assert!(w[0].dist(w[1]) <= 1.2 + 1e-9, "side {}", w[0].dist(w[1]));
             }
@@ -173,8 +205,36 @@ mod tests {
             .solid
             .edge_iter()
             .filter(|e| seen.insert(e.id()) && !matches!(e.curve(), mt::Curve::Line(_)))
-            .map(|e| sample_edge(&e, &thirty, f64::INFINITY).len() - 1)
+            .map(|e| sample_edge(&e, &thirty, f64::INFINITY).points.len() - 1)
             .sum();
         assert!(sides >= 24, "{sides}");
     }
+
+    /// An intersection edge's nodes carry their parameters on both faces: each face on one of
+    /// the edge's surfaces gets them, and the surface there is the node.
+    #[test]
+    fn intersection_nodes_carry_their_parameters_on_both_faces() {
+        let Ok(pin) = crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0) else { panic!("pin") };
+        let Ok(cross) = crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0) else { panic!("cross") };
+        let Ok(Some(body)) = crate::boolean(&pin, &cross, crate::BoolOp::Cut) else { panic!("cut") };
+        let rules = Rules { tol: 0.01, max_len: f64::INFINITY, turn: None };
+        let surfaces: Vec<EvalSurface> = body.solid.face_iter().filter_map(|f| super::super::surface::eval_surface(&f.surface())).collect();
+        let (mut edges, mut checked) = (0, 0);
+        let mut seen = std::collections::HashSet::new();
+        for e in body.solid.edge_iter().filter(|e| seen.insert(e.id()) && matches!(e.curve(), mt::Curve::IntersectionCurve(_))) {
+            let n = sample_edge(&e, &rules, f64::INFINITY);
+            let on: Vec<_> = surfaces.iter().filter_map(|s| n.params_on(s).map(|q| (s, q))).collect();
+            assert_eq!(on.len(), 2, "both faces of an intersection edge know its parameters");
+            edges += 1;
+            for (s, q) in on {
+                for (p, uv) in n.points.iter().zip(q).skip(1).take(n.points.len().saturating_sub(2)) {
+                    let Some(uv) = uv else { continue };
+                    let Some((at, _, _)) = s.derivs(uv[0], uv[1]) else { panic!("evaluates") };
+                    assert!(at.dist(*p) < 1e-7, "{at:?} vs {p:?}");
+                    checked += 1;
+                }
+            }
+        }
+        assert!(edges > 0 && checked > 20, "{edges} edges, {checked} nodes");
+    }
 }

crates/kernel/src/tess/mod.rs

@@ -21,6 +21,7 @@ pub(crate) mod coons;
 pub(crate) mod edges;
 pub(crate) mod grid;
 pub(crate) mod size;
+pub(crate) mod ssi;
 pub(crate) mod surface;
 pub(crate) mod uv;
 
@@ -83,13 +84,13 @@ pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>) -> Result<
     let sampled = crate::par::map("edge nodes", &jobs, solved >= PARALLEL_CURVES || edges.len() >= PARALLEL_EDGES, |(e, h)| {
         Ok(edges::sample_edge(e, &rules, *h))
     });
-    let mut nodes: HashMap<_, Vec<Vec3>> = HashMap::with_capacity(edges.len());
-    for (e, pts) in edges.iter().zip(sampled) {
-        nodes.insert(e.id(), pts?);
+    let mut nodes: HashMap<_, edges::Nodes> = HashMap::with_capacity(edges.len());
+    for (e, n) in edges.iter().zip(sampled) {
+        nodes.insert(e.id(), n?);
     }
     relax_crowded_faces(&faces, &jobs, &rules, &mut nodes);
     let oriented = |e: &mt::Edge| -> Option<Vec<Vec3>> {
-        let mut pts = nodes.get(&e.id())?.clone();
+        let mut pts = nodes.get(&e.id())?.points.clone();
         if !e.orientation() {
             pts.reverse();
         }
@@ -102,21 +103,32 @@ pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>) -> Result<
             Some((Some(eval), Some(size))) => (eval.clone(), *size),
             _ => return Err(unsupported(format!("face {fi}: surface not evaluable"))),
         };
-        // Boundary nodes, loop by loop (each edge's last node is the next edge's first).
+        // Boundary nodes, loop by loop (each edge's last node is the next edge's first), with
+        // their parameters on this face where the edge knows them (an intersection edge on
+        // this face's surface).
         let mut loops3: Vec<Vec<Vec3>> = Vec::new();
+        let mut known: Vec<Vec<Option<[f64; 2]>>> = Vec::new();
         for wire in face.absolute_boundaries() {
-            let mut pts = Vec::new();
+            let (mut pts, mut uvs) = (Vec::new(), Vec::new());
             for e in wire.edge_iter() {
-                let mut ep = oriented(e).ok_or_else(|| unsupported("edge without nodes"))?;
+                let n = nodes.get(&e.id()).ok_or_else(|| unsupported("edge without nodes"))?;
+                let (mut ep, mut eq) = (n.points.clone(), n.params_on(&eval).unwrap_or_else(|| vec![None; n.points.len()]));
+                if !e.orientation() {
+                    ep.reverse();
+                    eq.reverse();
+                }
                 ep.pop();
+                eq.pop();
                 pts.extend(ep);
+                uvs.extend(eq);
             }
             if pts.len() >= 2 {
                 loops3.push(pts);
+                known.push(uvs);
             }
         }
         let mut stats = PlanStats::default();
-        let fp = plan_face(fi, &eval, &loops3, &rules, size, &mut stats)?;
+        let fp = plan_face(fi, &eval, &loops3, &known, &rules, size, &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() },
@@ -146,7 +158,7 @@ pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>) -> Result<
 /// holes is not worth a hundred times the work to draw each hole round). An edge keeps one set
 /// of nodes for both its faces, at the coarser of their levels, so every rim of a crowded face
 /// gets the same share whatever the order of the faces.
-fn relax_crowded_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &size::Rules, nodes: &mut HashMap<mt::EdgeID, Vec<Vec3>>) {
+fn relax_crowded_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &size::Rules, nodes: &mut HashMap<mt::EdgeID, edges::Nodes>) {
     use truck_meshalgo::tessellation::Deflection;
     let Some(turn) = rules.turn else { return };
     let h_of: HashMap<mt::EdgeID, (&mt::Edge, f64)> = jobs.iter().map(|(e, h)| (e.id(), (*e, *h))).collect();
@@ -154,7 +166,7 @@ fn relax_crowded_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &s
         faces.iter().map(|f| f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().map(|e| e.id())).collect()).collect();
     let mut level: HashMap<mt::EdgeID, u8> = HashMap::new();
     for l in 1..=Deflection::MAX_RELAX {
-        let segments = |id: &mt::EdgeID| nodes.get(id).map_or(0, |n| n.len().saturating_sub(1));
+        let segments = |id: &mt::EdgeID| nodes.get(id).map_or(0, |n| n.points.len().saturating_sub(1));
         let mut crowded: Vec<mt::EdgeID> = boundaries
             .iter()
             .filter(|ids| ids.iter().map(segments).sum::<usize>() > Deflection::MAX_ANGLE_BOUNDARY)
@@ -184,6 +196,7 @@ fn plan_face(
     fi: usize,
     s: &EvalSurface,
     loops3: &[Vec<Vec3>],
+    known: &[Vec<Option<[f64; 2]>>],
     rules: &size::Rules,
     size: (f64, [f64; 2]),
     stats: &mut PlanStats,
@@ -197,11 +210,17 @@ fn plan_face(
     let mut fixed: Vec<Vec3> = Vec::new();
     let mut uvs: Vec<[f64; 2]> = Vec::new();
     let mut ranges: Vec<std::ops::Range<usize>> = Vec::new();
-    for pts in loops3 {
+    for (li, pts) in loops3.iter().enumerate() {
         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"))?;
+        for (k, p) in pts.iter().enumerate() {
+            // Known parameters (an intersection edge's), else a projection seeded by the
+            // previous node.
+            let given = known.get(li).and_then(|l| l.get(k)).copied().flatten();
+            let mut q = match given {
+                Some(q) => q,
+                None => 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]);

crates/kernel/src/tess/ssi.rs

@@ -0,0 +1,219 @@
+//! Points on a surface–surface intersection edge, on our own surface evaluators.
+//!
+//! truck stores an intersection edge as a leader curve (a polyline or B-spline near the
+//! intersection) and its two surfaces, and finds each point by projecting the leader point onto
+//! each surface in turn from no starting guess: a nearest-point search per surface per point.
+//! Here the intersection is solved exactly at a set of leader parameters, and a cubic B-spline
+//! is fitted through those points together with both surfaces' parameters there (a smooth
+//! leader with its curves on both surfaces). Each point is then refined from the B-spline by
+//! Newton's method on the four unknowns (u1, v1, u2, v2): S1(u1, v1) = S2(u2, v2), held on the
+//! plane through the B-spline point normal to the B-spline (Barnhill & Kersey, "A marching
+//! method for parametric surface/surface intersection", CAGD 7, 1990; Patrikalakis & Maekawa,
+//! "Shape Interrogation for CAD/CAM", 2002, ch. 5). Every point is independent once seeded. Where the surfaces are close to
+//! tangent along the edge (the system becomes singular) or Newton does not converge, the caller
+//! uses truck's evaluation.
+
+use solidcraft_geom::Vec3;
+use solidcraft_geom::bspline::CurveN;
+use solidcraft_geom::tess::EvalSurface;
+use truck_meshalgo::prelude::*;
+use truck_modeling as mt;
+
+use super::uv::{self, Domain};
+use crate::body::from_p3;
+
+/// Leader parameters at which both surfaces' parameters are found up front.
+const SEEDS: usize = 32;
+/// Newton iterations per point (two or three suffice from a seed).
+const MAX_ITER: usize = 8;
+/// Below this sine of the angle between the surface normals the intersection is treated as
+/// near-tangent: points there are left to truck.
+const MIN_SIN: f64 = 0.05;
+
+/// An intersection edge ready for point evaluation.
+pub(crate) struct Intersection {
+    surfaces: [EvalSurface; 2],
+    domains: [Domain; 2],
+    /// The smooth leader over the curve's parameter: the intersection point (x, y, z) and both
+    /// surfaces' parameters (u1, v1, u2, v2), periodic directions unwrapped along it.
+    leader: CurveN<7>,
+    /// Convergence tolerance (mm).
+    eps: f64,
+}
+
+impl Intersection {
+    /// `None` unless `curve` is an intersection curve whose surfaces have evaluators and whose
+    /// leader projects onto both.
+    pub(crate) fn new(curve: &mt::Curve) -> Option<Self> {
+        let mt::Curve::IntersectionCurve(ic) = curve else { return None };
+        let surfaces = [super::surface::eval_surface(ic.surface0())?, super::surface::eval_surface(ic.surface1())?];
+        let domains = [uv::domain(&surfaces[0]), uv::domain(&surfaces[1])];
+        let (t0, t1) = curve.range_tuple();
+        if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
+            return None;
+        }
+        let leader = |t: f64| match ic.leader() {
+            mt::Leader::Polyline(l) => (from_p3(l.subs(t)), from_v3(l.der(t))),
+            mt::Leader::BSpline(l) => (from_p3(l.subs(t)), from_v3(l.der(t))),
+        };
+        // The leader lies within the curve's tolerance of both surfaces.
+        let accept = (ic.tolerance() * 100.0).max(1e-6);
+        let scale = (0..=SEEDS).map(|k| leader(t0 + (t1 - t0) * k as f64 / SEEDS as f64).0.len()).fold(1.0, f64::max);
+        let mut it = Intersection { surfaces, domains, leader: CurveN { degree: 1, knots: Vec::new(), ctrl: Vec::new() }, eps: scale * 1e-10 };
+        let (mut params, mut points) = (Vec::with_capacity(SEEDS + 1), Vec::with_capacity(SEEDS + 1));
+        let mut prev: Option<[f64; 4]> = None;
+        for k in 0..=SEEDS {
+            let t = t0 + (t1 - t0) * k as f64 / SEEDS as f64;
+            let (p, d) = leader(t);
+            let mut q = [0.0; 4];
+            for (i, (s, dom)) in it.surfaces.iter().zip(&it.domains).enumerate() {
+                let hint = prev.map(|w| [w[2 * i], w[2 * i + 1]]);
+                let mut uv = uv::project(s, dom, p, hint, accept)?;
+                if let Some(h) = hint {
+                    for d in 0..2 {
+                        uv[d] = uv::unwrap(uv[d], h[d], dom.period[d]);
+                    }
+                }
+                q[2 * i] = uv[0];
+                q[2 * i + 1] = uv[1];
+            }
+            // Exactly on the intersection where Newton converges (near-tangent seeds keep the
+            // leader point and its projections).
+            let (p, q) = d.normalized().and_then(|dir| it.newton(p, dir, q)).unwrap_or((p, q));
+            params.push(t);
+            points.push([p.x, p.y, p.z, q[0], q[1], q[2], q[3]]);
+            prev = Some(q);
+        }
+        it.leader = CurveN::interpolate(&params, &points, 3).ok()?;
+        Some(it)
+    }
+
+    /// The point at curve parameter `t` and both surfaces' parameters there, or `None`
+    /// (near-tangent, or Newton failed).
+    pub(crate) fn eval(&self, t: f64) -> Option<(Vec3, [f64; 4])> {
+        let (x, dx) = self.leader.eval(t)?;
+        let dir = Vec3::new(dx[0], dx[1], dx[2]).normalized()?;
+        self.newton(Vec3::new(x[0], x[1], x[2]), dir, [x[3], x[4], x[5], x[6]])
+    }
+
+    /// The two surfaces, in the order of the parameters [`Intersection::eval`] gives.
+    pub(crate) fn surfaces(&self) -> &[EvalSurface; 2] {
+        &self.surfaces
+    }
+
+    /// Newton from the parameters `x` on S1 = S2, held on the plane through `l` normal to `dir`:
+    /// the point and its parameters, or `None` (near-tangent, singular, or not converged).
+    fn newton(&self, l: Vec3, dir: Vec3, mut x: [f64; 4]) -> Option<(Vec3, [f64; 4])> {
+        for _ in 0..MAX_ITER {
+            let (Some((p1, u1, v1)), Some((p2, u2, v2))) = (self.surfaces[0].derivs(x[0], x[1]), self.surfaces[1].derivs(x[2], x[3])) else {
+                return None;
+            };
+            let gap = p1 - p2;
+            let plane = dir.dot(p1 - l);
+            let n = u1.cross(v1).normalized()?.cross(u2.cross(v2).normalized()?);
+            if n.len() < MIN_SIN {
+                return None;
+            }
+            if gap.len() <= self.eps && plane.abs() <= self.eps {
+                return Some(((p1 + p2) * 0.5, x));
+            }
+            // J · dx = −F, F = (S1 − S2, dir · (S1 − L)).
+            let j = [[u1.x, v1.x, -u2.x, -v2.x], [u1.y, v1.y, -u2.y, -v2.y], [u1.z, v1.z, -u2.z, -v2.z], [dir.dot(u1), dir.dot(v1), 0.0, 0.0]];
+            let dx = solve4(j, [-gap.x, -gap.y, -gap.z, -plane])?;
+            for (xi, d) in x.iter_mut().zip(dx) {
+                *xi += d;
+            }
+            for i in 0..2 {
+                let [u, v] = self.domains[i].clamp([x[2 * i], x[2 * i + 1]]);
+                x[2 * i] = u;
+                x[2 * i + 1] = v;
+            }
+        }
+        None
+    }
+}
+
+fn from_v3(v: mt::Vector3) -> Vec3 {
+    Vec3::new(v.x, v.y, v.z)
+}
+
+/// `a · x = b` for a 4 × 4 system by Gaussian elimination with partial pivoting; `None` when
+/// singular.
+fn solve4(mut a: [[f64; 4]; 4], mut b: [f64; 4]) -> Option<[f64; 4]> {
+    let scale = a.iter().flatten().fold(0.0f64, |m, x| m.max(x.abs()));
+    if !(scale > 0.0) {
+        return None;
+    }
+    for c in 0..4 {
+        let p = (c..4).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs()))?;
+        if !(a[p][c].abs() > 1e-12 * scale) {
+            return None;
+        }
+        a.swap(c, p);
+        b.swap(c, p);
+        let pivot = a[c];
+        for r in c + 1..4 {
+            let f = a[r][c] / pivot[c];
+            for (x, p) in a[r].iter_mut().zip(pivot).skip(c) {
+                *x -= f * p;
+            }
+            b[r] -= f * b[c];
+        }
+    }
+    let mut x = [0.0; 4];
+    for c in (0..4).rev() {
+        let s: f64 = (c + 1..4).map(|k| a[c][k] * x[k]).sum();
+        x[c] = (b[c] - s) / a[c][c];
+    }
+    x.iter().all(|v| v.is_finite()).then_some(x)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn solves_a_four_by_four_system() {
+        let a = [[2.0, 1.0, 0.0, 0.0], [1.0, 3.0, 1.0, 0.0], [0.0, 1.0, 4.0, 1.0], [0.0, 0.0, 1.0, 5.0]];
+        let x = [1.0, -2.0, 0.5, 3.0];
+        let b = [0, 1, 2, 3].map(|r| (0..4).map(|c| a[r][c] * x[c]).sum::<f64>());
+        let Some(got) = solve4(a, b) else { panic!("solves") };
+        for (g, w) in got.iter().zip(x) {
+            assert!((g - w).abs() < 1e-12);
+        }
+        assert!(solve4([[1.0; 4]; 4], [1.0; 4]).is_none(), "singular");
+    }
+
+    /// Points on the intersection edges of a boolean (a cross hole: two cylinders) lie on both
+    /// surfaces, on truck's curve (the parameterisation follows our smooth leader, so the same
+    /// parameter may sit at a slightly different place along it), and advance along the curve.
+    #[test]
+    fn intersection_points_lie_on_both_surfaces_and_on_truck_s_curve() {
+        let Ok(pin) = crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0) else { panic!("pin") };
+        let Ok(cross) = crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0) else { panic!("cross") };
+        let Ok(Some(body)) = crate::boolean(&pin, &cross, crate::BoolOp::Cut) else { panic!("cut") };
+        let mut checked = 0;
+        for e in body.solid.edge_iter() {
+            let curve = e.curve();
+            let Some(ix) = Intersection::new(&curve) else { continue };
+            let (t0, t1) = curve.range_tuple();
+            let dense: Vec<Vec3> = (0..=1000).map(|k| from_p3(curve.subs(t0 + (t1 - t0) * f64::from(k) / 1000.0))).collect();
+            let near = |p: Vec3| dense.windows(2).map(|w| p.dist_to_segment(w[0], w[1])).fold(f64::MAX, f64::min);
+            let mut last: Option<Vec3> = None;
+            for k in 0..=50 {
+                let t = t0 + (t1 - t0) * f64::from(k) / 50.0;
+                let Some((p, _)) = ix.eval(t) else { continue };
+                for (s, dom) in ix.surfaces.iter().zip(&ix.domains) {
+                    assert!(uv::project(s, dom, p, None, 1e-8).is_some(), "t {t}: {p:?} off a surface");
+                }
+                assert!(near(p) < 1e-4, "t {t}: {p:?} is {} from truck's curve", near(p));
+                if let Some(q) = last {
+                    assert!(p.dist(q) > 0.0, "t {t}: the points advance");
+                }
+                last = Some(p);
+                checked += 1;
+            }
+        }
+        assert!(checked > 100, "{checked} points checked");
+    }
+}

crates/kernel/src/tess/uv.rs

@@ -33,7 +33,7 @@ pub(crate) fn domain(s: &EvalSurface) -> Domain {
 }
 
 impl Domain {
-    fn clamp(&self, uv: [f64; 2]) -> [f64; 2] {
+    pub(crate) 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]) {