diff --git a/crates/smlib/src/import/analytic.rs b/crates/smlib/src/import/analytic.rs index ce7a61af27..205d61e49a 100644 --- a/crates/smlib/src/import/analytic.rs +++ b/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); diff --git a/crates/smlib/src/import/tests.rs b/crates/smlib/src/import/tests.rs index 553a18d1e0..cf96de0d29 100644 --- a/crates/smlib/src/import/tests.rs +++ b/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 = (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); +}