diff --git a/crates/geom/src/fill_iso.rs b/crates/geom/src/fill_iso.rs index c318d50c0..b8da2d29e 100644 --- a/crates/geom/src/fill_iso.rs +++ b/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) { diff --git a/crates/geom/src/fill_iso_tests.rs b/crates/geom/src/fill_iso_tests.rs index 618870c56..84d9d0199 100644 --- a/crates/geom/src/fill_iso_tests.rs +++ b/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 = (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 = (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>, 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:?}"); + } +}