geom: patches about a fill's centre join G1 there, with the linear gluing function

Across each spoke the two patches' derivatives away from it were held opposed (parametric C1).
At the centre that asks the spokes' tangents d_{k-1} + d_{k+1} = 0, which for N other than a
multiple of four only d = 0 meets: every patch got a flat, singular corner at the centre (a
pentagon on a wavy rim: derivatives of 3.3e-5 there; the hexagonal rim's control net doubled
back and SMLib's check refused it as NonMonotonic).

The join is now G1 with the linear gluing function of Loop, Peters and Hahn:
d/du S_k + d/dv S_{k+1} = alpha(t) d/dt S along the spoke, alpha(t) = 2 cos(2 pi / N) t. At the
midpoint alpha is 0 (the edge's two halves leave it at the same speed, opposed); at the centre
it is the value any affine image of the regular N-star meets, d_{k-1} + d_{k+1} = 2 cos(2 pi/N)
d_k. N = 4 is the old row. The pentagon's centre corners now open 72° with derivatives of 16.2
(N = 5 to 8: 360°/N each).

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

crates/geom/src/fill_iso.rs

@@ -61,8 +61,15 @@
 //!   - across a patch's own C0 knot lines (a side's break where it is only C0 in its
 //!     parameter), the derivatives either side in the ratio of the crossing sides' speeds
 //!     ([`knot_line_rows`]): the fill is tangent continuous there too;
-//!   - across each spoke, the two patches' derivatives away from it opposed (`∂S_k/∂u + ∂S_{k+1}/∂v
-//!     = 0` on the shared parameter): the patches join tangent there;
+//!   - across each spoke, the two patches' derivatives away from it opposed up to the linear
+//!     gluing function of G1 joins at an N-valent point (C. Loop, J. Peters, J. Hahn):
+//!     `∂S_k/∂u + ∂S_{k+1}/∂v = α(t)·∂S/∂t` on the shared parameter t (0 at the edge's
+//!     midpoint, 1 at the centre), `α(t) = 2·cos(2π/N)·t`: the patches join tangent there. At
+//!     the midpoint `α = 0` (the edge's two halves leave it at the same speed, opposed); at the
+//!     centre it is the value the spokes' tangents `d_k` meet when they are any affine image of
+//!     the regular N-star, `d_{k−1} + d_{k+1} = 2·cos(2π/N)·d_k`. Parametric C1 (`α = 0`) asks
+//!     `d_{k−1} + d_{k+1} = 0` there, which for N other than a multiple of four only `d = 0`
+//!     meets: every patch would have a flat corner at the centre;
 //!   - points the fill passes through, located on the last solution (nearest point of any patch)
 //!     and fitted again (Hoschek's parameter correction).
 //!
@@ -1243,12 +1250,15 @@ fn solve(mut quilt: Quilt, stations: &[Station], spec: &IsoSpec) -> SplineResult
             // The span across the side: of v for sides 0 and 2, of u for 1 and 3.
             linkage_rows(&mut rows, patch, s, span(patch).get(1 - s.side % 2).copied().unwrap_or(1.0))?;
         }
-        // Across each spoke the two patches' derivatives away from it are opposed.
+        // Across each spoke the two patches' derivatives away from it are opposed, up to the
+        // gluing function's share of the derivative along it (module docs).
+        let gluing = 2.0 * (std::f64::consts::TAU / quilt.patches.len() as f64).cos();
         for &(left, right) in &quilt.spokes {
             let (Some(a), Some(b)) = (quilt.patches.get(left), quilt.patches.get(right)) else { continue };
             for i in 0..spoke_count {
                 let t = (i as f64 + 0.5) / spoke_count as f64;
-                on_each_axis(&mut rows, &[(left, a.row(1.0, t, 1, 0)?), (right, b.row(t, 1.0, 0, 1)?)], Vec3::ZERO, 1.0);
+                let along: Vec<(usize, usize, f64)> = a.row(1.0, t, 0, 1)?.into_iter().map(|(i, j, c)| (i, j, -gluing * t * c)).collect();
+                on_each_axis(&mut rows, &[(left, a.row(1.0, t, 1, 0)?), (right, b.row(t, 1.0, 0, 1)?), (left, along)], Vec3::ZERO, 1.0);
             }
         }
         for (p, (patch, uv)) in spec.points.iter().zip(&located) {

crates/geom/src/fill_iso_tests.rs

@@ -279,3 +279,51 @@ fn a_tube_end_capped_g1_holds_the_wall_and_does_not_fold() {
     assert!(fill.achieved.normal_excess <= s.angle_tol);
     assert!(!any_fold(&fill), "the cap folds");
 }
+
+/// The angle (degrees) between a patch's two derivatives at a corner of its square, and their
+/// lengths.
+fn corner_angle(p: &IsoPatch, uv: [f64; 2]) -> (f64, f64, f64) {
+    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 (_, su, sv) = s.derivs(uv[0], uv[1]).unwrap();
+    let cos = su.dot(sv) / (su.len() * sv.len());
+    (cos.clamp(-1.0, 1.0).acos().to_degrees(), su.len(), sv.len())
+}
+
+/// Loops of five to eight sides about a centre on a wavy rim: at the centre every
+/// patch has a proper corner (its two derivatives neither vanishing nor parallel). Parametric C1
+/// across the spokes cannot give that for N other than a multiple of four: it asks the spokes'
+/// tangents at the centre `d_{k−1} + d_{k+1} = 0`, which only `d = 0` meets.
+#[test]
+fn patches_about_a_centre_have_proper_corners_there() {
+    for n in 5..=8usize {
+        let corners: Vec<Vec3> = (0..n)
+            .map(|k| {
+                let a = 0.3 + k as f64 * std::f64::consts::TAU / n as f64;
+                Vec3::new(20.0 * a.cos(), 20.0 * a.sin(), 3.0 * (2.0 * a).sin())
+            })
+            .collect();
+        // Quadratic edges bowed up out of the rim's surface.
+        let edges: Vec<IsoEdge> = (0..n)
+            .map(|k| {
+                let (a, b) = (corners[k], corners[(k + 1) % n]);
+                IsoEdge { degree: 2, knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], points: vec![a, (a + b) * 0.5 + Vec3::Z * 2.0, b], weights: None }
+            })
+            .collect();
+        let mut g0 = |_k: usize, ts: &[f64]| -> Result<Vec<Option<Across>>, String> { Ok(vec![None; ts.len()]) };
+        let fill = iso_fill(&edges, &mut g0, &IsoSpec { tol: 2e-5, ..spec() }).unwrap();
+        let worst = fill
+            .patches
+            .iter()
+            .map(|p| corner_angle(p, [1.0, 1.0]))
+            .fold((180.0f64, f64::MAX), |(a, l), (b, x, y)| (a.min(b).min(180.0 - b), l.min(x).min(y)));
+        eprintln!("{n} sides: at the centre the patches' corners are at least {:.2}° from flat, derivatives at least {:.3e} long", worst.0, worst.1);
+        assert!(worst.0 > 0.5 && worst.1 > 1e-3, "{n} sides: a flat corner at the centre {worst:?}");
+    }
+}