diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index e5a76862..dfb516b7 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -106,7 +106,7 @@ in THIRD-PARTY-LICENSES.txt; no images or artwork from that project are copied. | Directory | Source | Licence | Changes | |---|---|---|---| -| vendor/truck-meshalgo | truck-meshalgo 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | meshing: bounded parameter searches, subdivision and provisional fixed-grid previews (vendor/README.md) | +| vendor/truck-meshalgo | truck-meshalgo 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | meshing: bounded parameter searches, collapsed-boundary projection, subdivision and provisional fixed-grid previews (vendor/README.md) | | vendor/truck-stepio | truck-stepio 0.3.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | STEP writer: intersection-curve surfaces written at unique ids (vendor/README.md) | | vendor/truck-shapeops | truck-shapeops 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | robust meshing and nearest-parameter fallback in booleans; degenerate intersection-curve rejection; face and edge history of booleans; trace macro (vendor/README.md) | | vendor/matext4cgmath | matext4cgmath 0.1.0, https://github.com/ricosjp/matext4cgmath (RICOS Co. Ltd.) | Apache-2.0 | root refinement stops after 100 Newton steps (vendor/README.md) | diff --git a/crates/kernel/src/blend_variable.rs b/crates/kernel/src/blend_variable.rs index c18742a8..e0d7a495 100644 --- a/crates/kernel/src/blend_variable.rs +++ b/crates/kernel/src/blend_variable.rs @@ -819,3 +819,7 @@ mod tests { assert!(rel(measure(&f).unwrap().volume, measure(&k).unwrap().volume) < 1e-9); } } + +#[cfg(test)] +#[path = "taper_mesh_tests.rs"] +mod taper_mesh_tests; diff --git a/crates/kernel/src/taper_mesh_tests.rs b/crates/kernel/src/taper_mesh_tests.rs new file mode 100644 index 00000000..25c1acea --- /dev/null +++ b/crates/kernel/src/taper_mesh_tests.rs @@ -0,0 +1,177 @@ +use super::*; +use mt::{ParameterDivision2D, ParametricSurface, ParametricSurface3D, SearchNearestParameter, SearchParameter}; +use truck_meshalgo::prelude::*; + +fn wedge() -> (Body, RadiusSpline) { + use solidcraft_geom::{Loop2, Plane, Region2, Vec2}; + let outer = Loop2::polygon(&[Vec2::new(0., 0.), Vec2::new(40., 0.), Vec2::new(0., 20.)]); + let body = crate::extrude(&Plane::XY, &[Region2 { outer, holes: vec![] }], 0., 15.).unwrap().pop().unwrap(); + let at = Vec3::new(20., 0., 15.); + let (filleted, tapers) = fillet_tapered(&body, &[at], 16.).unwrap(); + let taper = tapers.into_iter().flatten().next().unwrap(); + (filleted, taper.spline) +} + +fn expected(law: &RadiusSpline) -> f64 { + let n = 20_000; + let integral = (0..n).map(|i| law.eval((i as f64 + 0.5) / n as f64).powi(2)).sum::() / n as f64; + 6000. - (1. - std::f64::consts::FRAC_PI_4) * 40. * integral +} + +fn mesh_volume(polygon: &PolygonMesh) -> f64 { + let positions = polygon.positions(); + polygon + .faces() + .triangle_iter() + .map(|t| { + let [a, b, c] = t.map(|v| positions[v.pos].to_vec()); + a.dot(b.cross(c)) / 6. + }) + .sum() +} + +#[test] +fn collapsed_taper_mesh_stays_accurate_as_tolerance_decreases() { + let (body, law) = wedge(); + let wanted = expected(&law); + for tol in [0.002, 0.001, 0.0005, 0.0002] { + for refine in [false, true] { + let mesh = body.tessellate_with(tol, refine).unwrap(); + assert!(mesh.positions.iter().all(|p| p.x.is_finite() && p.y.is_finite() && p.z.is_finite())); + let got = mesh.measure().volume; + assert!((got - wanted).abs() / wanted < 1e-4, "tol={tol}, refine={refine}: {got} vs {wanted}"); + } + } +} + +/// A valid exact parameter search can choose any parameter along a collapsed +/// edge, even far outside its nominal range. Force that answer on every platform +/// instead of relying on the Windows Newton solve that originally exposed #796. +#[derive(Clone)] +struct ArbitraryPole(mt::Surface); +impl ParametricSurface for ArbitraryPole { + type Point = mt::Point3; + type Vector = mt::Vector3; + fn subs(&self, u: f64, v: f64) -> Self::Point { + self.0.subs(u, v) + } + fn uder(&self, u: f64, v: f64) -> Self::Vector { + self.0.uder(u, v) + } + fn vder(&self, u: f64, v: f64) -> Self::Vector { + self.0.vder(u, v) + } + fn uuder(&self, u: f64, v: f64) -> Self::Vector { + self.0.uuder(u, v) + } + fn uvder(&self, u: f64, v: f64) -> Self::Vector { + self.0.uvder(u, v) + } + fn vvder(&self, u: f64, v: f64) -> Self::Vector { + self.0.vvder(u, v) + } + fn parameter_range(&self) -> (mt::ParameterRange, mt::ParameterRange) { + self.0.parameter_range() + } + fn u_period(&self) -> Option { + self.0.u_period() + } + fn v_period(&self) -> Option { + self.0.v_period() + } +} +impl ParametricSurface3D for ArbitraryPole { + fn normal(&self, u: f64, v: f64) -> mt::Vector3 { + self.0.normal(u, v) + } +} +impl ParameterDivision2D for ArbitraryPole { + fn parameter_division(&self, range: ((f64, f64), (f64, f64)), tol: f64) -> (Vec, Vec) { + self.0.parameter_division(range, tol) + } +} +impl SearchParameter for ArbitraryPole { + type Point = mt::Point3; + fn search_parameter>(&self, point: mt::Point3, hint: H, trials: usize) -> Option<(f64, f64)> { + if matches!(self.0, mt::Surface::NurbsSurface(_)) && (point - mt::Point3::new(40., 0., 15.)).magnitude2() < 1e-20 { + Some((1., -1e14)) + } else { + self.0.search_parameter(point, hint, trials) + } + } +} +impl SearchNearestParameter for ArbitraryPole { + type Point = mt::Point3; + fn search_nearest_parameter>(&self, point: mt::Point3, hint: H, trials: usize) -> Option<(f64, f64)> { + self.0.search_nearest_parameter(point, hint, trials) + } +} + +#[test] +fn both_mesh_searches_keep_an_arbitrary_exact_pole_inside_the_domain() { + let (body, law) = wedge(); + let shell = body.solid.boundaries()[0].mapped(Clone::clone, Clone::clone, |s| ArbitraryPole(s.clone())); + for robust in [false, true] { + let mesh = if robust { shell.robust_triangulation(0.0005) } else { shell.triangulation(0.0005) }; + for (source, face) in shell.face_iter().zip(mesh.face_iter()) { + let polygon = face.surface().expect("valid pole must not discard its face"); + if matches!(source.surface().0, mt::Surface::NurbsSurface(_)) { + assert!(polygon.uv_coords().iter().all(|p| (-1e-6..=1. + 1e-6).contains(&p.x) && (-1e-6..=1. + 1e-6).contains(&p.y))); + } + } + let polygon = mesh.to_polygon(); + let volume = mesh_volume(&polygon); + assert!((volume - expected(&law)).abs() / expected(&law) < 1e-4, "robust={robust}: {volume}"); + } +} + +#[test] +fn bounded_projection_does_not_clip_a_planes_unbounded_support() { + let (body, law) = wedge(); + let shell = body.solid.boundaries()[0].mapped(Clone::clone, Clone::clone, |surface| { + if let mt::Surface::Plane(plane) = surface { + let origin = plane.subs(0., 0.); + mt::Surface::Plane(mt::Plane::new(origin, origin + plane.uder(0., 0.).normalize(), origin + plane.vder(0., 0.).normalize())) + } else { + surface.clone() + } + }); + for robust in [false, true] { + let mesh = if robust { shell.robust_triangulation(0.0005) } else { shell.triangulation(0.0005) }; + let mut outside = false; + for (source, face) in shell.face_iter().zip(mesh.face_iter()) { + let polygon = face.surface().unwrap(); + if matches!(source.surface(), mt::Surface::Plane(_)) { + outside |= polygon.uv_coords().iter().any(|p| p.x.abs() > 2. || p.y.abs() > 2.); + } + } + assert!(outside, "fixture must extend beyond the plane's nominal unit box"); + let got = mesh_volume(&mesh.to_polygon()); + assert!((got - expected(&law)).abs() / expected(&law) < 1e-4, "robust={robust}: {got}"); + } +} + +#[test] +fn provisional_preview_bounds_the_same_arbitrary_pole() { + use mt::{BoundedCurve, ParameterDivision1D, ParametricCurve}; + let (body, _) = wedge(); + let face = body.solid.face_iter().find(|f| matches!(f.surface(), mt::Surface::NurbsSurface(_))).unwrap(); + let loops: Vec> = face + .absolute_boundaries() + .iter() + .map(|wire| { + wire.edge_iter() + .flat_map(|edge| { + let curve = edge.oriented_curve(); + let (mut parameters, _) = curve.parameter_division(curve.range_tuple(), 0.01); + parameters.pop(); + parameters.into_iter().map(move |t| curve.subs(t)) + }) + .collect() + }) + .collect(); + let mesh = truck_meshalgo::tessellation::preview_face(&ArbitraryPole(face.surface()), &loops, 8).unwrap(); + assert!(mesh.faces().triangle_iter().next().is_some()); + assert!(mesh.positions().iter().all(|p| p.x.is_finite() && p.y.is_finite() && p.z.is_finite())); + assert!(mesh.uv_coords().iter().all(|p| (-1e-6..=1. + 1e-6).contains(&p.x) && (-1e-6..=1. + 1e-6).contains(&p.y))); +} diff --git a/vendor/README.md b/vendor/README.md index 5ec4577c..88638dff 100644 --- a/vendor/README.md +++ b/vendor/README.md @@ -52,6 +52,13 @@ Changes: It does not run adaptive subdivision or grid balancing. The kernel caller bounds projection complexity and total allocation before using it; ordinary modelling tessellation keeps its existing accuracy policy. +- `tessellation/triangulation.rs`: exact, robust and provisional-preview projections prefer an + equivalent in-domain parameter at collapsed surface edges, and snap to exact bounded + endpoints only at collapsed edges when that improves the spatial fit within the existing + search tolerance. Closed spline seams retain their distinct parameter endpoints. + This prevents huge arbitrary pole parameters from corrupting trimmed meshes (#796), + while retaining genuine extrapolation for unbounded planes. Work is bounded to four + endpoint evaluations per projected point; no search or tessellation tolerance is relaxed. - `tessellation/triangulation.rs`: the inside test of a face's boundary skips loops whose parameter box is clear of the point (it tested every segment of every loop, so a face with hundreds of holes meshed in quadratic time). diff --git a/vendor/truck-meshalgo/src/tessellation/triangulation.rs b/vendor/truck-meshalgo/src/tessellation/triangulation.rs index 0349d9bc..13ccb4ec 100644 --- a/vendor/truck-meshalgo/src/tessellation/triangulation.rs +++ b/vendor/truck-meshalgo/src/tessellation/triangulation.rs @@ -46,17 +46,47 @@ fn clamped(surface: &S, (u, v): (f64, f64)) -> (f64, f64) (c(u, ur, surface.u_period()), c(v, vr, surface.v_period())) } -/// SolidCraft: a nearest-parameter search can run off a B-spline's domain and extrapolate to -/// thousands of parameter units (meshing that grid never finishes). Outside the nominal domain -/// keep the result only if it is really closer to the point than the clamped one (planes report -/// a nominal unit range but extend for ever). -fn sane(surface: &S, point: Point3, uv: (f64, f64)) -> (f64, f64) { - if in_domain(surface, uv) { - return uv; - } +/// SolidCraft: exact and nearest searches can run off a B-spline's domain, especially at a +/// collapsed edge. Prefer a bounded equivalent within the existing spatial search tolerance; +/// otherwise retain extrapolation only when it is closer (planes have a nominal unit range +/// but extend for ever). At most four boundary candidates are evaluated, with no new search. +fn sane(surface: &S, point: Point3, uv: (f64, f64)) -> (f64, f64) { let c = clamped(surface, uv); let d = |(u, v): (f64, f64)| surface.subs(u, v).distance2(point); - if d(uv) <= d(c) { uv } else { c } + let (mut best, mut error) = (uv, d(uv)); + let dc = d(c); + let collapsed = |(u, v): (f64, f64), axis| { + let tangent = if axis == 0 { surface.vder(u, v) } else { surface.uder(u, v) }; + tangent.magnitude2() <= TOLERANCE2 + }; + if !in_domain(surface, uv) + && (dc < error || (dc <= TOLERANCE2 && (collapsed(c, 0) || collapsed(c, 1)))) { + (best, error) = (c, dc); + } + // SolidCraft: at a collapsed boundary an exact search can return arbitrary, + // enormous parameters, or stop just inside the pole. Prefer an equivalent + // in-domain point and the exact boundary when it improves the spatial fit. + // A plane's nominal unit box must not clip genuinely unbounded geometry. + let (ur, vr) = surface.try_range_tuple(); + for (axis, range, period) in [(0, ur, surface.u_period()), (1, vr, surface.v_period())] { + if period.is_some() { continue; } + if let Some((a, b)) = range { + for end in [a, b] { + if !end.is_finite() { continue; } + let coordinate = if axis == 0 { best.0 } else { best.1 }; + let other = if end == a { b } else { a }; + if (coordinate - end).abs() > (coordinate - other).abs() { continue; } + let candidate = if axis == 0 { (end, best.1) } else { (best.0, end) }; + let distance = d(candidate); + // A closed spline seam has equivalent endpoints too, but moving + // between them changes the trim loop. Only canonicalize poles. + if distance <= TOLERANCE2 && distance < error && collapsed(candidate, axis) { + (best, error) = (candidate, distance); + } + } + } + } + best } pub(super) fn by_search_parameter( @@ -83,8 +113,8 @@ where surface .search_parameter(point, hint, 100) .or_else(|| surface.search_parameter(point, None, 100)) - .or_else(|| surface.search_nearest_parameter(point, hint, 100).map(|uv| sane(surface, point, uv))) - .or_else(|| surface.search_nearest_parameter(point, None, 100).map(|uv| sane(surface, point, uv))) + .or_else(|| surface.search_nearest_parameter(point, hint, 100)) + .or_else(|| surface.search_nearest_parameter(point, None, 100)) } /// Tessellates faces @@ -393,7 +423,9 @@ impl PolyBoundaryPiece { .into_iter() .flat_map(|pt| { let (mut u, mut v) = match sp(surface, pt, previous) { - Some(hint) => hint, + // SolidCraft: validate exact and nearest results alike, once per + // boundary point, including the bounded provisional preview. + Some(hint) => sane(surface, pt, hint), None => return vec![None], }; if let (Some(up), Some((u0, _))) = (up, previous) { @@ -704,8 +736,8 @@ pub(super) fn preview_face(surface: &S, loops: &[Vec