Equal cylinders crossing each other no longer refuse a cut (#788) (#1179)

Equal cylinders crossing each other no longer refuse a cut (#788).

crates/kernel/src/boolean_guard.rs

@@ -19,7 +19,7 @@ use crate::{BoolOp, KernelError, Result};
 const MAX_SAMPLES: usize = 20_000;
 
 /// Distance from `p` to the triangle `abc`.
-fn point_triangle_dist(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> f64 {
+pub(crate) fn point_triangle_dist(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> f64 {
     let n = (b - a).cross(c - a);
     let Some(nn) = n.normalized() else {
         return [(a, b), (b, c), (c, a)].iter().map(|(u, v)| p.dist_to_segment(*u, *v)).fold(f64::INFINITY, f64::min);

crates/kernel/src/issue_788_tests.rs

@@ -0,0 +1,150 @@
+//! Two bores of one radius crossing each other (#788): the surfaces touch where the two ellipses
+//! of the crossing meet, which a mesh intersection cannot trace and a boolean cannot take in the
+//! middle of a face.
+use crate::tangent_cylinders;
+use crate::*;
+use solidcraft_geom::Vec3;
+use std::f64::consts::PI;
+
+const SIDE: f64 = 45.0;
+
+fn volume(b: &Body) -> f64 {
+    measure(b).unwrap().volume
+}
+
+fn block() -> Body {
+    let h = SIDE / 2.0;
+    box_solid(Vec3::new(-h, -h, -h), Vec3::new(h, h, h)).unwrap()
+}
+
+/// A cylinder of radius `r` along `axis` through `centre`, longer than the block.
+fn bore(centre: Vec3, axis: Vec3, r: f64, length: f64) -> Body {
+    cylinder(centre - axis * (length / 2.0), axis, r, length).unwrap()
+}
+
+/// The block with a bore along Z, then one of the same radius along `axis` through `centre`.
+fn crossed(r: f64, centre: Vec3, axis: Vec3) -> Result<Option<Body>> {
+    let one = boolean(&block(), &bore(Vec3::ZERO, Vec3::Z, r, 60.0), BoolOp::Cut)?.expect("first bore");
+    boolean(&one, &bore(centre, axis, r, 60.0), BoolOp::Cut)
+}
+
+fn within(got: f64, want: f64, what: &str) {
+    assert!((got - want).abs() < 2e-3 * want, "{what}: {got} vs {want}");
+}
+
+/// Block less two perpendicular bores through its middle: each takes `π r² L`, and the Steinmetz
+/// solid they share, 16r³/3, is counted twice.
+fn expected_perpendicular(r: f64) -> f64 {
+    SIDE.powi(3) - 2.0 * PI * r * r * SIDE + 16.0 * r.powi(3) / 3.0
+}
+
+#[test]
+fn equal_bores_crossing_at_right_angles() {
+    for r in [10.0, 17.5, 4.0] {
+        let body = crossed(r, Vec3::ZERO, Vec3::Y).unwrap_or_else(|e| panic!("r {r}: {e}")).unwrap();
+        assert_eq!(body.shell_count(), 1, "r {r}");
+        within(volume(&body), expected_perpendicular(r), &format!("r {r}"));
+    }
+}
+
+#[test]
+fn equal_bores_crossing_along_every_axis_pair() {
+    let r = 10.0;
+    let want = expected_perpendicular(r);
+    for (first, second) in [(Vec3::X, Vec3::Y), (Vec3::Y, Vec3::X), (Vec3::X, Vec3::Z), (Vec3::Z, Vec3::X), (Vec3::Y, Vec3::Z)] {
+        let one = boolean(&block(), &bore(Vec3::ZERO, first, r, 60.0), BoolOp::Cut).unwrap().unwrap();
+        let body = boolean(&one, &bore(Vec3::ZERO, second, r, 60.0), BoolOp::Cut).unwrap().unwrap();
+        within(volume(&body), want, &format!("{first:?} then {second:?}"));
+    }
+}
+
+#[test]
+fn equal_bores_crossing_at_an_angle() {
+    // Two infinite cylinders of radius r at angle t share 16r³ / (3 sin t).
+    let (r, t) = (5.0f64, PI / 3.0);
+    let axis = Vec3::new(0.0, t.sin(), t.cos());
+    let body = crossed(r, Vec3::ZERO, axis).unwrap().unwrap();
+    let first = PI * r * r * SIDE;
+    let second = PI * r * r * SIDE / t.sin();
+    let shared = 16.0 * r.powi(3) / (3.0 * t.sin());
+    within(volume(&body), SIDE.powi(3) - first - second + shared, "60 degrees");
+}
+
+#[test]
+fn equal_bores_crossing_off_the_middle_of_the_block() {
+    // Both axes meet at (3, -2, 5): the touching points are not at the block's centre.
+    let r = 8.0;
+    let c = Vec3::new(3.0, -2.0, 5.0);
+    let one = boolean(&block(), &bore(c, Vec3::Z, r, 60.0), BoolOp::Cut).unwrap().unwrap();
+    let body = boolean(&one, &bore(c, Vec3::Y, r, 80.0), BoolOp::Cut).unwrap().unwrap();
+    // The bores leave the block through its faces; the shared part lies inside it.
+    let removed = PI * r * r * SIDE * 2.0 - 16.0 * r.powi(3) / 3.0;
+    within(volume(&body), SIDE.powi(3) - removed, "off centre");
+}
+
+#[test]
+fn equal_cylinders_joined_and_intersected() {
+    let r = 10.0;
+    let (z, y) = (bore(Vec3::ZERO, Vec3::Z, r, 60.0), bore(Vec3::ZERO, Vec3::Y, r, 60.0));
+    let tee = boolean(&z, &y, BoolOp::Union).unwrap().unwrap();
+    within(volume(&tee), 2.0 * PI * r * r * 60.0 - 16.0 * r.powi(3) / 3.0, "join");
+    assert_eq!(tee.shell_count(), 1);
+    let steinmetz = boolean(&z, &y, BoolOp::Intersect).unwrap().unwrap();
+    within(volume(&steinmetz), 16.0 * r.powi(3) / 3.0, "intersect");
+}
+
+#[test]
+fn the_result_keeps_a_history_over_the_operands_as_given() {
+    let r = 10.0;
+    let one = boolean(&block(), &bore(Vec3::ZERO, Vec3::Z, r, 60.0), BoolOp::Cut).unwrap().unwrap();
+    let tool = bore(Vec3::ZERO, Vec3::Y, r, 60.0);
+    let body = boolean(&one, &tool, BoolOp::Cut).unwrap().unwrap();
+    let history = body.history().expect("the cut has a history");
+    assert_eq!(history.operands, 2);
+    assert_eq!(history.faces.len(), body.face_count());
+}
+
+#[test]
+fn faces_split_for_the_crossing_leave_the_operands_alone() {
+    // The split is on the working copies: the target keeps its faces.
+    let r = 10.0;
+    let one = boolean(&block(), &bore(Vec3::ZERO, Vec3::Z, r, 60.0), BoolOp::Cut).unwrap().unwrap();
+    let (faces, volume_before) = (one.face_count(), volume(&one));
+    let _ = boolean(&one, &bore(Vec3::ZERO, Vec3::Y, r, 60.0), BoolOp::Cut).unwrap();
+    assert_eq!(one.face_count(), faces);
+    assert_eq!(volume(&one), volume_before);
+}
+
+#[test]
+fn crossings_that_are_not_touching_are_left_alone() {
+    let r = 10.0;
+    let size = 100.0;
+    let tol = size * 1e-6;
+    let cyl = |face, origin: Vec3, axis: Vec3, radius| tangent_cylinders::Cylinder::new(face, origin, axis, radius);
+    let z = cyl(0, Vec3::ZERO, Vec3::Z, r);
+    // Different radius, parallel axes, skew axes: no touching points to split at.
+    assert!(tangent_cylinders::crossing(&z, &cyl(1, Vec3::ZERO, Vec3::Y, r * 0.9), tol).is_none());
+    assert!(tangent_cylinders::crossing(&z, &cyl(1, Vec3::new(5.0, 0.0, 0.0), Vec3::Z, r), tol).is_none());
+    assert!(tangent_cylinders::crossing(&z, &cyl(1, Vec3::new(0.0, 0.0, 0.0), Vec3::new(1e-6, 0.0, 1.0), r), tol).is_none());
+    assert!(tangent_cylinders::crossing(&z, &cyl(1, Vec3::new(1.0, 0.0, 0.0), Vec3::Y, r), tol).is_none());
+    // Axes that meet: two touching points, a radius either side of the crossing.
+    let found = tangent_cylinders::crossing(&z, &cyl(1, Vec3::ZERO, Vec3::Y, r), tol).unwrap();
+    assert!(found.points.iter().any(|p| p.dist(Vec3::new(r, 0.0, 0.0)) < 1e-9));
+    assert!(found.points.iter().any(|p| p.dist(Vec3::new(-r, 0.0, 0.0)) < 1e-9));
+}
+
+#[test]
+fn hostile_cylinders_are_refused() {
+    let tol = 1e-4;
+    let cyl = |axis: Vec3, radius, origin: Vec3| tangent_cylinders::Cylinder::new(0, origin, axis, radius);
+    let good = cyl(Vec3::Z, 10.0, Vec3::ZERO);
+    for bad in [
+        cyl(Vec3::new(f64::NAN, 0.0, 1.0), 10.0, Vec3::ZERO),
+        cyl(Vec3::Y, f64::NAN, Vec3::ZERO),
+        cyl(Vec3::Y, f64::INFINITY, Vec3::ZERO),
+        cyl(Vec3::Y, 10.0, Vec3::new(f64::INFINITY, 0.0, 0.0)),
+        cyl(Vec3::ZERO, 10.0, Vec3::ZERO),
+    ] {
+        assert!(tangent_cylinders::crossing(&good, &bad, tol).is_none(), "{bad:?}");
+    }
+}

crates/kernel/src/lib.rs

@@ -49,6 +49,8 @@ mod issue_674_tests;
 #[cfg(test)]
 mod issue_679_tests;
 #[cfg(test)]
+mod issue_788_tests;
+#[cfg(test)]
 mod issue_851_tests;
 #[cfg(test)]
 mod issue_871_tests;
@@ -89,6 +91,7 @@ mod surface_kind;
 mod surface_offset;
 mod surfaces;
 mod sweep_path;
+mod tangent_cylinders;
 mod tess;
 #[cfg(test)]
 mod thread_tests;

crates/kernel/src/ops.rs

@@ -300,7 +300,21 @@ thread_local! {
 /// The result carries its exact history ([`Body::history`]) over the operands `[a, b]`,
 /// whichever way it was made, or none where a step it took does not report one.
 pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
-    let (mut r, mut diagnostics) = capture_boolean_diagnostics(|| boolean_whole(a, b, op));
+    let r = boolean_via(a, b, op, &boolean_whole);
+    if r.is_ok() || !crate::tangent_cylinders::needed(a, b) {
+        return r;
+    }
+    // Equal cylinders crossing each other touch at two points the boolean cannot take in the
+    // middle of a face (#788): once more with the faces split along them.
+    let split = |a: &Body, b: &Body, op| {
+        crate::tangent_cylinders::boolean(a, b, op, &boolean_whole).unwrap_or_else(|| Err(KernelError::Failed("no split of faces helps".into())))
+    };
+    boolean_via(a, b, op, &split).or(r)
+}
+
+/// [`boolean`], the boolean of the whole bodies taken by `whole`.
+fn boolean_via(a: &Body, b: &Body, op: BoolOp, whole: &dyn Fn(&Body, &Body, BoolOp) -> Result<Option<Body>>) -> Result<Option<Body>> {
+    let (mut r, mut diagnostics) = capture_boolean_diagnostics(|| whole(a, b, op));
     // Bodies in several pieces: when the whole failed (or came out empty), a piece at a time.
     let retry = match &r {
         Err(_) => true,

crates/kernel/src/tangent_cylinders.rs

@@ -0,0 +1,200 @@
+//! Cylinders of one radius whose axes cross (#788: a second bore of the same size drilled across
+//! a bore, a pipe tee of equal pipes).
+//!
+//! Two such cylinders meet in two ellipses that cross each other at two points, where the
+//! surfaces touch along a common normal. A mesh intersection cannot trace the curves through
+//! there, and a boolean with a crossing in the middle of a face has no way to tell which of the
+//! curves' four arms lies in which sector of the face. With the touching points on the edges of
+//! the faces instead, each arm is a curve running from an edge to an edge, which the boolean
+//! takes. So before the boolean, each of the two faces is split along the line of its cylinder
+//! through the touching points (the solids are unchanged), and the result's history is composed
+//! back over the operands as they were given.
+
+use solidcraft_geom::{Plane, Vec3};
+use truck_modeling::{ParametricCurve, SearchNearestParameter};
+
+use crate::body::{Body, from_p3, p3};
+use crate::boolean_record as recorded;
+use crate::history::History;
+use crate::splitface::{SplitTool, split_any};
+use crate::{BoolOp, FaceSurface, Result, face_surface};
+
+/// Most faces split before the boolean (a face crossing several cylinders is split for each).
+const MAX_SPLITS: usize = 16;
+
+/// A cylindrical face and its cylinder.
+#[derive(Clone, Copy, Debug)]
+pub(crate) struct Cylinder {
+    face: usize,
+    origin: Vec3,
+    axis: Vec3,
+    radius: f64,
+}
+
+impl Cylinder {
+    #[cfg(test)]
+    pub(crate) fn new(face: usize, origin: Vec3, axis: Vec3, radius: f64) -> Self {
+        Self { face, origin, axis, radius }
+    }
+}
+
+fn cylinders(body: &Body) -> Vec<Cylinder> {
+    (0..body.face_count())
+        .filter_map(|face| match face_surface(body, face)? {
+            FaceSurface::Cylinder { origin, axis, radius } if origin.is_finite() && axis.is_finite() && radius.is_finite() => {
+                Some(Cylinder { face, origin, axis: axis.normalized()?, radius })
+            }
+            _ => None,
+        })
+        .collect()
+}
+
+/// Where two cylinders of one radius with crossing axes touch, and the planes that hold the
+/// touching points and one cylinder's axis each.
+pub(crate) struct Crossing {
+    pub(crate) points: [Vec3; 2],
+    plane_x: Plane,
+    plane_y: Plane,
+}
+
+pub(crate) fn crossing(x: &Cylinder, y: &Cylinder, tol: f64) -> Option<Crossing> {
+    let finite = [x, y].iter().all(|c| c.origin.is_finite() && c.axis.is_finite() && c.radius.is_finite() && c.radius > 0.0);
+    if !finite || !tol.is_finite() || (x.radius - y.radius).abs() > tol {
+        return None;
+    }
+    // Unit axes, as the cross product below and the planes take them.
+    let (x, y) = (&Cylinder { axis: x.axis.normalized()?, ..*x }, &Cylinder { axis: y.axis.normalized()?, ..*y });
+    let across = x.axis.cross(y.axis);
+    // Axes within a hair of parallel touch along a line, which is another matter.
+    let sine = across.len();
+    if !(sine > 1e-4) {
+        return None;
+    }
+    // The closest points of the two axes, which meet in a point.
+    let (b, w) = (x.axis.dot(y.axis), x.origin - y.origin);
+    let (d, e) = (x.axis.dot(w), y.axis.dot(w));
+    let denominator = 1.0 - b * b;
+    if !(denominator > 1e-12) {
+        return None;
+    }
+    let s = (b * e - d) / denominator;
+    let t = (e - b * d) / denominator;
+    let (on_x, on_y) = (x.origin + x.axis * s, y.origin + y.axis * t);
+    if !(on_x.dist(on_y) <= tol) {
+        return None;
+    }
+    let centre = (on_x + on_y) * 0.5;
+    // Both cylinders reach their common normal `c` at a radius from the crossing.
+    let c = across / sine;
+    let points = [centre + c * x.radius, centre - c * x.radius];
+    Some(Crossing { points, plane_x: Plane::from_normal(centre, x.axis.cross(c))?, plane_y: Plane::from_normal(centre, y.axis.cross(c))? })
+}
+
+/// Whether `p` is on face `face` of the body whose mesh is `mesh`: within `reach` of one of the
+/// face's triangles.
+fn on_face(mesh: &solidcraft_geom::Mesh, face: usize, p: Vec3, reach: f64) -> bool {
+    mesh.triangles
+        .iter()
+        .zip(&mesh.tri_face)
+        .filter(|(_, f)| **f as usize == face)
+        .filter_map(|(t, _)| mesh.tri(t))
+        .any(|[a, b, c]| crate::boolean_guard::point_triangle_dist(p, a, b, c) <= reach)
+}
+
+/// Whether `p` lies on an edge of face `face` (within `margin`).
+fn on_boundary(body: &Body, face: usize, p: Vec3, margin: f64) -> bool {
+    let Some(face) = body.solid.face_iter().nth(face) else { return false };
+    face.absolute_boundaries().iter().any(|wire| {
+        wire.edge_iter().any(|edge| {
+            let curve = edge.curve();
+            curve.search_nearest_parameter(p3(p), None, 100).is_some_and(|t| from_p3(curve.subs(t)).dist(p) <= margin)
+        })
+    })
+}
+
+/// The split the bodies need next: `(whether in the first body, face, plane)`, if a touching
+/// point of two cylinders lies on both their faces but on the edges of neither yet.
+fn next_split(a: &Body, b: &Body, tol: f64) -> Option<(bool, usize, Plane)> {
+    let (in_b, in_a) = (cylinders(b), cylinders(a));
+    if in_a.is_empty() || in_b.is_empty() {
+        return None;
+    }
+    let size = a.size().max(b.size());
+    let (margin, chord) = ((size * 1e-6).max(1e-9), (size * 1e-3).max(1e-6));
+    let mut meshes = None;
+    for x in &in_a {
+        for y in &in_b {
+            let Some(crossing) = crossing(x, y, tol) else { continue };
+            if meshes.is_none() {
+                meshes = Some((a.tessellate(chord).ok()?, b.tessellate(chord).ok()?));
+            }
+            let (ma, mb) = meshes.as_ref()?;
+            for p in crossing.points {
+                if !(on_face(ma, x.face, p, 3.0 * chord) && on_face(mb, y.face, p, 3.0 * chord)) {
+                    continue;
+                }
+                if !on_boundary(a, x.face, p, margin) {
+                    return Some((true, x.face, crossing.plane_x));
+                }
+                if !on_boundary(b, y.face, p, margin) {
+                    return Some((false, y.face, crossing.plane_y));
+                }
+            }
+        }
+    }
+    None
+}
+
+/// Whether the bodies have cylinders of one radius crossing each other that touch inside both
+/// their faces, so that [`boolean`] has faces to split.
+pub(crate) fn needed(a: &Body, b: &Body) -> bool {
+    !a.is_mesh() && !b.is_mesh() && !SPLITTING.with(|c| c.get()) && next_split(a, b, (a.size().max(b.size()) * 1e-7).max(1e-9) * 10.0).is_some()
+}
+
+thread_local! {
+    /// Set while the boolean of the split bodies runs (no second round of splitting).
+    static SPLITTING: std::cell::Cell<bool> = const { std::cell::Cell::new(false) };
+}
+
+/// The boolean `op` of `a` and `b` with their faces split along the touching points of equal
+/// cylinders crossing each other, `run` doing the boolean of the split bodies. `None` when the
+/// bodies have no such cylinders, or the boolean of the split bodies gave no (checked) result:
+/// the caller then takes its other ways.
+pub(crate) fn boolean(a: &Body, b: &Body, op: BoolOp, run: &dyn Fn(&Body, &Body, BoolOp) -> Result<Option<Body>>) -> Option<Result<Option<Body>>> {
+    if a.is_mesh() || b.is_mesh() || SPLITTING.with(|c| c.get()) {
+        return None;
+    }
+    let size = a.size().max(b.size());
+    let tol = (size * 1e-7).max(1e-9) * 10.0;
+    let ((mut x, mut hx), (mut y, mut hy)) = ((a.clone(), Some(History::identity(a))), (b.clone(), Some(History::identity(b))));
+    let mut split = false;
+    for _ in 0..MAX_SPLITS {
+        let Some((first, face, plane)) = next_split(&x, &y, tol) else { break };
+        let (body, history) = if first { (&x, &hx) } else { (&y, &hy) };
+        let (mut done, _) = split_any(body, &[face], &SplitTool::Plane(plane)).ok()?;
+        // The split is no authored boundary: a later feature may still merge what it separates.
+        done.regularization_protected = body.regularization_protected;
+        let through = recorded::through(&done, 0, history.as_ref());
+        if first {
+            (x, hx) = (done, through);
+        } else {
+            (y, hy) = (done, through);
+        }
+        split = true;
+    }
+    if !split {
+        return None;
+    }
+    SPLITTING.with(|c| c.set(true));
+    let result = crate::guard("boolean of split faces", || run(&x, &y, op));
+    SPLITTING.with(|c| c.set(false));
+    let Ok(result) = result else { return None };
+    // Held to what a boolean can do, as a retry's result is (see `boolean_guard`).
+    crate::boolean_guard::check_retry(a, b, op, result.as_ref()).ok()?;
+    let counts = crate::record::counts(&[a, b]);
+    Some(Ok(result.map(|body| {
+        // Each step's history over the operands `[a, b]`.
+        let history = recorded::of_pair(&body, &recorded::placed(hx.clone(), 0, 2), &recorded::placed(hy.clone(), 1, 2));
+        body.with_checked_history(history, &counts)
+    })))
+}