brep: a pcurve keeps its tangent where its curve is smooth through a full knot

SMLib's intersection curves carry knots of full multiplicity where the curve is in fact
differentiable (its derivative the same either side). The pcurve fit took every such knot
for a corner and fitted the two sides apart; on a thread's flank, whose u runs hundreds of
times as fast as its v, the pieces met at up to 1.4 deg in (u, v), and SMLib's check refused
the record read back (Test #7 Bad CurveKink, offset 0.9997 against cos 1 deg).

A knot that leaves a spline only continuous by its multiplicity now gets a cubic with a
tangent there when the spline's one-sided derivatives (EvalCurve::sides) agree within 1e-9
of their size. Knots where the curve does turn are corners as before.

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

crates/geom/src/brep/pcurve.rs

@@ -323,7 +323,7 @@ impl On {
         let reach = 0.25 * tol;
         // Interior knots: where the curve's own parameterisation is less smooth (a pcurve is no
         // smoother than its curve), else even.
-        let mut inner = breaks(record, t0, t1);
+        let mut inner = breaks(record, Some(curve), t0, t1);
         if inner.len() + 1 > MAX_SPANS / 4 {
             return Err(Error::TooLarge(format!("a pcurve of a curve with {} knots", inner.len())));
         }
@@ -499,7 +499,12 @@ impl On {
 /// The interior parameters where `curve` is less than twice continuously differentiable, with
 /// the knot multiplicity a cubic needs there: a spline's knots by their multiplicity, and the
 /// vertices of a polyline or an intersection curve's leader polyline (its parameter turns there).
-fn breaks(curve: &Curve, t0: f64, t1: f64) -> Vec<(f64, usize)> {
+/// A knot that leaves a spline only continuous by its multiplicity, where the spline's
+/// derivative is the same either side (`eval`, the curve's own evaluator), is once
+/// differentiable and gets a pcurve with a tangent there: two pieces fitted apart would meet at
+/// an angle, sharp in (u, v) where the surface's parameters run at very different speeds (a
+/// thread's flank, hundreds of mm per unit along the helix and one across it).
+fn breaks(curve: &Curve, eval: Option<&CurveEval>, t0: f64, t1: f64) -> Vec<(f64, usize)> {
     let spline = |knots: &[f64], count: usize| -> Vec<(f64, usize)> {
         let degree = knots.len().saturating_sub(count + 1);
         let mut out: Vec<(f64, usize)> = Vec::new();
@@ -509,7 +514,15 @@ fn breaks(curve: &Curve, t0: f64, t1: f64) -> Vec<(f64, usize)> {
                 _ => out.push((k, 1)),
             }
         }
-        out.into_iter().map(|(k, m)| (k, (3 - degree.saturating_sub(m).min(2)).clamp(1, 3))).collect()
+        out.into_iter()
+            .map(|(k, m)| {
+                let mut smooth = degree.saturating_sub(m);
+                if m == degree && eval.is_some_and(|e| differentiable_at(e, k)) {
+                    smooth = 1;
+                }
+                (k, (3 - smooth.min(2)).clamp(1, 3))
+            })
+            .collect()
     };
     let corners = |n: usize| (1..n.saturating_sub(1)).map(|i| (i as f64, 3)).filter(|k| k.0 > t0 && k.0 < t1).collect();
     match curve {
@@ -519,11 +532,23 @@ fn breaks(curve: &Curve, t0: f64, t1: f64) -> Vec<(f64, usize)> {
         Curve::RationalSpline { knots, homogeneous } => spline(knots, homogeneous.len()),
         Curve::SurfaceIntersection { leader, .. } => match &**leader {
             Curve::Polyline { points } => corners(points.len()),
-            other => breaks(other, t0, t1),
+            other => breaks(other, None, t0, t1),
         },
     }
 }
 
+/// Relative difference of a spline's derivatives either side of a knot within which it is
+/// differentiable there (rounding of derivatives made equal).
+const SAME_DERIVATIVE: f64 = 1e-9;
+
+/// Whether the spline `eval` has one derivative at its knot `k` (the same from either side).
+fn differentiable_at(eval: &CurveEval, k: f64) -> bool {
+    let CurveEval::Spline(c) = eval else { return false };
+    let Some((l, r)) = c.sides(k) else { return false };
+    let size = l.len().max(r.len());
+    size > 0.0 && (l - r).len() <= SAME_DERIVATIVE * size
+}
+
 /// Signed area in (u, v) of the loop with its uses moved by `shift` along `d`.
 fn area(uses: &[(Curve2, bool)], shift: &[f64], d: usize) -> f64 {
     let mut ring: Vec<[f64; 2]> = Vec::new();

crates/geom/src/brep/pcurve_tests.rs

@@ -528,6 +528,36 @@ fn edges_on_a_spline_winding_several_turns_stay_on_their_turn() {
     }
 }
 
+/// A curve that is smooth through a knot of full multiplicity (one smooth curve with the knot
+/// inserted three times, as SMLib's intersection curves come) gets a pcurve with one tangent
+/// there. Two pieces fitted apart meet at an angle in (u, v) (a tenth of a degree here, over the
+/// degree SMLib's check allows on a real thread), sharp where u runs hundreds of times as fast as
+/// v.
+#[test]
+fn a_smooth_curve_through_a_full_knot_gets_a_pcurve_with_a_tangent_there() {
+    let (surface, _) = thread_flank();
+    let s = SurfaceEval::new(&surface).unwrap();
+    // Along the flank as it winds, a third of a turn (u = t), weaving a little across it.
+    let [a, b] = [0.40625, 0.4234375];
+    let ts: Vec<f64> = (0..=12).map(|i| a + (b - a) * i as f64 / 12.).collect();
+    let v = |t: f64| 0.5 + 0.02 * (3. * TAU * (t - a) / (b - a)).sin();
+    let on: Vec<[f64; 3]> = ts.iter().map(|&t| s.point([t, v(t)]).unwrap()).map(|p| [p.x, p.y, p.z]).collect();
+    let c = crate::bspline::CurveN::interpolate(&ts, &on, 3).unwrap();
+    let k = 0.4167;
+    let (mut knots, mut ctrl) = (c.knots, c.ctrl.iter().map(|p| crate::bspline::HPoint::weighted(Vec3::from(*p), 1.)).collect::<Vec<_>>());
+    for _ in 0..3 {
+        (knots, ctrl) = crate::bspline::insert_knot(3, &knots, &ctrl, k).unwrap();
+    }
+    let weave = Curve::Spline { knots, points: ctrl.iter().map(|h| h.project().unwrap()).map(|p| [p.x, p.y, p.z]).collect() };
+    let pc = &pcurves(&surface, &one(&weave, [a, b]), 1e-3).unwrap()[0][0];
+    assert!(lift_error(&surface, &weave, [a, b], pc) <= 1e-3);
+    let h = 1e-7;
+    let (l, m, r) = (pc.eval(k - h).unwrap(), pc.eval(k).unwrap(), pc.eval(k + h).unwrap());
+    let (dl, dr) = ([m[0] - l[0], m[1] - l[1]], [r[0] - m[0], r[1] - m[1]]);
+    let cos = (dl[0] * dr[0] + dl[1] * dr[1]) / (dl[0].hypot(dl[1]) * dr[0].hypot(dr[1]));
+    assert!(cos > 1. - 1e-6, "the pcurve turns by {} degrees at the knot", cos.min(1.).acos().to_degrees());
+}
+
 /// Hostile records give errors, never a panic.
 #[test]
 fn hostile_records_give_errors() {

crates/geom/src/tess.rs

@@ -193,7 +193,31 @@ impl EvalCurve {
         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);
+        self.eval_on(span_of(&self.knots, p, n, t), t)
+    }
+
+    /// The first derivative at `t` from either side: on the knot spans that end and start
+    /// there where `t` is a knot (the polynomials either side of it), else the derivative twice.
+    pub fn sides(&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 right = span_of(&self.knots, p, n, t);
+        // The last non-empty span before `t`'s, when `t` is where it starts.
+        let left = match self.knots.get(right) {
+            Some(&k) if k == t => {
+                (p..right).rev().find(|&i| matches!((self.knots.get(i), self.knots.get(i + 1)), (Some(a), Some(b)) if a < b)).unwrap_or(right)
+            }
+            _ => right,
+        };
+        Some((self.eval_on(left, t)?.1, self.eval_on(right, t)?.1))
+    }
+
+    /// Point and first derivative at `t` of the polynomial on knot span `span`.
+    fn eval_on(&self, span: usize, t: f64) -> Option<(Vec3, Vec3)> {
+        let p = self.degree;
         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() {