diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index f3b0fcee..c49a473f 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -259,6 +259,7 @@ mod path_scope; mod primitive_points; mod robustness; mod robustness_goldens; +mod spiral_join; mod xz_plane; /// Plan area a round of radius `r` fills at the concave corner `x` where circles (c1, r1) and diff --git a/crates/engine/src/tests/spiral_join.rs b/crates/engine/src/tests/spiral_join.rs new file mode 100644 index 00000000..515a17fb --- /dev/null +++ b/crates/engine/src/tests/spiral_join.rs @@ -0,0 +1,44 @@ +//! Two interleaved spiral walls, joined (the case a of #844: 22 mm pitch, 0.75 turn). + +use std::f64::consts::TAU; + +use super::*; + +/// The join used to be an error: the walls' flat tops and bottoms are B-spline faces, so nothing +/// saw that they coincide. Now they are prisms, and the join is decided in 2D. +#[test] +fn two_interleaved_spiral_walls_join_into_one_body() { + let (pitch, turns, points) = (22.0, 0.75, 19usize); + let (thickness, height, start) = (12.0, 60.0, 60.0); + let (k, end) = (pitch / TAU, turns * TAU); + let arm = |sign: f64| -> Vec<[f64; 2]> { + (0..points) + .map(|i| { + let a = end * i as f64 / (points - 1) as f64; + let r = start + k * a; + [sign * r * a.cos(), sign * r * a.sin()] + }) + .collect() + }; + let mut s = Session::default(); + for (tag, sign) in [("A", 1.0), ("B", -1.0)] { + let x0 = sign * start; + run(&mut s, "sketch.create", json!({"plane": "XY", "name": format!("Path{tag}")})); + run(&mut s, "sketch.spline.fit_point", json!({"points": arm(sign)})); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "sketch.create", json!({"plane": "XZ", "name": format!("Prof{tag}")})); + run(&mut s, "sketch.rectangle.two_point", json!({"p0": [x0 - thickness / 2.0, -height], "p1": [x0 + thickness / 2.0, 0]})); + run(&mut s, "sketch.finish", json!({})); + let sweep = json!({ + "sketch": format!("Prof{tag}"), "path_sketch": format!("Path{tag}"), "path": ["s1"], + "operation": "new", "name": format!("Sweep{tag}"), "body_name": format!("Band{tag}"), + }); + run(&mut s, "solid.sweep", sweep); + } + run(&mut s, "solid.combine", json!({"target": "BandA", "tools": ["BandB"], "operation": "join", "name": "Joined"})); + assert_eq!(s.model.state().bodies.len(), 1); + // The two walls' plan areas united, by a plan-view raster (0.6 mm cells, flat ends), times the height. + let want = 450_878.0; + let got = volume(&mut s); + assert!(rel(got, want) < 5e-3, "{got} vs {want}"); +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 83d5bc0e..0df3d1aa 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -59,6 +59,9 @@ mod par; pub use par::{for_each_on_threads, map_on_threads}; #[cfg(test)] mod pinch_tests; +mod planar_faces; +#[cfg(test)] +mod planar_faces_tests; mod plane_fit; mod polybool; mod polyhedron; diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index ce69aa4e..eb6016c4 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -432,7 +432,26 @@ fn prism_result(a: &Body, b: &Body, op: BoolOp, size: f64, counts: &[(usize, usi }) }) } +/// Walls swept along splines have flat B-spline tops and bottoms. With those made planes the +/// operands are prisms, and two of them (same height, one direction) are decided in 2D: fast, and +/// with no intersection of spline surfaces, which is what fails on them. The outlines are fitted +/// by arcs within a micron (`prism_boolean_free_form`). `None` when the operands have no such +/// face, are not such prisms, or the answer is implausible: the general boolean then takes over. +fn swept_prism(a: &Body, b: &Body, op: BoolOp) -> Option { + let (pa, pb) = (crate::planar_faces::planarized(a), crate::planar_faces::planarized(b)); + if pa.is_none() && pb.is_none() { + return None; + } + let (ra, rb) = (pa.as_ref().unwrap_or(a), pb.as_ref().unwrap_or(b)); + let Some(Ok(Some(body))) = crate::prism::prism_boolean_free_form(ra, rb, op) else { return None }; + let (va, vb) = (volume(a), volume(b)); + plausible(volume(&body), op, va, vb, 2e-3 * (va + vb) + 1e-9).then_some(body) +} + fn boolean_whole(a: &Body, b: &Body, op: BoolOp) -> Result> { + if let Some(body) = swept_prism(a, b, op) { + return Ok(Some(body)); + } a.require_brep("a boolean")?; b.require_brep("a boolean")?; let size = a.size().max(b.size()); diff --git a/crates/kernel/src/planar_faces.rs b/crates/kernel/src/planar_faces.rs new file mode 100644 index 00000000..d4df8fa6 --- /dev/null +++ b/crates/kernel/src/planar_faces.rs @@ -0,0 +1,106 @@ +//! Flat faces that sit on B-spline surfaces become planes. +//! +//! A wall swept along a curve has a top and a bottom that are flat but are made of the same +//! B-spline surfaces as everything else the sweep makes. The coplanar handling of a boolean +//! (coincident faces pushed apart, faces glued back to back, the exact planar fallback) looks +//! for faces on true planes, so it never sees them: two such walls that share a top plane were +//! left to the general intersection, which fails on faces that coincide. + +use solidcraft_geom::Vec3; +use truck_modeling as mt; + +use crate::body::{Body, Solid, from_p3, p3}; +use crate::{KernelError, guard}; + +/// Most control points looked at in one surface. +const MAX_CONTROL_POINTS: usize = 1_000_000; + +/// The control points of a B-spline or NURBS surface, as points. +fn control_points(s: &mt::Surface) -> Option> { + match s { + mt::Surface::BSplineSurface(b) => { + let pts: Vec = b.control_points().iter().flatten().take(MAX_CONTROL_POINTS + 1).map(|p| from_p3(*p)).collect(); + (pts.len() <= MAX_CONTROL_POINTS).then_some(pts) + } + mt::Surface::NurbsSurface(n) => { + let mut pts = Vec::new(); + for v in n.control_points().iter().flatten() { + if !(v.w > 0.0) || pts.len() >= MAX_CONTROL_POINTS { + return None; + } + pts.push(Vec3::new(v.x / v.w, v.y / v.w, v.z / v.w)); + } + Some(pts) + } + _ => None, + } +} + +/// A point and the unit normal of the plane every point lies within `tol` of, if there is one. +fn common_plane(pts: &[Vec3], tol: f64) -> Option<(Vec3, Vec3)> { + let p0 = *pts.first()?; + let far = |from: Vec3, pick: &dyn Fn(Vec3) -> f64| pts.iter().copied().max_by(|a, b| pick(*a - from).total_cmp(&pick(*b - from))); + let p1 = far(p0, &|d| d.len())?; + let dir = (p1 - p0).normalized()?; + let p2 = far(p0, &|d| (d - dir * d.dot(dir)).len())?; + let n = (p1 - p0).cross(p2 - p0).normalized()?; + pts.iter().all(|p| n.dot(*p - p0).abs() <= tol).then_some((p0, n)) +} + +/// `b` with every face that is flat on a B-spline surface put on a plane; `None` when there is +/// no such face. The faces keep their order, their edges and their orientation. +pub(crate) fn planarized(b: &Body) -> Option { + if b.is_mesh() { + return None; + } + let tol = (b.size() * 1e-9).max(1e-9); + guard("planar faces", || { + use mt::{ParametricSurface, ParametricSurface3D}; + let mut changed = 0usize; + let mut shells = Vec::new(); + for sh in b.solid.boundaries() { + let mut faces: Vec = Vec::new(); + for f in sh.face_iter() { + let surface = f.surface(); + let plane = control_points(&surface).and_then(|pts| common_plane(&pts, tol)); + let Some((origin, mut n)) = plane else { + faces.push(f.clone()); + continue; + }; + // The plane faces the way the surface does; a surface with no normal at any of the + // places looked at is left as it is. + let (ur, vr) = surface.try_range_tuple(); + let ((u0, u1), (v0, v1)) = (ur.unwrap_or((0.0, 1.0)), vr.unwrap_or((0.0, 1.0))); + let old = [(0.5, 0.5), (0.25, 0.25), (0.75, 0.75), (0.25, 0.75), (0.75, 0.25)] + .iter() + .map(|(s, t)| surface.normal(u0 + s * (u1 - u0), v0 + t * (v1 - v0))) + .find(|v| v.x.is_finite() && v.y.is_finite() && v.z.is_finite() && v.x * v.x + v.y * v.y + v.z * v.z > 1e-18); + let Some(old) = old else { + faces.push(f.clone()); + continue; + }; + if n.dot(Vec3::new(old.x, old.y, old.z)) < 0.0 { + n = -n; + } + let u = n.any_perp().normalized().ok_or_else(|| KernelError::Failed("planar faces: no axis".into()))?; + let w = n.cross(u); + let plane = mt::Plane::new(p3(origin), p3(origin + u), p3(origin + w)); + let mut nf = + mt::Face::try_new(f.absolute_boundaries().to_vec(), mt::Surface::Plane(plane)).map_err(|e| KernelError::Failed(e.to_string()))?; + if !f.orientation() { + nf.invert(); + } + faces.push(nf); + changed += 1; + } + shells.push(faces.into()); + } + if changed == 0 { + return Ok(None); + } + let solid = Solid::try_new(shells).map_err(|e| KernelError::Failed(e.to_string()))?; + Ok(Some(Body::new(solid)?)) + }) + .ok() + .flatten() +} diff --git a/crates/kernel/src/planar_faces_tests.rs b/crates/kernel/src/planar_faces_tests.rs new file mode 100644 index 00000000..1c6788f4 --- /dev/null +++ b/crates/kernel/src/planar_faces_tests.rs @@ -0,0 +1,173 @@ +#![allow(clippy::unwrap_used, clippy::panic)] +//! A wall swept along a spiral has a flat top and bottom that are B-spline faces. They have to be +//! planes for the coplanar handling of a boolean to see them: a join of two such walls ran into +//! the general intersection with faces that coincide, and failed. With planes, the walls are prisms +//! and a join of two is decided in 2D (`prism.rs`). + +use std::f64::consts::{FRAC_PI_2, TAU}; + +use solidcraft_geom::{Loop2, Plane, Region2, Vec2, Vec3}; + +use crate::*; + +const WALL: f64 = 12.0; +const HEIGHT: f64 = 60.0; +/// How far the spiral moves out in one turn. +const PITCH: f64 = 22.0; + +/// The centreline of the spiral that starts `start` radians round from +X at `radius`, `t` radians along. +fn centre(radius: f64, start: f64, t: f64) -> Vec3 { + let (r, a) = (radius + PITCH / TAU * t, start + t); + Vec3::new(r * a.cos(), r * a.sin(), 0.0) +} + +/// A wall `WALL` thick and `HEIGHT` tall along that spiral for `turn` radians, in `stations` steps, +/// the profile square to the path. +fn spiral_wall(radius: f64, start: f64, turn: f64, stations: usize) -> Body { + let path: Vec = (0..=stations).map(|i| centre(radius, start, turn * i as f64 / stations as f64)).collect(); + let tangent = (path[1] - path[0]).normalized().unwrap(); + // x up the wall, y across it: x cross y runs along the path. + let plane = Plane { origin: path[0], x: Vec3::Z, y: tangent.cross(Vec3::Z) }; + let profile = Region2 { + outer: Loop2::polygon(&[ + Vec2::new(0.0, -WALL / 2.0), + Vec2::new(HEIGHT, -WALL / 2.0), + Vec2::new(HEIGHT, WALL / 2.0), + Vec2::new(0.0, WALL / 2.0), + ]), + holes: vec![], + }; + crate::sweep_path(&plane, &profile, &path).unwrap() +} + +/// Two such walls whose centrelines are 11 mm apart, so that they overlap by 1 mm; the second +/// starts a little further round, so that no end face lies in the other's. +fn two_walls(turn: f64, stations: usize) -> (Body, Body) { + (spiral_wall(60.0, 0.0, turn, stations), spiral_wall(71.0, 0.26, turn, stations)) +} + +/// Faces that lie on a `Plane` surface (the face-type counts of `measure` also count flat splines). +fn plane_faces(b: &Body) -> usize { + b.solid.face_iter().filter(|f| matches!(f.surface(), truck_modeling::Surface::Plane(_))).count() +} + +fn volume(b: &Body) -> f64 { + b.tessellate_with(1e-3, true).unwrap().measure().volume +} + +#[test] +fn the_flat_faces_of_a_swept_wall_become_planes() { + let wall = spiral_wall(60.0, 0.0, FRAC_PI_2, 8); + // Two end caps are planes already; the top and bottom are flat B-splines like the two sides. + assert_eq!(plane_faces(&wall), 2); + let flat = crate::planar_faces::planarized(&wall).expect("the top and the bottom are flat"); + assert_eq!(plane_faces(&flat), 4); + assert_eq!(flat.face_count(), wall.face_count()); + assert!(flat.validity().is_empty(), "{:?}", flat.validity()); + let (before, after) = (volume(&wall), volume(&flat)); + assert!((before - after).abs() < 1e-3 * before, "{before} vs {after}"); + // Nothing flat is left to change, and a box has nothing to begin with. + assert!(crate::planar_faces::planarized(&flat).is_none()); + assert!(crate::planar_faces::planarized(&box_solid(Vec3::ZERO, Vec3::new(1.0, 2.0, 3.0)).unwrap()).is_none()); +} + +#[test] +fn walls_that_share_a_top_and_a_bottom_coincide_only_once_those_are_planes() { + let (a, b) = two_walls(FRAC_PI_2, 8); + assert!(crate::coplanar::coincident_planes(&a, &b).is_empty()); + let (pa, pb) = (crate::planar_faces::planarized(&a).unwrap(), crate::planar_faces::planarized(&b).unwrap()); + let planes = crate::coplanar::coincident_planes(&pa, &pb); + assert_eq!(planes.len(), 2, "{planes:?}"); + let mut heights: Vec = planes.iter().map(|(n, d, _)| d * n.z).collect(); + heights.sort_by(f64::total_cmp); + assert!(heights[0].abs() < 1e-6 && (heights[1] - HEIGHT).abs() < 1e-6, "{heights:?}"); +} + +/// The plan area of the union of the walls (each given as `(radius, start)`), by sampling every +/// 0.4 mm: a point is in a wall if it is within half its thickness of the wall's centreline and +/// not beyond either end (flat ends). +fn union_plan_area(walls: &[(f64, f64)], turn: f64) -> f64 { + let lines: Vec> = walls.iter().map(|(r, s)| (0..=400).map(|i| centre(*r, *s, turn * i as f64 / 400.0)).collect()).collect(); + let inside = |line: &[Vec3], p: Vec3| -> bool { + let mut best = (f64::INFINITY, 0.0, 0); + for (j, w) in line.windows(2).enumerate() { + let d = w[1] - w[0]; + let s = ((p - w[0]).dot(d) / d.dot(d)).clamp(0.0, 1.0); + let dist = (p - (w[0] + d * s)).len(); + if dist < best.0 { + best = (dist, (p - w[0]).dot(d) / d.dot(d), j); + } + } + let beyond = (best.2 == 0 && best.1 < 0.0) || (best.2 + 2 == line.len() && best.1 > 1.0); + best.0 <= WALL / 2.0 && !beyond + }; + let h = 0.4; + let n = (2.0 * 100.0 / h) as usize; + let mut cells = 0usize; + for i in 0..n { + for j in 0..n { + let p = Vec3::new(-100.0 + (i as f64 + 0.5) * h, -100.0 + (j as f64 + 0.5) * h, 0.0); + if lines.iter().any(|l| inside(l, p)) { + cells += 1; + } + } + } + cells as f64 * h * h +} + +/// The faces of `b` that are neither planes nor cylinders (the face types `measure` reports). +fn free_form_faces(b: &Body) -> usize { + crate::measure(b).unwrap().merged.face_types.get("other").copied().unwrap_or(0) +} + +fn join_checked(a: &Body, b: &Body, walls: &[(f64, f64)], turn: f64) { + let joined = boolean(a, b, BoolOp::Union).unwrap().unwrap(); + assert_eq!(joined.solid.boundaries().len(), 1); + assert!(joined.validity().is_empty(), "{:?}", joined.validity()); + // Decided in 2D: the walls come back as arcs, not as splines. + assert_eq!(free_form_faces(&joined), 0); + let want = HEIGHT * union_plan_area(walls, turn); + let got = joined.tessellate_with(1e-2, true).unwrap().measure().volume; + assert!((got - want).abs() < 5e-3 * want, "{got} vs {want}"); +} + +/// The general intersection fails on these two walls (`no result`: their tops and bottoms coincide, +/// and as B-splines the coincident-face handling does not see them). With flat faces made planes +/// they are prisms, and the join is decided in 2D, with the walls fitted by arcs within a micron. +#[test] +fn a_join_of_swept_walls_is_decided_in_2d() { + let turn = 0.6; + let (a, b) = two_walls(turn, 16); + join_checked(&a, &b, &[(60.0, 0.0), (71.0, 0.26)], turn); +} + +/// Two interleaved spirals, 0.75 turn, as in the snail: the second is the first turned half a turn. +#[test] +fn interleaved_spiral_walls_join() { + let turn = 0.75 * TAU; + let half = std::f64::consts::PI; + let (a, b) = (spiral_wall(60.0, 0.0, turn, 36), spiral_wall(60.0, half, turn, 36)); + join_checked(&a, &b, &[(60.0, 0.0), (60.0, half)], turn); +} + +/// An extrusion of a spline outline has plane caps already: a cut of it by a cylinder keeps the +/// exact spline wall, whatever the 2D route does for swept walls. +#[test] +fn a_cut_of_an_extruded_spline_keeps_its_spline_wall() { + let p = Vec2::new; + let line = |a: Vec2, b: Vec2| solidcraft_geom::Seg2::Line { a, b }; + let outline = Loop2 { + segs: vec![ + line(p(0.0, 0.0), p(40.0, 0.0)), + line(p(40.0, 0.0), p(40.0, 20.0)), + solidcraft_geom::Seg2::Cubic { p0: p(40.0, 20.0), p1: p(30.0, 35.0), p2: p(10.0, 5.0), p3: p(0.0, 20.0) }, + line(p(0.0, 20.0), p(0.0, 0.0)), + ], + }; + let plate = crate::extrude(&Plane::XY, &[Region2 { outer: outline, holes: vec![] }], 0.0, 10.0).unwrap().remove(0); + assert_eq!(free_form_faces(&plate), 1); + let hole = crate::cylinder(Vec3::new(20.0, 10.0, -1.0), Vec3::Z, 3.0, 12.0).unwrap(); + let cut = boolean(&plate, &hole, BoolOp::Cut).unwrap().unwrap(); + assert_eq!(free_form_faces(&cut), 1); + assert!(cut.face_count() < 12, "{}", cut.face_count()); +} diff --git a/crates/kernel/src/prism.rs b/crates/kernel/src/prism.rs index 213879e1..3fc4c6eb 100644 --- a/crates/kernel/src/prism.rs +++ b/crates/kernel/src/prism.rs @@ -485,7 +485,11 @@ impl Index { for w in f.boundaries() { let (mut segs, mut tags) = (Vec::new(), Vec::new()); for e in w.edge_iter() { - let pieces = crate::prism_curve::segment(e, plane, tol, work)?; + let pieces = if FREE_FORM.with(|c| c.get()) { + crate::prism_curve::segment_fitted(e, plane, tol, PLANAR_FIT_TOL, work)? + } else { + crate::prism_curve::segment(e, plane, tol, work)? + }; *segments = segments.checked_sub(pieces.len())?; let tag = self.tag(e, cap); tags.extend((0..pieces.len()).map(|k| tag.piece(k, pieces.len()))); @@ -1919,6 +1923,28 @@ fn stack_indexed(a: &Body, b: &Body, op: BoolOp, [ia, ib]: [&LazyIndex; 2]) -> O sew_sections(&plane, d, &levels, &slabs, tol, coalesce, &history).map(|r| r.map(Some)) } +thread_local! { + /// Set while [`prism_boolean_free_form`] runs: spline outlines are read as arcs within + /// [`PLANAR_FIT_TOL`] instead of refusing the prism. + static FREE_FORM: std::cell::Cell = const { std::cell::Cell::new(false) }; +} + +/// [`prism_boolean`] for prisms whose walls are splines (swept along one): the outlines are +/// fitted by arcs within [`PLANAR_FIT_TOL`], so the result is exact to that, and the walls it +/// keeps are those arcs. Only for bodies the caller knows are such sweeps. +pub(crate) fn prism_boolean_free_form(a: &Body, b: &Body, op: BoolOp) -> Option>> { + /// Clears the flag however the boolean ends, a caught panic included. + struct Reset; + impl Drop for Reset { + fn drop(&mut self) { + FREE_FORM.with(|c| c.set(false)); + } + } + FREE_FORM.with(|c| c.set(true)); + let _reset = Reset; + prism_boolean(a, b, op) +} + /// Boolean of prisms along a common direction, in 2D (see the module notes); `None` when it /// doesn't apply. The result carries its history over `[a, b]` (see [`lineage`]). pub(crate) fn prism_boolean(a: &Body, b: &Body, op: BoolOp) -> Option>> { diff --git a/crates/kernel/src/prism_curve.rs b/crates/kernel/src/prism_curve.rs index 6e303a1c..cada5e37 100644 --- a/crates/kernel/src/prism_curve.rs +++ b/crates/kernel/src/prism_curve.rs @@ -198,6 +198,46 @@ fn circle_span(c: &[H], center: Vec2, radius: f64, tol: f64, source_error: f64, (a.signum() == b.signum()).then_some(a + b) } +/// [`segment`], except that a spline that is neither a line nor a circular arc comes back as lines +/// and arcs that stay within `fit` of it (the walls of a body swept along a spline). `None` where +/// [`segment`] refuses for another reason: not in the plane, not finite, over the work budget, +/// or a rational spline. +pub(super) fn segment_fitted(e: &mt::Edge, plane: &Plane, tol: f64, fit: f64, work: &mut usize) -> Option> { + if let Some(exact) = segment(e, plane, tol, work) { + return Some(exact); + } + if !tol.is_finite() || tol <= 0. || !fit.is_finite() || fit <= 0. { + return None; + } + let mut c = e.curve(); + preflight(&c, work)?; + if !e.orientation() { + c.invert(); + } + let (t0, t1) = c.range_tuple(); + if !(t0.is_finite() && t1.is_finite() && t1 > t0) { + return None; + } + let (a3, b3) = (from_p3(c.subs(t0)), from_p3(c.subs(t1))); + if !a3.is_finite() || !b3.is_finite() || a3.dist(from_p3(e.front().point())) > tol || b3.dist(from_p3(e.back().point())) > tol { + return None; + } + let (spans, _) = controls(&c, plane, plane.height(a3), tol)?; + let mut out = Vec::new(); + for span in spans { + if span.iter().any(|p| (p[2] - 1.).abs() > 1e-9) { + return None; + } + charge(work, 64)?; + let points: Vec = span.iter().map(|p| point(*p)).collect(); + out.extend(solidcraft_geom::bezier(&points)?.fit_arcs(fit)?); + if out.len() > MAX_CONTROLS { + return None; + } + } + Some(out) +} + pub(super) fn segment(e: &mt::Edge, plane: &Plane, tol: f64, work: &mut usize) -> Option> { if !tol.is_finite() || tol <= 0. { return None;