SMLib import: an extruded circle's circle comes from its widest span, not the face's range

Since 6cc759750 an extruded circle is a cylinder over the face's trimmed range, and the circle was
fitted through three points of that range. On the shallow strip of the extruded-ellipse test
(radius 3000, 6 mm wide: 0.002 rad round) the three points are 2 mm apart and the fit is
ill-conditioned: the centre came 2.5e-10 mm off the circle's plane (the frame's x axis tilted
8.3e-14 out of it) and the radius 8.1e-9 mm short. SMLib took that cylinder, but on the way back
(the import's SMLib round trip, then Body::smlib) SmCylinder::IsNurbSurfaceCylinder did not find
a cylinder in SMLib's own NURBS of it, and the strip was refused as needing repair
(SM_ERR_INVALID_INPUT, AnalyticOfNurbs). Before 6cc759750 the fit ran over the record's domain,
1.4e-6 rad wider on each side, and happened to come out level.

The patch's trims were right all along: they straddle the profile's seam (u = 4, the knot
between the spline's two quarter spans) without wrapping, [3.9992930395761980, 4.0007069604238020],
which is the exact 0.001 rad either side of it.

The circle is now fixed by three points on the arc's widest knot span (a quarter turn here),
where it is determined to rounding, and checked as before at 65 points over the face's range.
A revolution's profile circle is fitted as before.

Regression test (smlib import tests): a 0.002 rad patch of an extruded circle of radius 3000
across that knot is the cylinder of radius 3000 to 1e-10 (it was 2.8e-10 short), its axis
along z to 1e-15, over exactly the patch's range.

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

crates/smlib/src/import/analytic.rs

@@ -105,7 +105,10 @@ fn extruded_circle(u_knots: &[f64], v_knots: &[f64], net: &[Vec<[f64; 4]>], trim
         return None;
     }
     let arc = CurveEval::new(&b::Curve::RationalSpline { knots: arc_knots.to_vec(), homogeneous: first.clone() }).ok()?;
-    let (center, radius, normal) = circle(&arc, arc_range)?;
+    // The circle from the arc's widest span, not from the face's range: three points on a
+    // face a small part of a turn round (6 mm of a 3000 mm radius) put the centre off the
+    // circle's plane and the radius off by 1e-8 mm, which SMLib does not take back as a cylinder.
+    let (center, radius, normal) = circle(&arc, widest_span(arc_knots, first.len())?, arc_range)?;
     let z = d.normalized()?;
     if normal.cross(z).len() > 1e-9 {
         return None;
@@ -397,7 +400,7 @@ impl Revolution {
             return Ok((cone, [angle, v]));
         }
         // A circular profile.
-        let (center, radius, normal) = circle(&self.profile, t).ok_or(NOT_ANALYTIC)?;
+        let (center, radius, normal) = circle(&self.profile, t, t).ok_or(NOT_ANALYTIC)?;
         let slack = 1e-9 * radius.max(1.0);
         if normal.dot(self.z).abs() > 1e-9 || normal.dot(x).abs() > 1e-9 || (center - self.o).dot(y).abs() > slack {
             return Err(NOT_ANALYTIC);
@@ -506,11 +509,20 @@ fn agree((record, domain): (&b::Surface, [[f64; 2]; 2]), (analytic, range): (&b:
     Ok(d < 0.0)
 }
 
-/// The circle a curve runs along over `t`: centre, radius and unit normal (its turning sense),
-/// checked at 65 points within 10⁻⁹ of the radius (at least 1 mm).
-fn circle(c: &CurveEval, t: [f64; 2]) -> Option<(Vec3, f64, Vec3)> {
+/// The widest knot span of a B-spline with `points` control points on `knots`.
+fn widest_span(knots: &[f64], points: usize) -> Option<[f64; 2]> {
+    let degree = knots.len().checked_sub(points + 1)?;
+    knots.get(degree..=points)?.windows(2).map(|w| [w[0], w[1]]).max_by(|a, b| (a[1] - a[0]).total_cmp(&(b[1] - b[0])))
+}
+
+/// The circle a curve runs along over `t`, through its points at `fit` (a range of it on that
+/// circle a good part of a turn round, where three points fix the circle to rounding): centre,
+/// radius and unit normal (its turning sense over `fit`), checked at 65 points over `t` within
+/// 10⁻⁹ of the radius (at least 1 mm).
+fn circle(c: &CurveEval, fit: [f64; 2], t: [f64; 2]) -> Option<(Vec3, f64, Vec3)> {
+    let on = |s: f64| c.point(fit[0] + (fit[1] - fit[0]) * s);
     let at = |s: f64| c.point(t[0] + (t[1] - t[0]) * s);
-    let (p, q, r) = (at(0.0)?, at(1.0 / 3.0)?, at(2.0 / 3.0)?);
+    let (p, q, r) = (on(0.0)?, on(1.0 / 3.0)?, on(2.0 / 3.0)?);
     let (u, w) = (q - p, r - p);
     let normal = u.cross(w);
     let nn = normal.dot(normal);

crates/smlib/src/import/tests.rs

@@ -267,3 +267,54 @@ fn only_a_tube_clear_of_its_axis_is_a_ring_torus() {
         assert!(matches!(ring(major), Surface::PlacedRevolution { .. }), "major {major}");
     }
 }
+
+/// A 6 mm wide patch of an extruded circle of radius 3000 whose trims straddle the knot where
+/// the profile's seam lies (u = 4 of a half circle in two quarter spans): the cylinder takes the
+/// patch's own range, 0.002 rad round and 7 high, on the circle exactly. Fitted on the patch's
+/// range (three points 2 mm apart) the centre came off the circle's plane by 2.5e-10 mm and the
+/// radius off by 8e-9 mm, a cylinder SMLib did not take back from its own NURBS of it.
+#[test]
+fn a_patch_across_an_extruded_circles_seam_is_its_exact_cylinder() {
+    let (r, h, half) = (3000.0, 7.0, 0.001);
+    let c = [10.0 - r, 0.0, 0.0];
+    let w = std::f64::consts::FRAC_1_SQRT_2;
+    let profile = [([c[0], -r], 1.0), ([c[0] + r, -r], w), ([c[0] + r, 0.0], 1.0), ([c[0] + r, r], w), ([c[0], r], 1.0)];
+    let homogeneous = profile.iter().map(|&([x, y], wt)| vec![[x * wt, y * wt, 0.0, wt], [x * wt, y * wt, h * wt, wt]]).collect();
+    let surface = Surface::RationalSpline { u_knots: vec![3.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 5.0], v_knots: vec![0.0, 0.0, h, h], homogeneous };
+    // The quarter span after the seam reaches `half` at quarter_param(half); the one before it
+    // ends there, so the trims run symmetrically round the knot.
+    let (u0, u1) = (4.0 - quarter_param(half), 4.0 + quarter_param(half));
+    let corners = [[u0, 0.0], [u1, 0.0], [u1, h], [u0, h]];
+    let edges: Vec<OrientedEdge> = (0..4)
+        .map(|k| OrientedEdge {
+            edge: k as u32,
+            forward: true,
+            pcurve: Curve2::Spline { knots: vec![0.0, 0.0, 1.0, 1.0], points: vec![corners[k], corners[(k + 1) % 4]] },
+        })
+        .collect();
+    let line = Curve::Line { start: [0.0; 3], end: [1.0, 0.0, 0.0] };
+    let record = Record {
+        version: b::VERSION,
+        tolerance: 1e-6,
+        regions: Vec::new(),
+        shells: Vec::new(),
+        faces: vec![Face { surface, nurbs: None, domain: [[3.0, 5.0], [0.0, h]], loops: vec![Loop::Edges { edges }] }],
+        edges: (0..4)
+            .map(|_| Edge { nurbs: None, curve: line.clone(), interval: [0.0, 1.0], vertices: [0, 0], tolerance: 1e-6, uses: Vec::new() })
+            .collect(),
+        vertices: Vec::new(),
+    };
+    let made = analytic::analytic(&record).unwrap();
+    let f = &made.faces[0];
+    let Surface::Cylinder { frame, radius } = &f.surface else { panic!("{:?}", f.surface) };
+    assert!((radius - r).abs() <= 1e-10, "radius {radius}");
+    let z = cross(frame.x_axis, frame.y_axis);
+    assert!(z[0].abs() <= 1e-15 && z[1].abs() <= 1e-15 && frame.x_axis[2].abs() <= 1e-15, "axis {z:?}, x {:?}", frame.x_axis);
+    assert!(frame.origin[2].abs() <= 1e-12 && (frame.origin[0] - c[0]).abs() <= 1e-9 && frame.origin[1].abs() <= 1e-9, "centre {:?}", frame.origin);
+    // u from the patch's first edge, 2 * half round; v its height.
+    let [[a0, a1], [v0, v1]] = f.domain;
+    assert!(a0 == 0.0 && (a1 - 2.0 * half).abs() <= 1e-12, "{:?}", f.domain);
+    assert!(v0.abs() <= 1e-12 && (v1 - h).abs() <= 1e-12, "{:?}", f.domain);
+    let start = (-half).sin_cos();
+    assert!((frame.x_axis[0] - start.1).abs() <= 1e-12 && (frame.x_axis[1] - start.0).abs() <= 1e-12, "x {:?}", frame.x_axis);
+}