fix: STEP import reads unclamped (periodic) B-splines on the seam (#1916)

Thanks @NOWN-DEV.

crates/kernel/src/step_in/geom.rs

@@ -216,6 +216,57 @@ fn generated_knots(kind: &str, n_ctrl: usize, degree: usize) -> R<mt::KnotVec> {
     Ok(mt::KnotVec::from(v))
 }
 
+/// The ends of the part of a knot vector where a B-spline of `degree` is defined, when its end
+/// knots repeat fewer than `degree + 1` times (an unclamped or periodic knot vector, as
+/// OpenCascade writes closed B-splines and `UNIFORM_*` curves and surfaces have): only between
+/// knot `degree` and knot `len - degree - 1` do the basis functions sum to one. `None` for a
+/// clamped vector, whose whole range is the domain, and for one whose domain is a point.
+fn unclamped_ends(kv: &mt::KnotVec, degree: usize) -> Option<(Option<f64>, Option<f64>)> {
+    let n = kv.len();
+    if n < 2 * (degree + 1) || !(kv[n - 1 - degree] - kv[degree] > mt::TOLERANCE) {
+        return None;
+    }
+    let lo = (kv.multiplicity(0) <= degree).then(|| kv[degree]);
+    let hi = (kv.multiplicity(n - 1) <= degree).then(|| kv[n - 1 - degree]);
+    (lo.is_some() || hi.is_some()).then_some((lo, hi))
+}
+
+/// `c` restricted to its domain, with clamped ends (exact: knots are inserted at the domain's
+/// ends and the control points beyond them dropped). The kernel takes a B-spline's whole knot
+/// range as its parameter range, and an unclamped curve evaluated outside its domain leaves
+/// its control polygon (a rational one divides by zero at the outermost knots).
+fn clamped_curve<P: mt::control_point::ControlPoint<f64> + mt::Tolerance>(mut c: mt::BSplineCurve<P>) -> mt::BSplineCurve<P> {
+    let Some((lo, hi)) = unclamped_ends(c.knot_vec(), c.degree()) else { return c };
+    if let Some(lo) = lo {
+        c = c.cut(lo);
+    }
+    if let Some(hi) = hi {
+        let _ = c.cut(hi);
+    }
+    c
+}
+
+/// [`clamped_curve`] for a surface, in each direction whose knot vector is unclamped.
+fn clamped_surface<P: mt::control_point::ControlPoint<f64> + mt::Tolerance>(mut s: mt::BSplineSurface<P>) -> mt::BSplineSurface<P> {
+    if let Some((lo, hi)) = unclamped_ends(s.uknot_vec(), s.udegree()) {
+        if let Some(lo) = lo {
+            s = s.ucut(lo);
+        }
+        if let Some(hi) = hi {
+            let _ = s.ucut(hi);
+        }
+    }
+    if let Some((lo, hi)) = unclamped_ends(s.vknot_vec(), s.vdegree()) {
+        if let Some(lo) = lo {
+            s = s.vcut(lo);
+        }
+        if let Some(hi) = hi {
+            let _ = s.vcut(hi);
+        }
+    }
+    s
+}
+
 fn weights_list(p: Option<&Param>) -> R<Vec<f64>> {
     p.and_then(Param::as_list)
         .ok_or("rational B-spline without weights")?
@@ -267,9 +318,9 @@ fn bspline_curve(cx: &Ctx, e: &Entity) -> R<CurveGeo> {
         Some(r) => {
             let w = weights_list(r.first())?;
             let n = mt::NurbsCurve::<mt::Vector4>::try_from_bspline_and_weights(bsp, w).map_err(|e| format!("rational B-spline: {e}"))?;
-            Ok(CurveGeo::Nurbs(n))
+            Ok(CurveGeo::Nurbs(mt::NurbsCurve::new(clamped_curve(n.into_non_rationalized()))))
         }
-        None => Ok(CurveGeo::BSpline(bsp)),
+        None => Ok(CurveGeo::BSpline(clamped_curve(bsp))),
     }
 }
 
@@ -889,9 +940,9 @@ pub(crate) fn bspline_surface(cx: &Ctx, e: &Entity) -> R<mt::Surface> {
                 return Err("weights do not match the control net".into());
             }
             let n = mt::NurbsSurface::<mt::Vector4>::try_from_bspline_and_weights(bsp, w).map_err(|e| format!("rational surface: {e}"))?;
-            Ok(mt::Surface::NurbsSurface(n))
+            Ok(mt::Surface::NurbsSurface(mt::NurbsSurface::new(clamped_surface(n.into_non_rationalized()))))
         }
-        None => Ok(mt::Surface::BSplineSurface(bsp)),
+        None => Ok(mt::Surface::BSplineSurface(clamped_surface(bsp))),
     }
 }
 

crates/kernel/src/step_in/tests.rs

@@ -1255,3 +1255,120 @@ fn revolved_faces_skewed_never_panic_on_export() {
         }
     }
 }
+
+/// Weights of a rational quadratic full circle through nine control points.
+const CIRCLE_WEIGHTS: &str = "1.,0.707106781187,1.,0.707106781187,1.,0.707106781187,1.,0.707106781187,1.";
+/// Angle knots as OpenCascade writes a closed (periodic) quadratic circle: the end knots
+/// appear once, so the B-spline is defined from knot 2 to knot 9 only, -π to π.
+const PERIODIC_MULTS: &str = "1,2,2,2,2,2,1";
+const PERIODIC_KNOTS: &str = "-4.712388980385,-3.14159265359,-1.570796326795,0.,1.570796326795,3.14159265359,4.712388980385";
+
+/// Nine control points of a quadratic circle of radius `r` about (`x`, 0, `z`), from angle 0.
+fn circle_net(first: u64, x: f64, r: f64, z: f64) -> (String, Vec<String>) {
+    let (mut s, mut ids) = (String::new(), Vec::new());
+    for i in 0..9u64 {
+        let a = i as f64 * std::f64::consts::FRAC_PI_4;
+        let q = if i % 2 == 1 { r * std::f64::consts::SQRT_2 } else { r };
+        s += &format!("#{}=CARTESIAN_POINT('',({:?},{:?},{z:?}));\n", first + i, x + q * a.cos(), q * a.sin());
+        ids.push(format!("#{}", first + i));
+    }
+    (s, ids)
+}
+
+/// A tube face of radius `r` and length `h` about the Z axis through (`x`, 0, 0), written the
+/// way OpenCascade writes one: a rational degree (1, 2) B-spline surface with the periodic
+/// angle knots above, bounded by its two rims and a seam. Entities from `first` on; returns
+/// them and the face's id.
+fn periodic_knot_tube(first: u64, x: f64, r: f64, h: f64) -> (String, u64) {
+    let id = |k: u64| first + k;
+    let (net0, row0) = circle_net(id(40), x, r, 0.0);
+    let (net1, row1) = circle_net(id(50), x, r, h);
+    let (p0, p1, v0, v1) = (id(0), id(1), id(2), id(3));
+    let (c0, c1, up, down, ref_dir) = (id(4), id(5), id(6), id(7), id(8));
+    let (a0, a1, rim0, rim1) = (id(9), id(10), id(11), id(12));
+    let (e0, e1, along, seam_line, seam) = (id(13), id(14), id(15), id(16), id(17));
+    let (o0, o1, o2, o3, lp, bound, surface, face) = (id(18), id(19), id(20), id(21), id(22), id(23), id(24), id(25));
+    let mut s = net0 + &net1;
+    let xr = x + r;
+    s += &format!("#{p0}=CARTESIAN_POINT('',({xr:?},0.,0.));\n#{p1}=CARTESIAN_POINT('',({xr:?},0.,{h:?}));\n");
+    s += &format!("#{v0}=VERTEX_POINT('',#{p0});\n#{v1}=VERTEX_POINT('',#{p1});\n");
+    s += &format!("#{c0}=CARTESIAN_POINT('',({x:?},0.,0.));\n#{c1}=CARTESIAN_POINT('',({x:?},0.,{h:?}));\n");
+    s += &format!("#{up}=DIRECTION('',(0.,0.,1.));\n#{down}=DIRECTION('',(0.,0.,-1.));\n#{ref_dir}=DIRECTION('',(1.,0.,0.));\n");
+    s += &format!("#{a0}=AXIS2_PLACEMENT_3D('',#{c0},#{up},#{ref_dir});\n#{a1}=AXIS2_PLACEMENT_3D('',#{c1},#{down},#{ref_dir});\n");
+    s += &format!("#{rim0}=CIRCLE('',#{a0},{r:?});\n#{rim1}=CIRCLE('',#{a1},{r:?});\n");
+    s += &format!("#{e0}=EDGE_CURVE('',#{v0},#{v0},#{rim0},.T.);\n#{e1}=EDGE_CURVE('',#{v1},#{v1},#{rim1},.T.);\n");
+    s += &format!("#{along}=VECTOR('',#{up},1.);\n#{seam_line}=LINE('',#{p0},#{along});\n#{seam}=EDGE_CURVE('',#{v0},#{v1},#{seam_line},.T.);\n");
+    s += &format!(
+        "#{o0}=ORIENTED_EDGE('',*,*,#{e0},.F.);\n#{o1}=ORIENTED_EDGE('',*,*,#{seam},.T.);\n#{o2}=ORIENTED_EDGE('',*,*,#{e1},.F.);\n#{o3}=ORIENTED_EDGE('',*,*,#{seam},.F.);\n"
+    );
+    s += &format!("#{lp}=EDGE_LOOP('',(#{o0},#{o1},#{o2},#{o3}));\n#{bound}=FACE_BOUND('',#{lp},.T.);\n");
+    s += &format!(
+        "#{surface}=(BOUNDED_SURFACE()B_SPLINE_SURFACE(1,2,(({}),({})),.UNSPECIFIED.,.F.,.T.,.F.)B_SPLINE_SURFACE_WITH_KNOTS((2,2),({PERIODIC_MULTS}),(0.,{:?}),({PERIODIC_KNOTS}),.UNSPECIFIED.)GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_SURFACE((({CIRCLE_WEIGHTS}),({CIRCLE_WEIGHTS})))REPRESENTATION_ITEM('')SURFACE());\n",
+        row0.join(","),
+        row1.join(","),
+        h / r
+    );
+    s += &format!("#{face}=ADVANCED_FACE('',(#{bound}),#{surface},.T.);\n");
+    (s, face)
+}
+
+/// A B-spline whose end knots repeat fewer than degree + 1 times (an unclamped knot vector,
+/// as OpenCascade writes closed curves and surfaces) is defined only between knot `degree`
+/// and knot `len - degree - 1`. Read over the whole knot range, the circle and the tube went
+/// on past their domain, and the rational ones divided by zero at the outermost knots.
+#[test]
+fn periodic_knot_bsplines_are_read_on_their_domain() {
+    use mt::{BoundedCurve, ParametricCurve, ParametricSurface};
+    let (r, h) = (1.7, 1.5);
+    let (tube, face) = periodic_knot_tube(100, 0.0, r, h);
+    let (circle, row) = circle_net(10, 0.0, r, 0.0);
+    let text = format!(
+        "{HEAD}{tube}{circle}#30=(BOUNDED_CURVE()B_SPLINE_CURVE(2,({}),.CIRCULAR_ARC.,.T.,.F.)B_SPLINE_CURVE_WITH_KNOTS(({PERIODIC_MULTS}),({PERIODIC_KNOTS}),.UNSPECIFIED.)CURVE()GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_CURVE(({CIRCLE_WEIGHTS}))REPRESENTATION_ITEM(''));\n{TAIL}",
+        row.join(",")
+    );
+    let ex = p21::parse(&text).unwrap();
+    let cx = Ctx { ex: &ex, len: 1.0, ang: 1.0, tol: 1e-6 };
+    let pi = std::f64::consts::PI;
+    let on_tube = |p: mt::Point3| (p.x.hypot(p.y) - r).abs() < 1e-9 && p.z > -1e-9 && p.z < h + 1e-9;
+
+    let surface_id = ex.get(face).and_then(|e| e.params().get(2)).and_then(Param::as_ref_id).unwrap();
+    let surface = geom::bspline_surface(&cx, ex.get(surface_id).unwrap()).unwrap();
+    let (Some((u0, u1)), Some((v0, v1))) = surface.try_range_tuple() else { panic!("unbounded surface") };
+    assert!(u0.abs() < 1e-12 && (u1 - h / r).abs() < 1e-12, "u range {u0} .. {u1}");
+    assert!((v0 + pi).abs() < 1e-9 && (v1 - pi).abs() < 1e-9, "v range {v0} .. {v1}");
+    for i in 0..=8 {
+        for j in 0..=16 {
+            let (u, v) = (u0 + (u1 - u0) * i as f64 / 8.0, v0 + (v1 - v0) * j as f64 / 16.0);
+            assert!(on_tube(surface.subs(u, v)), "S({u}, {v}) = {:?}", surface.subs(u, v));
+        }
+    }
+
+    let geom::CurveGeo::Nurbs(curve) = geom::curve(&cx, 30, 0).unwrap() else { panic!("not a rational B-spline") };
+    let (t0, t1) = curve.range_tuple();
+    assert!((t0 + pi).abs() < 1e-9 && (t1 - pi).abs() < 1e-9, "curve range {t0} .. {t1}");
+    for j in 0..=16 {
+        let t = t0 + (t1 - t0) * j as f64 / 16.0;
+        assert!(on_tube(curve.subs(t)), "C({t}) = {:?}", curve.subs(t));
+    }
+}
+
+/// Two such tubes far apart in one body: the body's size sets a chord tolerance coarse beside
+/// the tubes, and meshing a face at it divided the surface's whole parameter box, whose
+/// corners lay outside the domain. The division never converged there (#1374's recursion),
+/// and the import grew by gigabytes without end.
+#[test]
+fn periodic_knot_tube_faces_import_and_mesh() {
+    let t0 = std::time::Instant::now();
+    let (r, h) = (1.7, 1.5);
+    let (a, fa) = periodic_knot_tube(100, 0.0, r, h);
+    let (b, fb) = periodic_knot_tube(200, 700.0, r, h);
+    let text = format!(
+        "{HEAD}{a}{b}#1=OPEN_SHELL('',(#{fa},#{fb}));\n#2=SHELL_BASED_SURFACE_MODEL('Tubes',(#1));\n#3=SHAPE_REPRESENTATION('',(#2),$);\n{TAIL}"
+    );
+    let imp = step_import(&text).unwrap();
+    assert_eq!(imp.bodies.len(), 1, "{:?}", imp.warnings);
+    let m = measure(&imp.bodies[0].body).unwrap();
+    let area = 2.0 * std::f64::consts::TAU * r * h;
+    assert!(rel(m.area, area) < 0.02, "area {} (expected {area})", m.area);
+    assert!(t0.elapsed().as_secs_f64() < 60.0, "took {:?}", t0.elapsed());
+}