diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 5021fd8b..7daefc8f 100644 --- a/crates/engine/src/tests.rs +++ b/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; diff --git a/crates/geom/src/mesh.rs b/crates/geom/src/mesh.rs index 792c6232..03aa7cc8 100644 --- a/crates/geom/src/mesh.rs +++ b/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 { + 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, mut test: impl FnMut(&BoxedTriangle, &BoxedTriangle) -> bool, ) -> Option { + 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 { - *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 { diff --git a/crates/geom/src/tests.rs b/crates/geom/src/tests.rs index 317a2697..3929e5ac 100644 --- a/crates/geom/src/tests.rs +++ b/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 { + 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 { + 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"); + } +}