Count the contact grid's own work, and test find_close_pair (#836)

find_close_pair's budget now counts the grid's own work too (each cell a triangle is filed under or visits), not only the pair tests, so a pathological pair of meshes can't spend unbounded time building and walking the grid. Adds tests for find_close_pair and the budget.

crates/engine/src/tests.rs

@@ -195,6 +195,63 @@ fn join_meeting_nothing_keeps_the_first_body() {
     assert!(rel(volume(&mut s), 1064.0) < 1e-6);
 }
 
+/// Two spheres of one size: the join's contact test tried every pair of their finely meshed
+/// triangles and refused the join as too much work (#232).
+#[test]
+fn join_of_two_overlapping_equal_spheres_is_not_refused() {
+    let mut s = Session::default();
+    run(&mut s, "solid.sphere", json!({"radius": 10, "center": [0, 0, 0], "body_name": "A"}));
+    run(&mut s, "solid.sphere", json!({"radius": 10, "center": [12, 0, 0], "body_name": "B"}));
+    run(&mut s, "solid.combine", json!({"target": "A", "tools": ["B"], "operation": "join"}));
+    assert_eq!(s.model.state().bodies.len(), 1);
+    // Two spheres less the lens they share: 2 * 4/3 pi r^3 - pi (4r + d) (2r - d)^2 / 12.
+    let (r, d) = (10.0, 12.0);
+    let want = 2.0 * 4.0 / 3.0 * PI * r * r * r - PI * (4.0 * r + d) * (2.0 * r - d) * (2.0 * r - d) / 12.0;
+    let got = volume(&mut s);
+    assert!(rel(got, want) < 1e-3, "{got} vs {want}");
+}
+
+/// Area two discs of radii `r1`, `r2` with centres `d` apart have in common.
+fn disc_overlap(r1: f64, r2: f64, d: f64) -> f64 {
+    if d >= r1 + r2 {
+        return 0.0;
+    }
+    if d <= (r1 - r2).abs() {
+        return PI * r1.min(r2) * r1.min(r2);
+    }
+    let a = ((d * d + r1 * r1 - r2 * r2) / (2.0 * d * r1)).acos();
+    let b = ((d * d + r2 * r2 - r1 * r1) / (2.0 * d * r2)).acos();
+    r1 * r1 * a + r2 * r2 * b - 0.5 * ((-d + r1 + r2) * (d + r1 - r2) * (d - r1 + r2) * (d + r1 + r2)).sqrt()
+}
+
+/// Two thick rings 11 mm apart: their walls cross, and the join was refused as too much work (#232).
+#[test]
+fn join_of_two_overlapping_rings_is_not_refused() {
+    let mut s = Session::default();
+    for (name, x) in [("RingA", 0.0), ("RingB", 11.0)] {
+        run(
+            &mut s,
+            "solid.cylinder",
+            json!({"diameter": 120, "height": 60, "base": [x, 0, 0], "operation": "new", "name": format!("{name}Outer"), "body_name": name}),
+        );
+        run(
+            &mut s,
+            "solid.cylinder",
+            json!({"diameter": 96, "height": 60, "base": [x, 0, 0], "operation": "cut", "targets": [name], "name": format!("{name}Bore")}),
+        );
+    }
+    run(&mut s, "solid.combine", json!({"target": "RingA", "tools": ["RingB"], "operation": "join"}));
+    assert_eq!(s.model.state().bodies.len(), 1);
+    // Union of the two rings' plan areas: both less what they share, which is the overlap of the
+    // two outer discs less the two bores (each inside the other's outer disc) plus the bores' own overlap.
+    let (outer, bore, d) = (60.0, 48.0, 11.0);
+    let ring = PI * (outer * outer - bore * bore);
+    let shared = disc_overlap(outer, outer, d) - 2.0 * disc_overlap(bore, outer, d) + disc_overlap(bore, bore, d);
+    let want = (2.0 * ring - shared) * 60.0;
+    let got = volume(&mut s);
+    assert!(rel(got, want) < 1e-3, "{got} vs {want}");
+}
+
 mod boolean_diagnostics;
 mod component_frames;
 mod mirror_planes;

crates/geom/src/mesh.rs

@@ -90,14 +90,32 @@ const MAX_CELLS_PER_TRIANGLE: usize = 64;
 /// `gap` of each other (each such pair once), until it returns true, and says whether it did.
 /// The cost follows the number of close pairs, not `a.len() * b.len()`: `b` is filed in a sparse
 /// grid of cells about its typical triangle's size, and each triangle of `a` meets only those in
-/// its cells. Every call of `test` takes one from `budget`; `None` when the budget runs out.
+/// its cells. The work counts against `budget`, one unit each: a grid cell a triangle of `b` is
+/// filed under, a cell a triangle of `a` visits, and a call of `test`. `None` when the budget runs out.
 pub fn find_close_pair(
     a: &[BoxedTriangle],
     b: &[BoxedTriangle],
     gap: f64,
     budget: &mut usize,
+    test: impl FnMut(&BoxedTriangle, &BoxedTriangle) -> bool,
+) -> Option<bool> {
+    let work = std::cell::Cell::new(*budget);
+    let found = find_close_pair_within(a, b, gap, &work, test);
+    *budget = work.get();
+    found
+}
+
+fn find_close_pair_within(
+    a: &[BoxedTriangle],
+    b: &[BoxedTriangle],
+    gap: f64,
+    work: &std::cell::Cell<usize>,
     mut test: impl FnMut(&BoxedTriangle, &BoxedTriangle) -> bool,
 ) -> Option<bool> {
+    let spend = |n: usize| -> Option<()> {
+        work.set(work.get().checked_sub(n)?);
+        Some(())
+    };
     let gap = if gap.is_finite() { gap.max(0.0) } else { 0.0 };
     let close = |x: &Aabb3, y: &Aabb3| {
         x.min.x <= y.max.x + gap
@@ -108,7 +126,7 @@ pub fn find_close_pair(
             && y.min.z <= x.max.z + gap
     };
     let mut call = |t: &BoxedTriangle, u: &BoxedTriangle| -> Option<bool> {
-        *budget = budget.checked_sub(1)?;
+        spend(1)?;
         Some(test(t, u))
     };
     let (Some(first), true) = (b.first(), !a.is_empty()) else { return Some(false) };
@@ -146,6 +164,7 @@ pub fn find_close_pair(
                 for ix in x.0..=x.1 {
                     for iy in y.0..=y.1 {
                         for iz in z.0..=z.1 {
+                            spend(1)?;
                             cells.entry((ix, iy, iz)).or_default().push(j as u32);
                         }
                     }
@@ -177,6 +196,7 @@ pub fn find_close_pair(
         for ix in x.0..=x.1 {
             for iy in y.0..=y.1 {
                 for iz in z.0..=z.1 {
+                    spend(1)?;
                     for &j in cells.get(&(ix, iy, iz)).into_iter().flatten() {
                         let (Some(mark), Some(u)) = (seen.get_mut(j as usize), b.get(j as usize)) else { continue };
                         if *mark == i + 1 {

crates/geom/src/tests.rs

@@ -319,3 +319,148 @@ fn meshes_cross_only_where_their_surfaces_pass_through_each_other() {
     assert!(!a.crosses(&cube(9.995, 3.0, 1), 0.01));
     assert!(!a.crosses(&cube(9.85, 3.0, 1), f64::NAN));
 }
+
+mod close_pairs {
+    use super::*;
+
+    /// A small deterministic generator, so the random cases below repeat from run to run.
+    struct Lcg(u64);
+
+    impl Lcg {
+        fn unit(&mut self) -> f64 {
+            self.0 = self.0.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1_442_695_040_888_963_407);
+            (self.0 >> 11) as f64 / (1u64 << 53) as f64
+        }
+
+        fn point(&mut self) -> Vec3 {
+            Vec3::new(self.unit(), self.unit(), self.unit())
+        }
+    }
+
+    fn boxed(tri: [Vec3; 3]) -> BoxedTriangle {
+        let mut bounds = Aabb3::EMPTY;
+        for p in tri {
+            bounds.add(p);
+        }
+        (tri, bounds)
+    }
+
+    /// `count` small triangles scattered over a cube of side `spread` at `origin`. On the lattice their
+    /// corners sit on whole numbers, so many touch exactly, share edges or lie in one plane.
+    fn scatter(rng: &mut Lcg, count: usize, origin: Vec3, spread: f64, size: f64, lattice: bool) -> Vec<BoxedTriangle> {
+        let mut out = Vec::with_capacity(count);
+        for _ in 0..count {
+            let centre = origin + rng.point() * spread;
+            let mut tri = [Vec3::ZERO; 3];
+            for corner in &mut tri {
+                let mut p = centre + (rng.point() - Vec3::new(0.5, 0.5, 0.5)) * (2.0 * size);
+                if lattice {
+                    p = Vec3::new(p.x.round(), p.y.round(), p.z.round());
+                }
+                *corner = p;
+            }
+            out.push(boxed(tri));
+        }
+        out
+    }
+
+    /// A pair that "meets" for these tests: two corners within `gap` of each other.
+    fn corners_meet(t: &[Vec3; 3], u: &[Vec3; 3], gap: f64) -> bool {
+        t.iter().any(|p| u.iter().any(|q| p.dist(*q) <= gap))
+    }
+
+    fn boxes_close(x: &Aabb3, y: &Aabb3, gap: f64) -> bool {
+        x.min.x <= y.max.x + gap
+            && y.min.x <= x.max.x + gap
+            && x.min.y <= y.max.y + gap
+            && y.min.y <= x.max.y + gap
+            && x.min.z <= y.max.z + gap
+            && y.min.z <= x.max.z + gap
+    }
+
+    /// Every pair tried, which is how the contact tests worked before they used a grid.
+    fn every_pair(a: &[BoxedTriangle], b: &[BoxedTriangle], gap: f64) -> bool {
+        a.iter().any(|(x, xb)| b.iter().any(|(y, yb)| boxes_close(xb, yb, gap) && corners_meet(x, y, gap)))
+    }
+
+    /// A flat grid of `n` x `n` squares, two triangles each, at height `z`.
+    fn flat_grid(n: usize, z: f64) -> Vec<BoxedTriangle> {
+        let p = |i: usize, j: usize| Vec3::new(i as f64, j as f64, z);
+        let mut out = Vec::with_capacity(2 * n * n);
+        for i in 0..n {
+            for j in 0..n {
+                out.push(boxed([p(i, j), p(i + 1, j), p(i + 1, j + 1)]));
+                out.push(boxed([p(i, j), p(i + 1, j + 1), p(i, j + 1)]));
+            }
+        }
+        out
+    }
+
+    /// The grid hands `test` the pairs that trying every pair finds: on random triangles, on
+    /// triangles that touch exactly (whole-number corners) and far from the origin.
+    #[test]
+    fn the_grid_finds_the_pairs_that_trying_every_pair_finds() {
+        let mut rng = Lcg(0x232);
+        let (mut met, mut apart) = (0, 0);
+        for case in 0..600usize {
+            let lattice = case % 3 == 0;
+            let origin = if case % 5 == 0 { Vec3::new(1.0e6, -2.0e6, 3.0e6) } else { Vec3::ZERO };
+            let gap = [1e-9, 0.01, 0.2][(case / 3) % 3];
+            let (na, nb) = (1 + (rng.unit() * 150.0) as usize, 1 + (rng.unit() * 150.0) as usize);
+            let shift = (rng.point() - Vec3::new(0.5, 0.5, 0.5)) * 8.0;
+            let size = 0.3 + rng.unit();
+            let a = scatter(&mut rng, na, origin, 5.0, size, lattice);
+            let b = scatter(&mut rng, nb, origin + shift, 5.0, size, lattice);
+            let want = every_pair(&a, &b, gap);
+            let mut budget = usize::MAX;
+            let got = find_close_pair(&a, &b, gap, &mut budget, |(x, _), (y, _)| corners_meet(x, y, gap));
+            assert_eq!(got, Some(want), "case {case}: {na} x {nb} triangles, gap {gap}, lattice {lattice}");
+            if want {
+                met += 1;
+            } else {
+                apart += 1;
+            }
+        }
+        assert!(met > 50 && apart > 50, "the cases should cover both answers: {met} meet, {apart} apart");
+        let mut budget = usize::MAX;
+        assert_eq!(find_close_pair(&[], &flat_grid(2, 0.0), 0.1, &mut budget, |_, _| true), Some(false));
+        assert_eq!(find_close_pair(&flat_grid(2, 0.0), &[], 0.1, &mut budget, |_, _| true), Some(false));
+    }
+
+    /// Two flat meshes of 16,200 triangles each, one above the other: all their boxes overlap, which
+    /// is 262 million pairs and far past a budget of 40 million.
+    #[test]
+    fn big_meshes_with_overlapping_boxes_are_searched_without_trying_every_pair() {
+        let (low, high) = (flat_grid(90, 0.0), flat_grid(90, 0.5));
+        assert!(low.len() * high.len() > 200_000_000);
+        let mut budget = 40_000_000;
+        assert_eq!(find_close_pair(&low, &high, 0.1, &mut budget, |(x, _), (y, _)| corners_meet(x, y, 0.1)), Some(false));
+        let mut budget = 40_000_000;
+        assert_eq!(find_close_pair(&low, &high, 0.6, &mut budget, |(x, _), (y, _)| corners_meet(x, y, 0.6)), Some(true));
+    }
+
+    /// A pile of triangles in one place is still a product of pairs: the budget refuses it.
+    #[test]
+    fn a_pile_of_triangles_in_one_place_runs_out_of_budget() {
+        let pile = |height: f64| {
+            let up = Vec3::Z * height;
+            vec![boxed([up, Vec3::X + up, Vec3::Y + up]); 7000]
+        };
+        let (low, high) = (pile(0.0), pile(0.5));
+        let mut budget = 40_000_000;
+        // Their boxes are 0.5 apart, so within a gap of 1 every pair is a candidate: 49 million tests.
+        assert_eq!(find_close_pair(&low, &high, 1.0, &mut budget, |_, _| false), None);
+    }
+
+    /// The grid's own work counts: a budget too small to file a big mesh is refused even when no
+    /// pair is ever tested.
+    #[test]
+    fn filing_and_visiting_cells_count_against_the_budget() {
+        let (low, high) = (flat_grid(30, 0.0), flat_grid(30, 5.0));
+        let mut budget = 100;
+        assert_eq!(find_close_pair(&low, &high, 0.1, &mut budget, |_, _| false), None);
+        let mut budget = 10_000_000;
+        assert_eq!(find_close_pair(&low, &high, 0.1, &mut budget, |_, _| false), Some(false));
+        assert!(budget < 10_000_000, "work was counted");
+    }
+}