From 172782705a343716cebb89d0dce8571407218090 Mon Sep 17 00:00:00 2001 From: NOWN-DEV <257771410+NOWN-DEV@users.noreply.github.com> Date: Sat, 17 Oct 2026 11:30:13 +0200 Subject: [PATCH] fix: STEP import reads unclamped (periodic) B-splines on the seam (#1916) Thanks @NOWN-DEV. --- crates/kernel/src/step_in/geom.rs | 59 ++++++++++++++- crates/kernel/src/step_in/tests.rs | 117 +++++++++++++++++++++++++++++ 2 files changed, 172 insertions(+), 4 deletions(-) diff --git a/crates/kernel/src/step_in/geom.rs b/crates/kernel/src/step_in/geom.rs index ef2d3a28..91dcb1a1 100644 --- a/crates/kernel/src/step_in/geom.rs +++ b/crates/kernel/src/step_in/geom.rs @@ -216,6 +216,57 @@ fn generated_knots(kind: &str, n_ctrl: usize, degree: usize) -> R { 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, Option)> { + 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 + mt::Tolerance>(mut c: mt::BSplineCurve

) -> mt::BSplineCurve

{ + 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 + mt::Tolerance>(mut s: mt::BSplineSurface

) -> mt::BSplineSurface

{ + 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> { 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 { Some(r) => { let w = weights_list(r.first())?; let n = mt::NurbsCurve::::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 { return Err("weights do not match the control net".into()); } let n = mt::NurbsSurface::::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))), } } diff --git a/crates/kernel/src/step_in/tests.rs b/crates/kernel/src/step_in/tests.rs index 6dd7f4ea..8a195af6 100644 --- a/crates/kernel/src/step_in/tests.rs +++ b/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) { + 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()); +}