From 9be8bb5e6b58832c11b9a7381593ac741efdf020 Mon Sep 17 00:00:00 2001 From: Moritz Date: Sat, 17 Oct 2026 19:33:17 +0200 Subject: [PATCH] Fix Steam Frame STEP trims and spline projection (#1974) (#2043) Thanks @moritzuehling. --- crates/kernel/src/step_in/brep.rs | 74 +++++++++- crates/kernel/src/step_in/tests.rs | 117 +++++++++++++++ .../kernel/tests/trimmed_spline_projection.rs | 92 ++++++++++++ vendor/README.md | 7 + .../truck-meshalgo/src/tessellation/closed.rs | 121 ++++++++++++++++ vendor/truck-meshalgo/src/tessellation/mod.rs | 13 ++ .../src/tessellation/triangulation.rs | 136 +++++++++++++++--- 7 files changed, 536 insertions(+), 24 deletions(-) create mode 100644 crates/kernel/tests/trimmed_spline_projection.rs create mode 100644 vendor/truck-meshalgo/src/tessellation/closed.rs diff --git a/crates/kernel/src/step_in/brep.rs b/crates/kernel/src/step_in/brep.rs index 2ae1dbb1..cea5a334 100644 --- a/crates/kernel/src/step_in/brep.rs +++ b/crates/kernel/src/step_in/brep.rs @@ -4,7 +4,7 @@ use std::collections::HashMap; -use mt::{BoundedCurve, InnerSpace, MetricSpace, ParametricCurve}; +use mt::{BoundedCurve, InnerSpace, MetricSpace, ParametricCurve, ParametricSurface}; use truck_modeling as mt; use super::Ctx; @@ -50,6 +50,62 @@ fn rid(p: Option<&Param>, what: &str) -> R { p.and_then(Param::as_ref_id).ok_or_else(|| format!("{what}: expected a reference")) } +/// A contractible outer loop whose winding opposes the surface's parameter normal. +/// Samples follow the wire through periodic seams; a loop that winds around a surface is +/// not a disk in parameter space and must retain the file's sense. +fn reversed_outer(surface: &mt::Surface, wire: &mt::Wire, tol: f64) -> bool { + if let mt::Surface::RevolutedCurve(proc) = surface { + let rc = proc.entity(); + let curve = rc.entity_curve(); + let on_axis = |p: mt::Point3| (p - rc.origin()).cross(rc.axis()).magnitude() <= tol; + // A compact support has two legitimate sides of a closed trim (a sphere's + // hemispheres, for example); its winding must choose the side, not repair it. + if (on_axis(curve.front()) && on_axis(curve.back())) || curve.front().distance(curve.back()) <= tol { + return false; + } + } + let periods = truck_meshalgo::tessellation::surface_periods(surface); + let mut first = None::<[f64; 2]>; + let mut previous = None::<[f64; 2]>; + let (mut area, mut magnitude) = (0.0, 0.0); + for edge in wire.edge_iter() { + let curve = edge.oriented_curve(); + let (a, b) = curve.range_tuple(); + for k in 0..8 { + let point = curve.subs(a + (b - a) * f64::from(k) / 8.0); + let hint = previous.map(|p| (p[0], p[1])); + let valid = |&(u, v): &(f64, f64)| u.is_finite() && v.is_finite() && surface.subs(u, v).distance(point) <= tol; + let Some((u, v)) = truck_meshalgo::tessellation::project_trim_point(surface, point, hint, tol).filter(valid) else { + return false; + }; + // A pole has no unique parameter in its collapsed direction. Its arbitrary + // search result cannot determine the winding of an outer loop. + let (du, dv) = (surface.uder(u, v), surface.vder(u, v)); + if du.cross(dv).magnitude2() <= 1e-12 * (du.magnitude2() + dv.magnitude2()).powi(2) { + return false; + } + let mut p = [u, v]; + let origin = *first.get_or_insert(p); + if let Some(prev) = previous { + for (d, x) in p.iter_mut().enumerate() { + if let Some(period) = periods.get(d).copied().flatten().filter(|p| *p > 0.0) { + *x -= ((*x - prev[d]) / period).round() * period; + } + } + let cross = (prev[0] - origin[0]) * (p[1] - origin[1]) - (p[0] - origin[0]) * (prev[1] - origin[1]); + area += cross; + magnitude += cross.abs(); + } + previous = Some(p); + } + } + let (Some(first), Some(last)) = (first, previous) else { return false }; + if periods.iter().enumerate().any(|(d, p)| p.is_some_and(|p| (first[d] - last[d]).abs() > p * 0.5)) { + return false; + } + area < -16.0 * f64::EPSILON * magnitude +} + impl Builder { fn vertex(&mut self, cx: &Ctx, id: u64) -> R { if let Some(v) = self.vertices.get(&id) { @@ -205,9 +261,13 @@ impl Builder { return Err(format!("a face with more than {MAX_BOUNDS} boundary loops")); } let mut wires = Vec::new(); + let mut outer = None; let mut edges = 0usize; for b in bounds { if let Some(w) = self.bound(cx, b)? { + if entity(cx, b)?.name() == "FACE_OUTER_BOUND" { + outer = Some(wires.len()); + } edges = edges.saturating_add(w.len()); if edges > MAX_LOOP { return Err(format!("a face with more than {MAX_LOOP} boundary edges")); @@ -243,7 +303,7 @@ impl Builder { (0..=8).map(move |i| c.subs(t0 + (t1 - t0) * i as f64 / 8.0)).collect::>() }) .collect(); - let surface = if e.name() == "FACE" { + let mut surface = if e.name() == "FACE" { crate::plane_fit::plane_through(&extent).map(mt::Surface::Plane).ok_or("polygonal face without a plane")? } else { let sid = rid(p.get(2), "face geometry")?; @@ -251,6 +311,16 @@ impl Builder { guard("step surface", || geom::surface(cx, sid, &extent, &seam, sense, 0).map_err(crate::KernelError::Failed)) .map_err(|e| e.to_string())? }; + // An explicit outer bound encloses the face. Some exporters give its surface the + // opposite sense; treating that loop as a hole meshes the rest of the supporting + // cylinder instead. Match the surface to the outer loop without reversing + // shared topological edges. Winding loops (bands and pole seams) are left alone. + if e.name() != "FACE" && seam.is_empty() && outer.and_then(|i| wires.get(i)).is_some_and(|w| reversed_outer(&surface, w, cx.tol)) { + let sid = rid(p.get(2), "face geometry")?; + let sense = p.get(3).and_then(Param::as_bool).unwrap_or(true); + surface = guard("step surface orientation", || geom::surface(cx, sid, &extent, &seam, !sense, 0).map_err(crate::KernelError::Failed)) + .map_err(|e| e.to_string())?; + } // Loops on a surface of revolution start away from the axis: a pole (cone apex, sphere // pole) has no definite rotation parameter, and the mesher takes the first point's // parameter as the anchor for the rest of the loop. diff --git a/crates/kernel/src/step_in/tests.rs b/crates/kernel/src/step_in/tests.rs index 41d841ba..4f04dea2 100644 --- a/crates/kernel/src/step_in/tests.rs +++ b/crates/kernel/src/step_in/tests.rs @@ -1483,3 +1483,120 @@ fn periodic_knot_tube_faces_import_and_mesh() { assert!(rel(m.area, area) < 0.02, "area {} (expected {area})", m.area); assert!(t0.elapsed().as_secs_f64() < 60.0, "took {:?}", t0.elapsed()); } + +/// An original analytic STEP patch: a 7 mm high, 0.002 radian strip on a 3 m +/// cylinder. Its explicit outer bound encloses the strip even when the surface's +/// supplied sense disagrees with the loop (#1974). +fn shallow_cylinder_patch(sense: bool, outer: bool) -> String { + shallow_cylinder_patch_at_radius(sense, outer, 3000.0) +} + +fn shallow_cylinder_patch_at_radius(sense: bool, outer: bool, r: f64) -> String { + let (angle, h) = (6.0 / r, 7.0); + let (x, y) = (10.0 + r * (angle.cos() - 1.0), r * angle.sin()); + let origin = 10.0 - r; + let bound = if outer { "FACE_OUTER_BOUND" } else { "FACE_BOUND" }; + let sense = if sense { "T" } else { "F" }; + format!( + "{HEAD} +#1=CARTESIAN_POINT('',({origin},0.,0.)); +#2=DIRECTION('',(0.,0.,1.)); +#3=DIRECTION('',(1.,0.,0.)); +#4=AXIS2_PLACEMENT_3D('',#1,#2,#3); +#5=CYLINDRICAL_SURFACE('',#4,{r}); +#6=CARTESIAN_POINT('',(10.,0.,0.)); +#7=CARTESIAN_POINT('',({x},{y},0.)); +#8=CARTESIAN_POINT('',({x},{y},{h})); +#9=CARTESIAN_POINT('',(10.,0.,{h})); +#10=VERTEX_POINT('',#6); +#11=VERTEX_POINT('',#7); +#12=VERTEX_POINT('',#8); +#13=VERTEX_POINT('',#9); +#14=CIRCLE('',#4,{r}); +#15=CARTESIAN_POINT('',({origin},0.,{h})); +#16=AXIS2_PLACEMENT_3D('',#15,#2,#3); +#17=CIRCLE('',#16,{r}); +#18=VECTOR('',#2,1.); +#19=LINE('',#7,#18); +#20=LINE('',#6,#18); +#21=EDGE_CURVE('',#10,#11,#14,.T.); +#22=EDGE_CURVE('',#11,#12,#19,.T.); +#23=EDGE_CURVE('',#12,#13,#17,.F.); +#24=EDGE_CURVE('',#13,#10,#20,.F.); +#25=ORIENTED_EDGE('',*,*,#21,.T.); +#26=ORIENTED_EDGE('',*,*,#22,.T.); +#27=ORIENTED_EDGE('',*,*,#23,.T.); +#28=ORIENTED_EDGE('',*,*,#24,.T.); +#29=EDGE_LOOP('',(#25,#26,#27,#28)); +#30={bound}('',#29,.T.); +#31=ADVANCED_FACE('',(#30),#5,.{sense}.); +#32=OPEN_SHELL('',(#31)); +#33=SHELL_BASED_SURFACE_MODEL('Strip',(#32)); +#34=SHAPE_REPRESENTATION('',(#33),$); +{TAIL}" + ) +} + +#[test] +fn explicit_outer_bound_keeps_shallow_cylinder_patch_trimmed() { + for (sense, radius) in [(true, 3000.), (false, 3000.), (true, 188_980.), (false, 188_980.)] { + let imp = step_import(&shallow_cylinder_patch_at_radius(sense, true, radius)).unwrap(); + assert_eq!(imp.bodies.len(), 1, "{:?}", imp.warnings); + let b = &imp.bodies[0]; + assert_eq!(b.file_faces, 1); + assert!(!b.closed); + for mesh in [ + b.body.tessellate(0.01).unwrap(), + b.body.reference_display_mesh_angular(0.01, crate::DISPLAY_ANGLE).unwrap(), + b.body.display_mesh_angular(0.01, crate::DISPLAY_ANGLE).unwrap(), + ] { + assert!(!mesh.triangles.is_empty()); + assert!(mesh.bounds().diagonal() < 10.0, "sense={sense}: {:?}", mesh.bounds()); + assert!((mesh.measure().area - 42.0).abs() < 0.01, "sense={sense}: area {}", mesh.measure().area); + assert!(mesh.normals.iter().all(|n| n.x > 0.99), "the outer loop determines the face normal"); + } + } +} + +#[test] +fn implicit_bound_can_still_describe_the_complement_of_a_patch() { + let imp = step_import(&shallow_cylinder_patch(false, false)).unwrap(); + let mesh = imp.bodies[0].body.tessellate(0.1).unwrap(); + assert!(mesh.bounds().diagonal() > 5000.0, "without an explicit outer bound the loop is a hole"); +} + +#[test] +fn explicit_outer_bound_crossing_a_closed_nurbs_seam_stays_trimmed() { + // The same analytic strip, with its cylinder expressed as a rational B-spline. + // Put the parameter seam halfway through the strip, independently of its edges. + let mut points = String::new(); + let mut rows = Vec::new(); + let mut weights = Vec::new(); + let w = std::f64::consts::FRAC_1_SQRT_2; + let circle = [(1., 0., 1.), (1., 1., w), (0., 1., 1.), (-1., 1., w), (-1., 0., 1.), (-1., -1., w), (0., -1., 1.), (1., -1., w), (1., 0., 1.)]; + let (sin, cos) = 0.001_f64.sin_cos(); + for (i, (x, y, weight)) in circle.into_iter().enumerate() { + let id = 1000 + i * 2; + let (x, y) = (10.0 - 3000.0 + 3000.0 * (x * cos - y * sin), 3000.0 * (x * sin + y * cos)); + points += &format!("#{id}=CARTESIAN_POINT('',({x:.16},{y:.16},0.));\n#{}=CARTESIAN_POINT('',({x:.16},{y:.16},7.));\n", id + 1); + rows.push(format!("(#{id},#{})", id + 1)); + weights.push(format!("({weight:.16},{weight:.16})")); + } + let surface = format!( + "#5=(BOUNDED_SURFACE() B_SPLINE_SURFACE(2,1,({}),.UNSPECIFIED.,.T.,.F.,.F.) B_SPLINE_SURFACE_WITH_KNOTS((3,2,2,2,3),(2,2),(0.,1.,2.,3.,4.),(0.,1.),.UNSPECIFIED.) GEOMETRIC_REPRESENTATION_ITEM() RATIONAL_B_SPLINE_SURFACE(({})) REPRESENTATION_ITEM('') SURFACE());", + rows.join(","), + weights.join(",") + ); + for sense in [true, false] { + let text = shallow_cylinder_patch(sense, true); + let old = "#5=CYLINDRICAL_SURFACE('',#4,3000);"; + assert!(text.contains(old)); + let text = text.replacen(old, &surface, 1).replace(TAIL, &format!("{points}{TAIL}")); + let imp = step_import(&text).unwrap(); + let body = &imp.bodies[0].body; + let mesh = body.reference_display_mesh_angular(0.01, crate::DISPLAY_ANGLE).unwrap(); + assert!(!mesh.triangles.is_empty(), "{:?}", imp.warnings); + assert!(mesh.bounds().diagonal() < 10.0, "sense={sense}: {:?}", mesh.bounds()); + assert!((mesh.measure().area - 42.0).abs() < 0.01, "sense={sense}: area {}", mesh.measure().area); + } +} diff --git a/crates/kernel/tests/trimmed_spline_projection.rs b/crates/kernel/tests/trimmed_spline_projection.rs new file mode 100644 index 00000000..26a63090 --- /dev/null +++ b/crates/kernel/tests/trimmed_spline_projection.rs @@ -0,0 +1,92 @@ +//! Original analytic surfaces exercising trim projection (#1974), without an external CAD fixture. +use truck_meshalgo::prelude::*; +use truck_modeling::{self as mt, BoundedCurve, InnerSpace}; + +fn mesh(surface: mt::Surface, curves: Vec) -> PolygonMesh { + let vertices: Vec<_> = curves.iter().map(|c| mt::Vertex::new(c.front())).collect(); + let count = vertices.len(); + let edges: Vec<_> = curves.into_iter().enumerate().map(|(i, c)| mt::Edge::new_unchecked(&vertices[i], &vertices[(i + 1) % count], c)).collect(); + let face = mt::Face::new_unchecked(vec![edges.into()], surface); + let shell = mt::Shell::from(vec![face]).robust_triangulation(0.001); + let polygon = shell.face_iter().next().unwrap().surface().unwrap(); + assert!(polygon.faces().triangle_iter().next().is_some(), "the trimmed face must remain visible"); + polygon +} + +#[test] +fn approximate_boundary_does_not_escape_a_bounded_spline() { + // S(u,v) = (50u + 10v, .001v, 10u² + 4v), 0 <= u,v <= 1. + // Near u=1 its tangents are almost parallel. A slightly approximate top edge + // sends an undamped search hinted from the preceding point to u>1. + let p = mt::Point3::new; + let surface = mt::BSplineSurface::new( + (mt::KnotVec::bezier_knot(2), mt::KnotVec::bezier_knot(1)), + vec![vec![p(0., 0., 0.), p(10., 0.001, 4.)], vec![p(25., 0., 0.), p(35., 0.001, 4.)], vec![p(50., 0., 10.), p(60., 0.001, 14.)]], + ); + let curves = vec![ + mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), vec![p(60., 0.001, 14.), p(35., 0., 4.), p(10., 0.001, 4.)]).into(), + mt::Line(p(10., 0.001, 4.), p(0., 0., 0.)).into(), + mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), vec![p(0., 0., 0.), p(25., 0., 0.), p(50., 0., 10.)]).into(), + mt::Line(p(50., 0., 10.), p(60., 0.001, 14.)).into(), + ]; + let mesh = mesh(mt::Surface::BSplineSurface(surface), curves); + for uv in mesh.uv_coords() { + assert!((-1e-6..=1. + 1e-6).contains(&uv.x) && (-1e-6..=1. + 1e-6).contains(&uv.y), "parameter {uv:?}"); + } + assert!(mesh.positions().iter().all(|p| p.x.is_finite() && p.y.is_finite() && p.z.is_finite())); +} + +#[test] +fn a_trim_crossing_a_closed_nurbs_seam_is_a_small_cylinder_patch() { + let (radius, angle, bottom, top) = (10.0_f64, 0.1_f64, 2.0, 5.0); + let w = std::f64::consts::FRAC_1_SQRT_2; + let circle = [(1., 0., 1.), (1., 1., w), (0., 1., 1.), (-1., 1., w), (-1., 0., 1.), (-1., -1., w), (0., -1., 1.), (1., -1., w), (1., 0., 1.)]; + let ctrl = circle + .iter() + .map(|&(x, y, w)| vec![mt::Vector4::new(radius * x * w, radius * y * w, 0., w), mt::Vector4::new(radius * x * w, radius * y * w, 10. * w, w)]) + .collect(); + let surface = mt::NurbsSurface::new(mt::BSplineSurface::new( + (mt::KnotVec::from(vec![0., 0., 0., 1., 1., 2., 2., 3., 3., 4., 4., 4.]), mt::KnotVec::bezier_knot(1)), + ctrl, + )); + let at = |theta: f64, z| mt::Point3::new(radius * theta.cos(), radius * theta.sin(), z); + let arc = |z, reverse| { + let w = angle.cos(); + let mut ctrl = vec![at(-angle, z).to_homogeneous(), mt::Vector4::new(radius, 0., z * w, w), at(angle, z).to_homogeneous()]; + if reverse { + ctrl.reverse(); + } + mt::NurbsCurve::new(mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), ctrl)).into() + }; + let curves = vec![ + arc(bottom, false), + mt::Line(at(angle, bottom), at(angle, top)).into(), + arc(top, true), + mt::Line(at(-angle, top), at(-angle, bottom)).into(), + ]; + let mesh = mesh(mt::Surface::NurbsSurface(surface), curves); + for p in mesh.positions() { + assert!(p.x >= radius * angle.cos() - 1e-6 && p.x <= radius + 1e-6, "{p:?}"); + assert!(p.y.abs() <= radius * angle.sin() + 1e-6 && (bottom - 1e-6..=top + 1e-6).contains(&p.z), "{p:?}"); + } + let area: f64 = mesh + .faces() + .triangle_iter() + .map(|t| { + let [a, b, c] = [0, 1, 2].map(|i| mesh.positions()[t[i].pos]); + (b - a).cross(c - a).magnitude() * 0.5 + }) + .sum(); + assert!((area - radius * 2.0 * angle * (top - bottom)).abs() < 0.01, "area {area}"); +} + +#[test] +fn an_affine_planes_unit_chart_does_not_clip_a_trim() { + let p = mt::Point3::new; + let points = [p(10., 20., 0.), p(13., 20., 0.), p(13., 22., 0.), p(10., 22., 0.)]; + let curves = (0..4).map(|i| mt::Line(points[i], points[(i + 1) % 4]).into()).collect(); + let mesh = mesh(mt::Plane::new(p(0., 0., 0.), p(1., 0., 0.), p(0., 1., 0.)).into(), curves); + for p in mesh.positions() { + assert!((10. - 1e-6..=13. + 1e-6).contains(&p.x) && (20. - 1e-6..=22. + 1e-6).contains(&p.y), "{p:?}"); + } +} diff --git a/vendor/README.md b/vendor/README.md index f1770570..41521026 100644 --- a/vendor/README.md +++ b/vendor/README.md @@ -83,6 +83,13 @@ Changes: ## truck-meshalgo 0.4.0 (Apache-2.0, © RICOS Co. Ltd., https://github.com/ricosjp/truck) Changes: +- `tessellation/triangulation.rs`, `closed.rs`: retry nearest trim projections outside a + spline's knot domain without the preceding point's hint, preferring a bounded fit within + the caller's chord tolerance, with a bounded damped least-squares fallback where + Newton's method diverges. Affine plane charts remain unbounded; modelling surfaces may + retain genuine continuation when no bounded fit is accurate enough. Geometrically closed, + smooth spline seams get a local periodic chart during meshing, with wrapped evaluation, + so a trim crossing the seam stays a small patch. The stored geometry is unchanged (#1974). - `tessellation/mod.rs`, `triangulation.rs`: a separate provisional `preview_face` path accepts at most 4,096 already sampled boundary points and a fixed grid of at most 16 subdivisions per direction, including periodic seam connections. diff --git a/vendor/truck-meshalgo/src/tessellation/closed.rs b/vendor/truck-meshalgo/src/tessellation/closed.rs new file mode 100644 index 00000000..ee8891c1 --- /dev/null +++ b/vendor/truck-meshalgo/src/tessellation/closed.rs @@ -0,0 +1,121 @@ +//! SolidCraft: a local periodic chart for geometrically closed spline surfaces. +//! The stored spline still has a finite knot domain; only meshing wraps its parameters. +use super::*; + +#[derive(Clone, Debug)] +pub(super) struct Closed<'a, S> { + pub base: &'a S, + ranges: [Option<(f64, f64)>; 2], + periods: [Option; 2], + inferred: [bool; 2], +} + +impl<'a, S: PreMeshableSurface> Closed<'a, S> { + pub fn new(base: &'a S) -> Self { + let (u, v) = base.try_range_tuple(); + let mut out = Self { base, ranges: [u, v], periods: [base.u_period(), base.v_period()], inferred: [false; 2] }; + let (Some((u0, u1)), Some((v0, v1))) = (u, v) else { return out }; + if ![u0, u1, v0, v1].iter().all(|v| v.is_finite()) || u1 <= u0 || v1 <= v0 { + return out; + } + let origin = base.subs(u0, v0); + let scale = (0..=4) + .flat_map(|i| (0..=4).map(move |j| (i, j))) + .map(|(i, j)| base.subs(u0 + (u1 - u0) * i as f64 / 4., v0 + (v1 - v0) * j as f64 / 4.).distance(origin)) + .fold(0.0_f64, f64::max); + let tol2 = (scale * 1e-8).max(1e-9).powi(2); + for d in 0..2 { + if out.periods[d].is_some() { + continue; + } + let (lo, hi, a, b) = if d == 0 { (u0, u1, v0, v1) } else { (v0, v1, u0, u1) }; + let mut noncollapsed = false; + let closed = (0..=8).all(|k| { + let t = a + (b - a) * k as f64 / 8.; + let at = |x| if d == 0 { (x, t) } else { (t, x) }; + let (p, q, m) = (at(lo), at(hi), at((lo + hi) * 0.5)); + let first = base.subs(p.0, p.1); + let (second, middle) = (base.subs(q.0, q.1), base.subs(m.0, m.1)); + if ![first.x, first.y, first.z, second.x, second.y, second.z, middle.x, middle.y, middle.z].iter().all(|x| x.is_finite()) { + return false; + } + noncollapsed |= first.distance2(middle) > tol2; + if first.distance2(second) > tol2 { + return false; + } + let (n0, n1) = (base.normal(p.0, p.1), base.normal(q.0, q.1)); + !n0.magnitude2().is_finite() || !n1.magnitude2().is_finite() || n0.dot(n1) > 1.0 - 1e-6 + }); + if closed && noncollapsed { + out.periods[d] = Some(hi - lo); + out.inferred[d] = true; + } + } + out + } + + pub fn periods(&self) -> [Option; 2] { + self.periods + } + + pub fn wrap(&self, (u, v): (f64, f64)) -> (f64, f64) { + let mut p = [u, v]; + for (d, x) in p.iter_mut().enumerate() { + if let (true, Some((lo, hi)), Some(period)) = (self.inferred[d], self.ranges[d], self.periods[d]) { + if *x < lo || *x > hi { + *x = lo + (*x - lo).rem_euclid(period); + } + } + } + (p[0], p[1]) + } +} + +impl ParametricSurface for Closed<'_, S> { + type Point = Point3; + type Vector = Vector3; + fn subs(&self, u: f64, v: f64) -> Point3 { + let (u, v) = self.wrap((u, v)); + self.base.subs(u, v) + } + fn uder(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.uder(u, v) + } + fn vder(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.vder(u, v) + } + fn uuder(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.uuder(u, v) + } + fn uvder(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.uvder(u, v) + } + fn vvder(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.vvder(u, v) + } + fn parameter_range(&self) -> (ParameterRange, ParameterRange) { + self.base.parameter_range() + } + fn u_period(&self) -> Option { + self.periods[0] + } + fn v_period(&self) -> Option { + self.periods[1] + } +} +impl ParametricSurface3D for Closed<'_, S> { + fn normal(&self, u: f64, v: f64) -> Vector3 { + let (u, v) = self.wrap((u, v)); + self.base.normal(u, v) + } +} +impl ParameterDivision2D for Closed<'_, S> { + fn parameter_division(&self, range: ((f64, f64), (f64, f64)), tol: f64) -> (Vec, Vec) { + triangulation::bounded_surface_division(self, range, Deflection::chord(tol), triangulation::MAX_GRID) + } +} diff --git a/vendor/truck-meshalgo/src/tessellation/mod.rs b/vendor/truck-meshalgo/src/tessellation/mod.rs index 896de00f..6c60fdb3 100644 --- a/vendor/truck-meshalgo/src/tessellation/mod.rs +++ b/vendor/truck-meshalgo/src/tessellation/mod.rs @@ -25,6 +25,18 @@ pub fn preview_face(surface: &S, loops: &[Vec] triangulation::preview_face(surface, loops, steps) } +/// Periods of a surface, including smooth geometric closure of bounded splines. +/// Meshing wraps inferred periods into the original knot domain before evaluation. +pub fn surface_periods(surface: &S) -> [Option; 2] { + closed::Closed::new(surface).periods() +} + +/// Project a trim point with the same bounded, tolerance-aware search used by the mesher. +/// A preceding parameter is a hint, never proof that the search converged correctly. +pub fn project_trim_point(surface: &S, point: Point3, hint: Option<(f64, f64)>, tolerance: f64) -> Option<(f64, f64)> { + triangulation::nearest_projection(surface, point, hint, tolerance, 100) +} + /// SolidCraft: an edge's polyline as the meshers divide it (at most 20 000 points). Exposed for /// the kernel's regressions in the workspace CI. #[doc(hidden)] @@ -508,6 +520,7 @@ impl RobustMeshableShape } } +mod closed; mod triangulation; #[cfg(feature = "test-hooks")] #[doc(hidden)] diff --git a/vendor/truck-meshalgo/src/tessellation/triangulation.rs b/vendor/truck-meshalgo/src/tessellation/triangulation.rs index 61b62402..cc09e76a 100644 --- a/vendor/truck-meshalgo/src/tessellation/triangulation.rs +++ b/vendor/truck-meshalgo/src/tessellation/triangulation.rs @@ -16,9 +16,9 @@ type MeshedShell = Shell>; type MeshedCShell = CompressedShell>; pub(super) trait SP: - Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> + Parallelizable { + Fn(&S, Point3, Option<(f64, f64)>, f64) -> Option<(f64, f64)> + Parallelizable { } -impl SP for F where F: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> + Parallelizable {} +impl SP for F where F: Fn(&S, Point3, Option<(f64, f64)>, f64) -> Option<(f64, f64)> + Parallelizable {} /// SolidCraft: a parameter inside the surface's bounded, non-periodic ranges. fn in_domain(surface: &S, (u, v): (f64, f64)) -> bool { @@ -93,6 +93,7 @@ pub(super) fn by_search_parameter( surface: &S, point: Point3, hint: Option<(f64, f64)>, + _accept: f64, ) -> Option<(f64, f64)> where S: MeshableSurface, @@ -106,15 +107,97 @@ pub(super) fn by_search_nearest_parameter( surface: &S, point: Point3, hint: Option<(f64, f64)>, + accept: f64, ) -> Option<(f64, f64)> -where - S: RobustMeshableSurface, -{ - surface - .search_parameter(point, hint, 100) - .or_else(|| surface.search_parameter(point, None, 100)) - .or_else(|| surface.search_nearest_parameter(point, hint, 100)) - .or_else(|| surface.search_nearest_parameter(point, None, 100)) +where S: RobustMeshableSurface { + nearest_projection(surface, point, hint, accept, 100) +} + +/// A hint can converge to a spline's continuation outside its knot domain. Prefer a +/// bounded fit within the caller's chord tolerance. Some modelling surfaces really are +/// extended beyond their stored domain; retain that continuation when no bounded fit is +/// accurate enough, rather than losing a valid blend face. +pub(super) fn nearest_projection(surface: &S, point: Point3, hint: Option<(f64,f64)>, accept: f64, trials: usize) -> Option<(f64,f64)> { + // An affine plane's unit parameter box is only a chart, not a finite support. + // Preserve its exact projection even when the trim extends beyond that box. + if affine_chart(surface) { + return surface.search_parameter(point, hint, trials) + .or_else(|| surface.search_parameter(point, None, trials)) + .or_else(|| surface.search_nearest_parameter(point, hint, trials)) + .or_else(|| surface.search_nearest_parameter(point, None, trials)); + } + let candidate = |uv| { + let uv = sane(surface, point, uv); + let error = surface.subs(uv.0, uv.1).distance2(point); + (uv.0.is_finite() && uv.1.is_finite() && error.is_finite()).then_some((uv, error)) + }; + let fits = |&(uv, error): &((f64, f64), f64)| in_domain(surface, uv) && error <= accept * accept; + let exact = surface.search_parameter(point, hint, trials) + .or_else(|| surface.search_parameter(point, None, trials)).and_then(candidate); + if let Some((uv, _)) = exact.filter(fits) { return Some(uv); } + let hinted = surface.search_nearest_parameter(point, hint, trials).and_then(candidate); + if let Some((uv, _)) = hinted.filter(fits) { return Some(uv); } + let fresh = surface.search_nearest_parameter(point, None, trials).and_then(candidate); + if let Some((uv, _)) = fresh.filter(fits) { return Some(uv); } + let best = [exact, hinted, fresh].into_iter().flatten().min_by(|a, b| a.1.total_cmp(&b.1)); + if let Some(uv) = best.map(|x| x.0).or(hint) { + if let Some(bounded) = bounded_projection(surface, point, uv).filter(fits) { + return Some(bounded.0); + } + } + // Even the unhinted search can fail for shallow patches of very large cylinders. + // Their domain midpoint supplies a finite start without trusting a previous failure. + if let (Some((a, b)), Some((c, d))) = surface.try_range_tuple() { + if [a, b, c, d].iter().all(|x| x.is_finite()) { + let start = (a * 0.5 + b * 0.5, c * 0.5 + d * 0.5); + if let Some(bounded) = bounded_projection(surface, point, start).filter(fits) { + return Some(bounded.0); + } + } + } + best.map(|x| x.0) +} + +fn affine_chart(surface: &S) -> bool { + let (Some((a, b)), Some((c, d))) = surface.try_range_tuple() else { return false }; + [a, b, c, d].iter().all(|x| x.is_finite()) && [0.0, 0.5, 1.0].into_iter().all(|t| { + let (u, v) = (a + (b - a) * t, c + (d - c) * t); + [surface.uuder(u, v), surface.uvder(u, v), surface.vvder(u, v)] + .iter().all(|x| x.magnitude2() <= 1e-24) + }) +} + +/// Damped least squares constrained to the finite knot domain. An unconstrained Newton +/// step is unreliable where the two tangents almost coincide. This fallback has fixed work +/// limits and is used only after the usual searches failed to find a bounded fit. +fn bounded_projection(surface: &S, point: Point3, start: (f64, f64)) -> Option<((f64, f64), f64)> { + let mut uv = clamped(surface, start); + let mut error = surface.subs(uv.0, uv.1).distance2(point); + let mut lambda = 1e-6_f64; + for _ in 0..60 { + let q = surface.subs(uv.0, uv.1); + let (du, dv) = (surface.uder(uv.0, uv.1), surface.vder(uv.0, uv.1)); + let r = point - q; + let (a, b, c) = (du.magnitude2(), du.dot(dv), dv.magnitude2()); + let scale = a + c; + if !scale.is_finite() || scale <= 0.0 { break; } + let (aa, cc) = (a + lambda * scale, c + lambda * scale); + let det = aa * cc - b * b; + if !det.is_finite() || det <= 0.0 { break; } + let (g, h) = (du.dot(r), dv.dot(r)); + let next = clamped(surface, (uv.0 + (cc * g - b * h) / det, uv.1 + (aa * h - b * g) / det)); + let nd = surface.subs(next.0, next.1).distance2(point); + if nd.is_finite() && nd <= error { + let moved = (next.0 - uv.0).abs() + (next.1 - uv.1).abs(); + (uv, error) = (next, nd); + lambda = (lambda * 0.1).max(1e-12); + if moved <= 1e-15 * (1.0 + uv.0.abs() + uv.1.abs()) { break; } + } else { + lambda *= 10.0; + if lambda > 1e8 { break; } + } + } + (error.is_finite() && in_domain(surface, uv)).then_some((uv, error)) } /// Tessellates faces @@ -341,14 +424,18 @@ where let edges: Vec<_> = shell.edges.iter().map(tessellate_edge).collect(); let tessellate_face = |face: &CompressedFace| { let boundaries = face.boundaries.clone(); - let surface = &face.surface; + let wrapped = closed::Closed::new(&face.surface); + let search = |s: &closed::Closed<'_, S>, point, hint: Option<(f64, f64)>, accept| { + sp(s.base, point, hint.map(|p| s.wrap(p)), accept) + }; + let surface = &wrapped; let create_edge = |edge_idx: &CompressedEdgeIndex| match edge_idx.orientation { true => Some(edges.get(edge_idx.index)?.curve.clone()), false => Some(edges.get(edge_idx.index)?.curve.inverse()), }; let create_boundary = |wire: &Vec| { let wire_iter = wire.iter().filter_map(create_edge); - PolyBoundaryPiece::try_new(surface, wire_iter, &sp) + PolyBoundaryPiece::try_new(surface, wire_iter, &search, tol.chord) }; let preboundary: Option> = boundaries.iter().map(create_boundary).collect(); let polygon: Option = (|| { @@ -380,11 +467,16 @@ fn shell_create_polygon( sp: impl SP, cap: usize, ) -> Face> { + let wrapped = closed::Closed::new(surface); + let search = |s: &closed::Closed<'_, S>, point, hint: Option<(f64, f64)>, accept| { + sp(s.base, point, hint.map(|p| s.wrap(p)), accept) + }; + let surface = &wrapped; let preboundary = wires .iter() .map(|wire: &Wire<_, _>| { let wire_iter = wire.iter().map(Edge::oriented_curve); - PolyBoundaryPiece::try_new(surface, wire_iter, &sp) + PolyBoundaryPiece::try_new(surface, wire_iter, &search, tol.chord) }) .collect::>>(); let polygon: Option = (|| { @@ -451,6 +543,7 @@ impl PolyBoundaryPiece { surface: &S, wire: impl Iterator, sp: impl SP, + accept: f64, ) -> Option { let (up, vp) = (surface.u_period(), surface.v_period()); let (urange, vrange) = surface.try_range_tuple(); @@ -466,7 +559,7 @@ impl PolyBoundaryPiece { bdry3d .into_iter() .flat_map(|pt| { - let (mut u, mut v) = match sp(surface, pt, previous) { + let (mut u, mut v) = match sp(surface, pt, previous, accept) { // SolidCraft: validate exact and nearest results alike, once per // boundary point, including the bounded provisional preview. Some(hint) => sane(surface, pt, hint), @@ -502,7 +595,7 @@ impl PolyBoundaryPiece { let mut vec = match walk(bdry3d.clone()) { Some(vec) => vec, None => { - let start = bdry3d.iter().skip(1).position(|pt| sp(surface, *pt, None).is_some())? + 1; + let start = bdry3d.iter().skip(1).position(|pt| sp(surface, *pt, None, accept).is_some())? + 1; let mut rotated = bdry3d; rotated.rotate_left(start); walk(rotated)? @@ -845,14 +938,13 @@ pub(super) fn preview_face(surface: &S, loops: &[Vec

() > 4096 { return None; } - let search = |surface: &S, point, hint| { - surface.search_parameter(point, hint, 16) - .or_else(|| surface.search_parameter(point, None, 16)) - .or_else(|| surface.search_nearest_parameter(point, hint, 16)) - .or_else(|| surface.search_nearest_parameter(point, None, 16)) + let wrapped = closed::Closed::new(surface); + let search = |s: &closed::Closed<'_, S>, point, hint: Option<(f64, f64)>, accept| { + nearest_projection(s.base, point, hint.map(|p| s.wrap(p)), accept, 16) }; + let surface = &wrapped; let pieces = loops.iter().map(|w| { - PolyBoundaryPiece::try_new(surface, std::iter::once(PolylineCurve(w.clone())), search) + PolyBoundaryPiece::try_new(surface, std::iter::once(PolylineCurve(w.clone())), search, 1e-3) }).collect::>>()?; let boundary = PolyBoundary::with_steps(pieces, surface, Deflection::chord(1e-3), Some(steps)); // Projection can duplicate points at singularities; periodic seams add at most @@ -995,7 +1087,7 @@ pub(super) fn bounded_polyline(curve: &C, range: (f64, f64 /// `cap` points or after 24 rounds; cells that do not evaluate (NaN) are not refined. A /// direction within the chord tolerance is still halved where a cell turns more than the /// deflection's angle along it, while the grid stays within `Deflection::MAX_ANGLE_GRID`. -fn bounded_surface_division(surface: &S, (ur, vr): ((f64, f64), (f64, f64)), d: Deflection, cap: usize) -> (Vec, Vec) { +pub(super) fn bounded_surface_division(surface: &S, (ur, vr): ((f64, f64), (f64, f64)), d: Deflection, cap: usize) -> (Vec, Vec) { let (tol, turn2) = (d.chord, d.turn2()); let (mut udiv, mut vdiv) = (vec![ur.0, ur.1], vec![vr.0, vr.1]); for _ in 0..24 {