kernel: exact pcurves of lines and splines on planes

A plane's parameters are affine in space, so the pcurve of a line or a
(rational) B-spline on a plane is the curve's own knots and weights with
its control points mapped into (u, v): no fit. Curved surfaces and
intersection curves still return "not supported yet". A test helper checks
every pcurve the way the schema validator does, its interval against the
curve's range, the lift error and the wire's closure in (u, v); the box's
line pcurves lift onto their lines' ends bit for bit.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QWhMwdqJYQKQq8RCLSL1LT

crates/kernel/src/native.rs

@@ -352,4 +352,105 @@ mod tests {
         let restored = decode_native_body(&encode_native_body(&body).unwrap()).unwrap();
         equal(&body, &restored);
     }
+
+    /// Every pcurve of every face of `body` validates as the schema validator does, spans its
+    /// curve's range, lifts onto the surface within `tol` (64 samples and every span midpoint),
+    /// and each wire closes in (u, v): consecutive uses meet, or meet across a pole (the gap
+    /// lifts to one point). Returns the wires whose closure is a whole number of periods (a loop
+    /// round a pole or the surface), as (face, wire).
+    fn check_pcurves(body: &Body, tol: f64) -> Vec<(usize, usize)> {
+        use truck_modeling::{BoundedCurve, MetricSpace, ParametricCurve, ParametricSurface};
+        let mut winding = Vec::new();
+        for (f, face) in body.solid.face_iter().enumerate() {
+            let surface = face.surface();
+            let all = pcurve::pcurves(face, tol).unwrap_or_else(|e| panic!("face {f}: {e}"));
+            let wires = face.boundaries();
+            assert_eq!(all.len(), wires.len());
+            let lift = |uv: [f64; 2]| surface.subs(uv[0], uv[1]);
+            for (w, (pcs, wire)) in all.iter().zip(&wires).enumerate() {
+                assert_eq!(pcs.len(), wire.len());
+                let mut ends = Vec::new();
+                for (i, (pc, edge)) in pcs.iter().zip(wire.edge_iter()).enumerate() {
+                    let at = format!("face {f} wire {w} edge {i}");
+                    let (knots, count, values): (&Vec<f64>, usize, Vec<f64>) = match pc {
+                        solidcraft_geom::brep::Curve2::Spline { knots, points } => (knots, points.len(), points.iter().flatten().copied().collect()),
+                        solidcraft_geom::brep::Curve2::RationalSpline { knots, homogeneous } => {
+                            assert!(homogeneous.iter().all(|h| h[2] >= 1e-12), "{at}: weight");
+                            (knots, homogeneous.len(), homogeneous.iter().flatten().copied().collect())
+                        }
+                    };
+                    let degree = knots.len() - count - 1;
+                    assert!((2..=4096).contains(&count) && (1..=16).contains(&degree) && degree < count, "{at}: size");
+                    assert!(knots.windows(2).all(|k| k[0] <= k[1]) && knots[degree] < knots[count], "{at}: knots");
+                    assert!(values.iter().chain(knots).all(|x| x.is_finite() && x.abs() <= 1e12), "{at}: numbers");
+                    let curve = edge.curve();
+                    let (t0, t1) = curve.range_tuple();
+                    assert_eq!(pc.interval(), Some([t0, t1]), "{at}: interval");
+                    let mut ts: Vec<f64> = (0..=64).map(|k| t0 + (t1 - t0) * k as f64 / 64.).collect();
+                    ts.extend(knots.windows(2).filter(|k| k[0] < k[1]).map(|k| 0.5 * (k[0] + k[1])));
+                    for t in ts {
+                        let uv = pcurve::eval(pc, t).unwrap_or_else(|| panic!("{at}: eval at {t}"));
+                        let d = lift(uv).distance(curve.subs(t));
+                        assert!(d <= tol, "{at}: lift error {d} at {t}");
+                    }
+                    let (a, b) = (pcurve::eval(pc, t0).unwrap(), pcurve::eval(pc, t1).unwrap());
+                    ends.push(if edge.orientation() { (a, b) } else { (b, a) });
+                }
+                let mut wraps = false;
+                for k in 0..ends.len() {
+                    let (end, start) = (ends[k].1, ends[(k + 1) % ends.len()].0);
+                    let gap = (end[0] - start[0]).abs().max((end[1] - start[1]).abs());
+                    if gap <= 1e-9 {
+                        continue;
+                    }
+                    // A pole: the gap's ends and middle are one point in space.
+                    let mid = [0.5 * (end[0] + start[0]), 0.5 * (end[1] + start[1])];
+                    let one = lift(end).distance(lift(start)).max(lift(end).distance(lift(mid)));
+                    if one <= tol {
+                        continue;
+                    }
+                    let periods = [surface.u_period(), surface.v_period()];
+                    let whole = (0..2).all(|d| {
+                        let g = start[d] - end[d];
+                        g.abs() <= 1e-9 || periods[d].is_some_and(|p| (g / p - (g / p).round()).abs() <= 1e-9)
+                    });
+                    assert!(
+                        whole && lift(end).distance(lift(start)) <= tol,
+                        "face {f} wire {w}: open between uses {k} and {}: {end:?} {start:?}",
+                        (k + 1) % ends.len()
+                    );
+                    wraps = true;
+                }
+                if wraps {
+                    winding.push((f, w));
+                }
+            }
+        }
+        winding
+    }
+
+    #[test]
+    fn pcurves_on_planes_are_the_curves_mapped_exactly() {
+        use truck_modeling::{BoundedCurve, ParametricCurve, ParametricSurface};
+        let cube = crate::box_solid(Vec3::new(-3., 1., 2.), Vec3::new(40., 30., 10.)).unwrap();
+        assert!(check_pcurves(&cube, 1e-9).is_empty());
+        // Bit for bit: each line's pcurve ends lift to the line's own ends.
+        for face in cube.solid.face_iter() {
+            let surface = face.surface();
+            for (pcs, wire) in pcurve::pcurves(face, 1e-9).unwrap().iter().zip(face.boundaries()) {
+                for (pc, edge) in pcs.iter().zip(wire.edge_iter()) {
+                    let solidcraft_geom::brep::Curve2::Spline { points, .. } = pc else { panic!("a line is a degree-1 spline") };
+                    let curve = edge.curve();
+                    let (t0, t1) = curve.range_tuple();
+                    assert_eq!(surface.subs(points[0][0], points[0][1]), curve.subs(t0));
+                    assert_eq!(surface.subs(points[1][0], points[1][1]), curve.subs(t1));
+                }
+            }
+        }
+        // Arcs (rational) on the cylinder's caps.
+        let cyl = crate::cylinder(Vec3::ZERO, Vec3::Z, 5., 40.).unwrap();
+        for face in cyl.solid.face_iter().filter(|f| matches!(f.surface(), truck_modeling::Surface::Plane(_))) {
+            assert!(pcurve::pcurves(face, 1e-9).unwrap().iter().flatten().all(|c| matches!(c, solidcraft_geom::brep::Curve2::RationalSpline { .. })));
+        }
+    }
 }

crates/kernel/src/native/pcurve.rs

@@ -3,13 +3,117 @@
 //! that curve's direction. Exact where the edge lies on its surface in closed form, otherwise a
 //! 2D B-spline fitted within the tolerance.
 use super::Result;
+use crate::KernelError;
+use crate::body::from_p3;
+use crate::tess::uv;
+use mt::BoundedCurve;
 use solidcraft_geom::brep as b;
+use solidcraft_geom::tess::EvalSurface;
 use truck_modeling as mt;
 
 /// For each of `face`'s boundary wires, in order, the pcurve of each of its edges, in order, in
 /// the (u, v) of the face's surface as truck evaluates it (`ParametricSurface::subs`), within
-/// `tolerance` mm when lifted onto the surface.
+/// `tolerance` mm when lifted onto the surface. Wires and edges run as `face.boundaries()` lists
+/// them (the order the adapter writes); each pcurve runs in its curve's own direction.
 pub(super) fn pcurves(face: &mt::Face, tolerance: f64) -> Result<Vec<Vec<b::Curve2>>> {
-    let _ = (face, tolerance);
-    Err(crate::KernelError::Failed("not supported yet: pcurves of truck faces".into()))
+    if !(tolerance.is_finite() && tolerance > 0.0) {
+        return Err(KernelError::Invalid(format!("pcurve tolerance {tolerance}")));
+    }
+    crate::guard("pcurves", || {
+        let surface = face.surface();
+        let mut out = Vec::new();
+        for (w, wire) in face.boundaries().iter().enumerate() {
+            let mut uses = Vec::new();
+            for (i, edge) in wire.edge_iter().enumerate() {
+                let curve = edge.curve();
+                let pc = pcurve(&surface, &curve).map_err(|e| KernelError::Failed(format!("pcurve of wire {w} edge {i}: {e}")))?;
+                uses.push(pc);
+            }
+            out.push(uses);
+        }
+        Ok(out)
+    })
+}
+
+fn pcurve(surface: &mt::Surface, curve: &mt::Curve) -> Result<b::Curve2> {
+    match surface {
+        mt::Surface::Plane(_) => on_plane(surface, curve),
+        _ => Err(KernelError::Failed("not supported yet: pcurves on curved surfaces".into())),
+    }
+}
+
+/// A line or (rational) B-spline on a plane: the plane's parameters are affine in space, so the
+/// pcurve is the curve's own knots and weights with its control points mapped.
+fn on_plane(surface: &mt::Surface, curve: &mt::Curve) -> Result<b::Curve2> {
+    let Some(plane @ EvalSurface::Plane { .. }) = crate::tess::surface::eval_surface(surface) else {
+        return Err(KernelError::Failed("plane parameters".into()));
+    };
+    let dom = uv::domain(&plane);
+    let map = |p: mt::Point3| uv::project(&plane, &dom, from_p3(p), None, 0.0).ok_or_else(|| KernelError::Failed("degenerate plane".into()));
+    Ok(match curve {
+        mt::Curve::Line(l) => {
+            let (t0, t1) = curve.range_tuple();
+            b::Curve2::Spline { knots: vec![t0, t0, t1, t1], points: vec![map(l.0)?, map(l.1)?] }
+        }
+        mt::Curve::BSplineCurve(c) => b::Curve2::Spline {
+            knots: c.knot_vec().iter().copied().collect(),
+            points: c.control_points().iter().map(|p| map(*p)).collect::<Result<_>>()?,
+        },
+        mt::Curve::NurbsCurve(c) => b::Curve2::RationalSpline {
+            knots: c.knot_vec().iter().copied().collect(),
+            homogeneous: c
+                .control_points()
+                .iter()
+                .map(|h| {
+                    let [u, v] = map(mt::Point3::new(h.x / h.w, h.y / h.w, h.z / h.w))?;
+                    Ok([u * h.w, v * h.w, h.w])
+                })
+                .collect::<Result<_>>()?,
+        },
+        mt::Curve::IntersectionCurve(_) => return Err(KernelError::Failed("not supported yet: intersection pcurves".into())),
+    })
+}
+
+/// The point of `c` at `t` (`None` outside its interval or for malformed records).
+pub(super) fn eval(c: &b::Curve2, t: f64) -> Option<[f64; 2]> {
+    let (knots, h): (&[f64], Vec<[f64; 3]>) = match c {
+        b::Curve2::Spline { knots, points } => (knots, points.iter().map(|&[u, v]| [u, v, 1.0]).collect()),
+        b::Curve2::RationalSpline { knots, homogeneous } => (knots, homogeneous.clone()),
+    };
+    let p = knots.len().checked_sub(h.len() + 1)?;
+    let span = solidcraft_geom::bspline::find_span(knots, p, h.len(), t)?;
+    let n = basis(knots, p, span, t)?;
+    let mut s = [0.0; 3];
+    for (j, nj) in n.iter().enumerate() {
+        let q = h.get(span + j - p)?;
+        for (a, x) in s.iter_mut().zip(q) {
+            *a += nj * x;
+        }
+    }
+    (s[2].abs() > 0.0).then(|| [s[0] / s[2], s[1] / s[2]])
+}
+
+/// The `p + 1` basis functions non-zero on knot span `span` at `t` (NURBS Book A2.2).
+fn basis(knots: &[f64], p: usize, span: usize, t: f64) -> Option<Vec<f64>> {
+    if span < p || span + p >= knots.len() {
+        return None;
+    }
+    let k = |i: usize| knots.get(i).copied().unwrap_or(0.0);
+    let mut n = vec![0.0; p + 1];
+    let (mut left, mut right) = (vec![0.0; p + 1], vec![0.0; p + 1]);
+    *n.first_mut()? = 1.0;
+    for j in 1..=p {
+        *left.get_mut(j)? = t - k(span + 1 - j);
+        *right.get_mut(j)? = k(span + j) - t;
+        let mut saved = 0.0;
+        for r in 0..j {
+            let (rr, lj) = (*right.get(r + 1)?, *left.get(j - r)?);
+            let den = rr + lj;
+            let tmp = if den != 0.0 { n.get(r)? / den } else { 0.0 };
+            *n.get_mut(r)? = saved + rr * tmp;
+            saved = lj * tmp;
+        }
+        *n.get_mut(j)? = saved;
+    }
+    Some(n)
 }