geom: the patch fill's cross-boundary targets agree with its corners

A tensor patch's cross derivative along a side is, at a corner, the tangent of the side it
meets there (the two share the corner's next control point): the compatibility condition of
Coons and Gregory interpolation. The weak rows asked t·c′ = 0 and w·c′ = s right up to every
corner, which no patch can meet where the next side is not square to the edge; at a smooth
vertex (a tube's end in four quarter arcs) the next side runs along the edge itself. The least
squares spread the contradiction into the first spans and folded them: the tube end's worst
normal went 13.8°, 17.9°, 89.8° (folded), 85.6°, ... 25.1° (folded) with refinement, and the
symmetric loop came back asymmetric.

The target is now the Coons start's own cross derivative, compatible at every corner by
construction, turned about the edge into the neighbour's tangent plane (its part along the
edge kept, the rest laid along the leave direction at its length). The tube end converges
21°, 3.8°, 0.76°, 0.18°, 0.045°, symmetric, and does not fold.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/geom/src/fill_iso.rs

@@ -44,10 +44,16 @@
 //!   formulation ([`crate::linkage`]): at stations along an edge, the iso-parameter curve `c`
 //!   of the fill leaving across it (σ growing into the fill) must have its height above the
 //!   neighbour vanish to order n:
-//!   - order 1: `n·c′ = 0` (`n` the neighbour's normal), with weak rows `w·c′ = s`, `t·c′ = 0`
-//!     (`w` the neighbour's leave direction, `t = n × w` along the edge, `s` how fast the Coons
-//!     start leaves along `w`) that keep the fill leaving into itself, about as fast; they act
-//!     in the tangent plane only, so they never pull against the G1 row;
+//!   - order 1: `n·c′ = 0` (`n` the neighbour's normal), with weak rows `w·c′ = s`, `t·c′ = a`
+//!     (`w` the neighbour's leave direction, `t = n × w` along the edge) that keep the fill
+//!     leaving into itself, about as fast: `(a, s)` is the Coons start's own cross derivative
+//!     `d` turned about the edge into the neighbour's tangent plane (`a = t·d`, `s` the length of
+//!     the rest of `d`). The Coons start is compatible at every corner (its cross derivative there
+//!     is the tangent of the side it meets, as any tensor patch's must be: the two share the
+//!     corner's next control point), so the target is too wherever that side lies in the
+//!     neighbour's tangent plane; a target that is not (`t·c′ = 0` up to the corner) is a
+//!     contradiction the least squares can only fold the corner's spans over. They act in the
+//!     tangent plane only, so they never pull against the G1 row;
 //!   - order k = 2 … n: `n·c⁽ᵏ⁾ = k!·[σᵏ] Σ_m H_m(x(σ))/m!` ([`crate::linkage::Jet::targets`]):
 //!     the normal curvature (G2), its rate of change (G3) and that rate's (G4) are the
 //!     neighbour's, in the fill's own leave direction. The right-hand side reads `c′ … c⁽ᵏ⁻¹⁾`
@@ -813,7 +819,10 @@ struct Station {
     tilt: f64,
     held: bool,
     clear: bool,
+    /// The weak rows' target for the cross derivative `c′` ([`stations`]): `w·c′ = speed`,
+    /// `t·c′ = slide`.
     speed: f64,
+    slide: f64,
 }
 
 impl Station {
@@ -1180,13 +1189,28 @@ fn stations(
             let ask = Asked::near(asked, [(corners.get(k).copied(), f), (corners.get((k + 1) % n).copied(), 1.0 - f)]);
             let clear = f > CORNER_SHARE && f < 1.0 - CORNER_SHARE && ask.full == asked;
             let (uv, (du, dv, sign)) = station_on(side, t);
-            // How fast the Coons start leaves into the leave direction: it blends the sides
-            // across, so it leaves as fast as the opposite sides lie apart; the side's own speed
-            // where it does not leave into the patch at all.
-            let along = p.eval(uv[0], uv[1], dv, du).map(|d| d.len()).unwrap_or(0.0);
-            let speed = p.eval(uv[0], uv[1], du, dv).map(|d| (d * sign).dot(c.leave)).ok().filter(|x| *x > 0.0).unwrap_or(along);
+            // The Coons start's cross derivative turned about the edge into the neighbour's
+            // tangent plane: its part along the edge kept, the rest laid along the leave
+            // direction at its own length. At a corner it is the next side's tangent, which the
+            // patch's cross derivative is there whatever its inside (module docs).
+            let d = p.eval(uv[0], uv[1], du, dv)? * sign;
+            let slide = d.dot(c.normal.cross(c.leave));
+            let speed = (d.dot(d) - slide * slide).max(0.0).sqrt();
             let tilt = edges.get(k).and_then(|e| e.tangent(param)).map_or(0.0, |d| d.dot(c.normal).abs().clamp(0.0, 1.0).asin());
-            out.push(Station { patch, side, t, across: c, order: ask.order(), full: ask.full, fade: ask.fade, tilt, held: ask.held, clear, speed });
+            out.push(Station {
+                patch,
+                side,
+                t,
+                across: c,
+                order: ask.order(),
+                full: ask.full,
+                fade: ask.fade,
+                tilt,
+                held: ask.held,
+                clear,
+                speed,
+                slide,
+            });
         }
     }
     Ok(out)
@@ -1296,7 +1320,7 @@ fn knot_line_rows(quilt: &Quilt, angle_tol: f64) -> SplineResult<Vec<Row>> {
 }
 
 /// A station's rows (module docs), linearised about the patch as it stands:
-/// - G1: `n·c′ = 0`, and the weak Hermite rows `w·c′ = s`, `t·c′ = 0` in the neighbour's tangent
+/// - G1: `n·c′ = 0`, and the weak Hermite rows `w·c′ = s`, `t·c′ = a` in the neighbour's tangent
 ///   plane only (`t = n × w`), so they never pull against the G1 row;
 /// - orders k = 2 … n: `n·c⁽ᵏ⁾ = T_k(x₁ … x_{k−1})` ([`crate::linkage::Jet::linearised`]) by
 ///   Newton's method: `T_k` is replaced by its tangent at the last solution's `x⁰`, so the row is
@@ -1318,7 +1342,7 @@ fn linkage_rows(rows: &mut Vec<Row>, patch: &Patch, s: &Station, span: f64) -> S
         .collect::<SplineResult<_>>()?;
     let Some(first) = derivative.first() else { return Ok(()) };
     rows.push(Row { terms: along(s.patch, first, n), value: 0.0, weight: 1.0 });
-    for (dir, value) in [(w, s.speed), (t, 0.0)] {
+    for (dir, value) in [(w, s.speed), (t, s.slide)] {
         rows.push(Row { terms: along(s.patch, first, dir), value, weight: HERMITE_WEIGHT });
     }
     let Some(jet) = s.across.jet else { return Ok(()) };

crates/geom/src/fill_iso_tests.rs

@@ -217,3 +217,65 @@ fn a_rational_patch_s_derivatives_are_its_differences() {
         assert!(exact.dist(differences) < 1e-5 * exact.len().max(1.0), "∂({du}, {dv}): {exact:?} vs {differences:?}");
     }
 }
+
+/// The four quarter arcs of the circle of radius `r` about z at height `z`, on their NURBS form.
+fn circle_quarters(z: f64, r: f64) -> Vec<IsoEdge> {
+    let f = frame(Vec3::Z * z, Vec3::X, Vec3::Y);
+    let q = std::f64::consts::FRAC_PI_2;
+    (0..4).map(|k| arc(f.clone(), r, [k as f64 * q, (k + 1) as f64 * q], false)).collect()
+}
+
+/// Neighbours given by a normal and a leave direction at each point of the loop.
+fn beside_by<'a>(
+    edges: &'a [IsoEdge],
+    frame_at: impl Fn(Vec3) -> (Vec3, Vec3) + 'a,
+) -> impl FnMut(usize, &[f64]) -> Result<Vec<Option<Across>>, String> + 'a {
+    move |k, ts| {
+        ts.iter()
+            .map(|t| {
+                let p = eval_edge_h(&edges[k], *t).map_err(|e| e.to_string())?.project().ok_or("a point at infinity")?;
+                let (normal, leave) = frame_at(p);
+                Ok(Some(Across { leave, normal, bend: None, jet: None }))
+            })
+            .collect()
+    }
+}
+
+/// Whether any patch of the fill folds (its normal reverses between neighbouring samples).
+fn any_fold(fill: &IsoFill) -> bool {
+    fill.patches.iter().any(|p| {
+        let s = SurfaceEval::new(&Surface::RationalSpline {
+            u_knots: p.surface.knots()[0].clone(),
+            v_knots: p.surface.knots()[1].clone(),
+            homogeneous: (0..p.surface.count()[0])
+                .map(|i| (0..p.surface.count()[1]).map(|j| p.surface.point(i, j)).map(|h| [h.x, h.y, h.z, h.w]).collect())
+                .collect(),
+        })
+        .unwrap();
+        let n = 48;
+        let normal = |i: usize, j: usize| {
+            let (_, su, sv) = s.derivs(i as f64 / n as f64, j as f64 / n as f64).unwrap();
+            su.cross(sv).normalized()
+        };
+        (0..=n).any(|i| {
+            (0..=n).any(|j| {
+                let Some(a) = normal(i, j) else { return false };
+                [(i + 1, j), (i, j + 1)].iter().any(|&(k, l)| k <= n && l <= n && normal(k, l).is_some_and(|b| a.dot(b) <= 0.0))
+            })
+        })
+    })
+}
+
+/// A tube's end: the circle in four quarter arcs, G1 to the tube's wall (its normal radial, the
+/// wall going on upwards). Every corner of the loop is a smooth vertex, where the cross derivative
+/// of one side is the next side's tangent, along the edge itself.
+#[test]
+fn a_tube_end_capped_g1_holds_the_wall_and_does_not_fold() {
+    let edges = circle_quarters(30.0, 10.0);
+    let radial = |p: Vec3| (Vec3::new(p.x, p.y, 0.0).normalized().unwrap(), Vec3::Z);
+    let s = IsoSpec { tol: 2e-5, ..spec() };
+    let fill = iso_fill(&edges, &mut beside_by(&edges, radial), &s).unwrap();
+    eprintln!("tube end g1: {} patches, {:?}", fill.patches.len(), fill.achieved);
+    assert!(fill.achieved.normal_excess <= s.angle_tol);
+    assert!(!any_fold(&fill), "the cap folds");
+}