sew: an edge folded across a bend gets a Hermite pcurve on the bend's form

A flat-pattern edge folded across a bend is a helix, which no NURBS is: sheet metal already
gives its 3D curve as degree-7 Hermite pieces of the fold, within 1e-7 mm. Its pcurve on the
bend is now built the same way, on the edge's own parameter: degree-7 Hermite pieces matching
the point and first three derivatives of its angle about the axis (through the form's quarter
parameter, by the chain rule) and of its height along it, broken where the angle crosses a
quarter (the form's parameter is only once differentiable there), halved until they lift onto
the edge within 1e-7 of its size. An approximation by construction from exact derivatives, like
the edge; nothing is fitted. B-spline and NURBS edges alike (a NURBS's derivatives from its
homogeneous form). Breaks within rounding of each other are one.

With the pcurve fit turned off, all 44 sew cases that needed it pass; with it on, all 44 pass.

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

crates/kernel/src/sew.rs

@@ -359,6 +359,216 @@ fn chain(e: &EdgeSpec, a: Vec3, b: Vec3) -> Result<(Vec<f64>, Vec<[f64; 3]>)> {
     Ok((knots, pts))
 }
 
+/// Most Hermite pieces of a folded pcurve.
+const FOLD_PIECES: usize = 1024;
+/// How close, against the edge's size, a folded pcurve lifts onto its edge: as close as the
+/// folded edge's own chain is to the fold (sheet metal's `BEND_CHAIN_TOL`).
+const FOLD_TOL: f64 = 1e-7;
+
+/// The point and first three derivatives at `t` of a B-spline or NURBS curve (homogeneous
+/// `[x·w, y·w, z·w, w]` points; weight 1 for a B-spline): the homogeneous derivatives `A⁽ᵏ⁾`,
+/// `w⁽ᵏ⁾`, then `C⁽ᵏ⁾ = (A⁽ᵏ⁾ − Σⱼ₌₁ᵏ C(k, j) w⁽ʲ⁾ C⁽ᵏ⁻ʲ⁾) / w`.
+fn curve_derivs(knots: &[f64], points: &[[f64; 4]], t: f64) -> Option<[Vec3; 4]> {
+    let degree = knots.len().checked_sub(points.len() + 1)?;
+    let (first, ders) = solidcraft_geom::bspline::basis_derivatives(knots, degree, points.len(), t, 3)?;
+    let (mut a, mut w) = ([Vec3::ZERO; 4], [0.0f64; 4]);
+    for (i, d) in ders.iter().enumerate() {
+        let h = points.get(first + i)?;
+        for k in 0..4 {
+            let n = d.get(k).copied().unwrap_or(0.0);
+            a[k] = a[k] + Vec3::new(h[0], h[1], h[2]) * n;
+            w[k] += h[3] * n;
+        }
+    }
+    if !(w[0].abs() > 0.0) {
+        return None;
+    }
+    let mut c = [Vec3::ZERO; 4];
+    for k in 0..4 {
+        let mut num = a[k];
+        for j in 1..=k {
+            let binom = [[1.0, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0], [1.0, 2.0, 1.0, 0.0], [1.0, 3.0, 3.0, 1.0]][k][j];
+            num = num - c[k - j] * (binom * w[j]);
+        }
+        c[k] = num * (1.0 / w[0]);
+    }
+    Some(c)
+}
+
+/// A complex number as `(re, im)`.
+type Cx = (f64, f64);
+
+fn cmul(p: Cx, q: Cx) -> Cx {
+    (p.0 * q.0 - p.1 * q.1, p.0 * q.1 + p.1 * q.0)
+}
+
+fn cdiv(p: Cx, q: Cx) -> Cx {
+    let n = q.0 * q.0 + q.1 * q.1;
+    ((p.0 * q.0 + p.1 * q.1) / n, (p.1 * q.0 - p.0 * q.1) / n)
+}
+
+/// The Bézier control points of the degree-7 Hermite piece on `[s0, s1]` with the values and
+/// first three derivatives `p` at its start and `r` at its end (the k-th derivative of a Bézier of
+/// degree 7 at an end is 7!/(7 − k)! times the k-th difference of the control points there,
+/// scaled by `(s1 − s0)^k`).
+fn hermite7(p: [[f64; 2]; 4], r: [[f64; 2]; 4], h: f64) -> [[f64; 2]; 8] {
+    let (k1, k2, k3) = (h / 7.0, h * h / 42.0, h * h * h / 210.0);
+    let at = |d: [[f64; 2]; 4], m: [f64; 4]| [0, 1].map(|i| d[0][i] * m[0] + d[1][i] * m[1] + d[2][i] * m[2] + d[3][i] * m[3]);
+    [
+        p[0],
+        at(p, [1.0, k1, 0.0, 0.0]),
+        at(p, [1.0, 2.0 * k1, k2, 0.0]),
+        at(p, [1.0, 3.0 * k1, 3.0 * k2, k3]),
+        at(r, [1.0, -3.0 * k1, 3.0 * k2, -k3]),
+        at(r, [1.0, -2.0 * k1, k2, 0.0]),
+        at(r, [1.0, -k1, 0.0, 0.0]),
+        r[0],
+    ]
+}
+
+/// The pcurve on a cylinder or cone face of a B-spline edge with no closed form there: a
+/// flat-pattern edge folded across a bend, a helix, which no NURBS is (its 3D curve is already
+/// degree-7 Hermite pieces of the fold). The same construction in (u, v): degree-7 Hermite pieces
+/// on the edge's own parameter, matching the point and first three derivatives of its angle about
+/// the axis (through the form's quarter parameter when `form` is given) and of its height along
+/// the axis at each piece's ends. The pieces break where the angle crosses a quarter (the form's
+/// parameter is only once differentiable there) and are halved until they lift onto the edge
+/// within [`FOLD_TOL`]. An approximation by construction from exact derivatives, like the edge
+/// itself; nothing is fitted. `None` for another face or curve, or past the piece budget.
+fn folded_pcurve(label: &b::Surface, form: Option<&b::Surface>, curve: &b::Curve, [t0, t1]: [f64; 2]) -> Option<Curve2> {
+    use solidcraft_geom::bspline::{HPoint, bezier_point};
+    use std::f64::consts::{FRAC_PI_2, FRAC_PI_8, TAU};
+    let (b::Surface::Cylinder { frame, .. } | b::Surface::Cone { frame, .. }) = label else { return None };
+    let (knots, points): (&[f64], Vec<[f64; 4]>) = match curve {
+        b::Curve::Spline { knots, points } => (knots, points.iter().map(|p| [p[0], p[1], p[2], 1.0]).collect()),
+        b::Curve::RationalSpline { knots, homogeneous } => (knots, homogeneous.clone()),
+        _ => return None,
+    };
+    let points = points.as_slice();
+    let (o, x, y) = (v3(frame.origin), v3(frame.x_axis).normalized()?, v3(frame.y_axis).normalized()?);
+    let z = x.cross(y).normalized()?;
+    let on = SurfaceEval::new(form.unwrap_or(label)).ok()?;
+    let edge = solidcraft_geom::brep::eval::CurveEval::new(curve).ok()?;
+    let raw_angle = |t: f64| -> Option<f64> {
+        let d = curve_derivs(knots, points, t)?[0] - o;
+        Some(d.dot(y).atan2(d.dot(x)))
+    };
+    let angle_near = |t: f64, near: f64| raw_angle(t).map(|a| a + ((near - a) / TAU).round() * TAU);
+    // The angle's branch at the start: inside the form's u range when there is one.
+    let mut start = raw_angle(t0)?;
+    if let Some(b::Surface::RationalSpline { u_knots, .. }) = form {
+        let (lo, hi) = (*u_knots.first()?, *u_knots.last()?);
+        let inside = |a: f64| (lo - 1e-9..=hi + 1e-9).contains(&solidcraft_geom::brep::form::param(a));
+        start = (-2..=2).map(|k| start + f64::from(k) * TAU).find(|a| inside(*a))?;
+    }
+    // The angle along the edge, unwrapped from its start, on a grid fine against the edge's spans.
+    let grid = 64 * knots.len().max(1);
+    let ts: Vec<f64> = (0..=grid).map(|i| t0 + (t1 - t0) * i as f64 / grid as f64).collect();
+    let mut angles = Vec::with_capacity(ts.len());
+    let mut prev = start;
+    for &t in &ts {
+        prev = angle_near(t, prev)?;
+        angles.push(prev);
+    }
+    // Breaks: the ends, and on a form each crossing of a quarter, found by bisection.
+    let mut breaks = vec![(t0, *angles.first()?)];
+    if form.is_some() {
+        for i in 1..ts.len() {
+            let (a0, a1) = (*angles.get(i - 1)?, *angles.get(i)?);
+            let (q0, q1) = ((a0 / FRAC_PI_2).floor(), (a1 / FRAC_PI_2).floor());
+            if q0 == q1 {
+                continue;
+            }
+            let target = q0.max(q1) * FRAC_PI_2;
+            let (mut lo, mut hi) = (*ts.get(i - 1)?, *ts.get(i)?);
+            for _ in 0..80 {
+                let mid = 0.5 * (lo + hi);
+                if (angle_near(mid, a0)? - target).signum() == (a0 - target).signum() {
+                    lo = mid;
+                } else {
+                    hi = mid;
+                }
+            }
+            breaks.push((0.5 * (lo + hi), target));
+        }
+    }
+    breaks.push((t1, *angles.last()?));
+    // A crossing within rounding of another break (or of an end) is that break.
+    let close = 1e-12 * (t1 - t0).abs().max(1.0);
+    let mut kept: Vec<(f64, f64)> = Vec::with_capacity(breaks.len());
+    for (i, b) in breaks.iter().enumerate() {
+        let is_end = i + 1 == breaks.len();
+        match kept.last_mut() {
+            Some(last) if (b.0 - last.0).abs() <= close => {
+                if is_end {
+                    *last = *b;
+                }
+            }
+            _ => kept.push(*b),
+        }
+    }
+    let breaks = kept;
+    // u and v with their first three derivatives at `t`, angle `a`, in quarter `q` of the form.
+    let tq = FRAC_PI_8.tan();
+    let uv = |t: f64, a: f64, q: f64| -> Option<[[f64; 2]; 4]> {
+        let c = curve_derivs(knots, points, t)?;
+        let w = |d: Vec3| (d.dot(x), d.dot(y));
+        let (w0, r1, r2, r3) = (w(c[0] - o), w(c[1]), w(c[2]), w(c[3]));
+        let (r1, r2, r3) = (cdiv(r1, w0), cdiv(r2, w0), cdiv(r3, w0));
+        // θ = Im log w: (log w)' = w'/w, (log w)'' = w''/w − (w'/w)², (log w)''' = w'''/w −
+        // 3 w''w'/w² + 2 (w'/w)³.
+        let a1 = r1.1;
+        let a2 = r2.1 - cmul(r1, r1).1;
+        let a3 = r3.1 - 3.0 * cmul(r2, r1).1 + 2.0 * cmul(cmul(r1, r1), r1).1;
+        let h = [(c[0] - o).dot(z), c[1].dot(z), c[2].dot(z), c[3].dot(z)];
+        let u = if form.is_some() {
+            // The quarter parameter q + ½(1 + g/tan(π/8)), g = tan(ψ/2 − π/8), and its derivatives
+            // in θ, through to t by the chain rule.
+            let g = (0.5 * (a - q * FRAC_PI_2) - FRAC_PI_8).tan();
+            let g1 = 0.5 * (1.0 + g * g);
+            let g2 = g * g1;
+            let g3 = g1 * g1 + g * g2;
+            let (f1, f2, f3) = (g1 / (2.0 * tq), g2 / (2.0 * tq), g3 / (2.0 * tq));
+            [q + 0.5 * (1.0 + g / tq), f1 * a1, f2 * a1 * a1 + f1 * a2, f3 * a1 * a1 * a1 + 3.0 * f2 * a1 * a2 + f1 * a3]
+        } else {
+            [a, a1, a2, a3]
+        };
+        Some([0, 1, 2, 3].map(|k| [u[k], h[k]]))
+    };
+    let size = points.iter().fold(1.0f64, |m, h| [h[0], h[1], h[2]].iter().fold(m, |m, c| m.max((c / h[3]).abs())));
+    let mut pieces = 1usize;
+    while pieces * breaks.len() <= FOLD_PIECES {
+        let mut knots2 = vec![t0; 8];
+        let mut ctrl: Vec<[f64; 2]> = Vec::new();
+        let mut worst = 0.0f64;
+        for w in breaks.windows(2) {
+            let (&(a, aa), &(e, ea)) = (w.first()?, w.get(1)?);
+            let q = (0.5 * (aa + ea) / FRAC_PI_2).floor();
+            let mut last = aa;
+            for k in 0..pieces {
+                let (s0, s1) = (a + (e - a) * k as f64 / pieces as f64, a + (e - a) * (k + 1) as f64 / pieces as f64);
+                let (b0, b1) = (last, if k + 1 == pieces { ea } else { angle_near(s1, last)? });
+                last = b1;
+                let piece = hermite7(uv(s0, b0, q)?, uv(s1, b1, q)?, s1 - s0);
+                ctrl.extend_from_slice(if ctrl.is_empty() { &piece[..] } else { &piece[1..] });
+                knots2.extend(std::iter::repeat_n(s1, 7));
+                let h = piece.map(|c| HPoint { x: c[0], y: c[1], z: 0.0, w: 1.0 });
+                for j in 1..8 {
+                    let f = f64::from(j) / 8.0;
+                    let q2 = bezier_point(&h, f).project()?;
+                    worst = worst.max(on.point([q2.x, q2.y])?.dist(edge.point(s0 + (s1 - s0) * f)?));
+                }
+            }
+        }
+        knots2.push(t1);
+        if worst <= FOLD_TOL * size {
+            return Some(Curve2::Spline { knots: knots2, points: ctrl });
+        }
+        pieces *= 2;
+    }
+    None
+}
+
 /// The exact pcurve of an edge of a ruled face `(1 − v)·A(u) + v·B(u)` (`rows` its chains A and
 /// B, of `degree` on `knots`): a chain lying on A or B either way round is `(u(t), 0 or 1)` with
 /// `u` affine in its parameter, and a line from a corner of A to the corner of B at the same `u`
@@ -1191,16 +1401,29 @@ pub(crate) fn assemble(parts: Parts, solid: bool) -> Result<Assembled> {
     let mut face_forms: Vec<Option<b::Surface>> = Vec::with_capacity(nf);
     for (fi, (f, ls)) in faces.iter().zip(&loops).enumerate() {
         let failed = |e: pcurve::Error| KernelError::Failed(format!("sew: face {fi}: {e}"));
+        // Edges folded across a bend: B-spline edges with no given pcurve on a cylinder or cone.
+        let folded_on = |form: Option<&b::Surface>| -> Vec<(u32, Curve2)> {
+            ls.iter()
+                .flatten()
+                .filter(|(_, _, given)| given.is_none())
+                .filter_map(|&(e, _, _)| {
+                    let edge = edges.get(e as usize)?;
+                    Some((e, folded_pcurve(&f.surface, form, &edge.curve, edge.interval)?))
+                })
+                .collect()
+        };
+        let folded_label = folded_on(None);
         let form = match f.surface {
             b::Surface::Cylinder { .. } | b::Surface::Cone { .. } | b::Surface::Sphere { .. } | b::Surface::Torus { .. } => {
-                let label = pcurve::pcurves(&f.surface, &uses_of(&edges, &forms, ls, false, tol)?, tol).map_err(failed)?;
+                let label = pcurve::pcurves(&f.surface, &uses_of(&edges, &forms, ls, false, &folded_label, tol)?, tol).map_err(failed)?;
                 let extent = domain_of(&f.surface, &label).ok_or_else(|| fail(format!("face {fi} has no area in its parameters")))?;
                 Some(b::form::surface(&f.surface, extent).ok_or_else(|| fail(format!("face {fi}: more than one turn")))?.0)
             }
             _ => None,
         };
         let on = form.as_ref().unwrap_or(&f.surface);
-        let pc = pcurve::pcurves(on, &uses_of(&edges, &forms, ls, true, tol)?, tol).map_err(failed)?;
+        let folded_form = if form.is_some() { folded_on(form.as_ref()) } else { Vec::new() };
+        let pc = pcurve::pcurves(on, &uses_of(&edges, &forms, ls, true, &folded_form, tol)?, tol).map_err(failed)?;
         domains.push(domain_of(on, &pc).ok_or_else(|| fail(format!("face {fi} has no area in its parameters")))?);
         pcs.push(pc);
         face_forms.push(form);
@@ -1353,12 +1576,14 @@ pub(crate) fn assemble(parts: Parts, solid: bool) -> Result<Assembled> {
 /// held inside the surface's natural domain (one turn where periodic). None when it has no
 /// width in either direction.
 /// A face's edge uses for `pcurve::pcurves`: on each edge's NURBS form where it has one and
-/// `formed`, else on its curve as given.
+/// `formed`, else on its curve as given; `folded` gives the pcurves of edges folded across the
+/// face (see [`folded_pcurve`]).
 fn uses_of<'a>(
     edges: &'a [PartEdge],
     forms: &'a [Option<(b::Curve, [f64; 2])>],
     ls: &[Vec<(u32, bool, Option<&'a Curve2>)>],
     formed: bool,
+    folded: &'a [(u32, Curve2)],
     tol: f64,
 ) -> Result<Vec<Vec<pcurve::EdgeUse<'a>>>> {
     ls.iter()
@@ -1370,6 +1595,7 @@ fn uses_of<'a>(
                         Some((c, i)) if formed => (c, *i),
                         _ => (&edge.curve, edge.interval),
                     };
+                    let given = given.or_else(|| folded.iter().find(|(f, _)| *f == e).map(|(_, c)| c));
                     Ok(pcurve::EdgeUse { edge: e, curve, interval, forward, tolerance: tol, given })
                 })
                 .collect()