From b6e3a6c3663cf9b794b4380b8cf5d87de8d075b8 Mon Sep 17 00:00:00 2001 From: Mqrius Date: Tue, 13 Oct 2026 22:33:14 +0200 Subject: [PATCH] fix: contact and shell-grouping checks use a shared broad phase (#838) Contact and shell-grouping checks find close triangle pairs through one shared grid instead of testing every pair. Refs #232. --- crates/geom/src/lib.rs | 2 +- crates/geom/src/mesh.rs | 162 +++++++++++++++++++-------- crates/kernel/src/body_lump_tests.rs | 12 +- crates/kernel/src/shell_groups.rs | 44 +++----- crates/kernel/src/tests.rs | 11 ++ 5 files changed, 150 insertions(+), 81 deletions(-) diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 4486ddbb..eb31da69 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -21,7 +21,7 @@ pub mod tspline; mod vec; pub use arcfit::MAX_FIT_PIECES; -pub use mesh::{Aabb3, InsideIndex, Mesh, MeshMeasure}; +pub use mesh::{Aabb3, BoxedTriangle, InsideIndex, Mesh, MeshMeasure, find_close_pair}; pub use plane::Plane; pub use profile::{Loop2, MAX_BEZIER_DEGREE, Region2, Seg2, bezier}; pub use vec::{Vec2, Vec3}; diff --git a/crates/geom/src/mesh.rs b/crates/geom/src/mesh.rs index d0b07add..792c6232 100644 --- a/crates/geom/src/mesh.rs +++ b/crates/geom/src/mesh.rs @@ -80,6 +80,120 @@ fn boxes_meet(a: &Aabb3, b: &Aabb3) -> bool { a.min.x <= b.max.x && b.min.x <= a.max.x && a.min.y <= b.max.y && b.min.y <= a.max.y && a.min.z <= b.max.z && b.min.z <= a.max.z } +/// A triangle with its bounding box, as [`find_close_pair`] takes them. +pub type BoxedTriangle = ([Vec3; 3], Aabb3); + +/// Most grid cells one triangle is filed under; a larger one is tested against every triangle. +const MAX_CELLS_PER_TRIANGLE: usize = 64; + +/// Calls `test` on pairs of triangles, one from `a` and one from `b`, whose boxes come within +/// `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. +pub fn find_close_pair( + a: &[BoxedTriangle], + b: &[BoxedTriangle], + gap: f64, + budget: &mut usize, + mut test: impl FnMut(&BoxedTriangle, &BoxedTriangle) -> bool, +) -> Option { + 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 + && 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 + }; + let mut call = |t: &BoxedTriangle, u: &BoxedTriangle| -> Option { + *budget = budget.checked_sub(1)?; + Some(test(t, u)) + }; + let (Some(first), true) = (b.first(), !a.is_empty()) else { return Some(false) }; + // Cell size: the median extent of `b`'s boxes (so most of them span a few cells), at least `gap`. + let mut extents: Vec = b + .iter() + .map(|(_, bx)| { + let d = bx.size(); + d.x.max(d.y).max(d.z) + }) + .filter(|e| e.is_finite()) + .collect(); + let mid = extents.len() / 2; + let median = if extents.is_empty() { 0.0 } else { *extents.select_nth_unstable_by(mid, f64::total_cmp).1 }; + let origin = first.1.min; + let cell = median.max(gap).max(1e-9); + if !(origin.is_finite() && cell.is_finite()) { + return Some(false); + } + let index = |v: f64, o: f64| ((v - o) / cell).floor() as i64; + let span = |bx: &Aabb3, pad: f64| -> Option<[(i64, i64); 3]> { + let r = [ + (index(bx.min.x - pad, origin.x), index(bx.max.x + pad, origin.x)), + (index(bx.min.y - pad, origin.y), index(bx.max.y + pad, origin.y)), + (index(bx.min.z - pad, origin.z), index(bx.max.z + pad, origin.z)), + ]; + let count = r.iter().try_fold(1usize, |n, (lo, hi)| n.checked_mul(usize::try_from(hi.checked_sub(*lo)?.checked_add(1)?).ok()?))?; + (count <= MAX_CELLS_PER_TRIANGLE).then_some(r) + }; + let mut cells: std::collections::HashMap<(i64, i64, i64), Vec> = std::collections::HashMap::new(); + let mut big = Vec::new(); + for (j, (_, bx)) in b.iter().enumerate() { + match span(bx, gap) { + Some([x, y, z]) => { + for ix in x.0..=x.1 { + for iy in y.0..=y.1 { + for iz in z.0..=z.1 { + cells.entry((ix, iy, iz)).or_default().push(j as u32); + } + } + } + } + None => big.push(j), + } + } + // Per triangle of `b`: the last triangle of `a` tested against it, plus one. + let mut seen = vec![0usize; b.len()]; + for (i, t) in a.iter().enumerate() { + for &j in &big { + if let Some(u) = b.get(j).filter(|u| close(&t.1, &u.1)) { + seen[j] = i + 1; + if call(t, u)? { + return Some(true); + } + } + } + let Some([x, y, z]) = span(&t.1, 0.0) else { + // A large triangle of `a`: against all of `b`. + for (j, u) in b.iter().enumerate() { + if seen[j] != i + 1 && close(&t.1, &u.1) && call(t, u)? { + return Some(true); + } + } + continue; + }; + for ix in x.0..=x.1 { + for iy in y.0..=y.1 { + for iz in z.0..=z.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 { + continue; + } + *mark = i + 1; + if close(&t.1, &u.1) && call(t, u)? { + return Some(true); + } + } + } + } + } + } + Some(false) +} + /// Two triangles that pass through each other (see [`Mesh::crosses`]). fn triangles_cross(t: &[Vec3; 3], u: &[Vec3; 3], tol: f64) -> bool { // Signed distances of one triangle's corners from the other's plane, when they lie more @@ -347,52 +461,8 @@ impl Mesh { if ta.is_empty() || tb.is_empty() { return false; } - // A grid over the shared box holds the other mesh's triangles. - let n = ((tb.len() as f64).cbrt().ceil() as usize).clamp(1, 24); - let size = region.size(); - let cell = |v: f64, lo: f64, s: f64| -> usize { if s > 0.0 { (((v - lo) / s * n as f64).floor().max(0.0) as usize).min(n - 1) } else { 0 } }; - let span = |bx: &Aabb3| -> [(usize, usize); 3] { - [ - (cell(bx.min.x, region.min.x, size.x), cell(bx.max.x, region.min.x, size.x)), - (cell(bx.min.y, region.min.y, size.y), cell(bx.max.y, region.min.y, size.y)), - (cell(bx.min.z, region.min.z, size.z), cell(bx.max.z, region.min.z, size.z)), - ] - }; - let mut cells: Vec> = vec![Vec::new(); n * n * n]; - for (i, (_, bx)) in tb.iter().enumerate() { - let [x, y, z] = span(bx); - for ix in x.0..=x.1 { - for iy in y.0..=y.1 { - for iz in z.0..=z.1 { - if let Some(c) = cells.get_mut((ix * n + iy) * n + iz) { - c.push(i as u32); - } - } - } - } - } - // Per triangle of `o`: the last triangle of `self` tested against it, plus one. - let mut seen = vec![0usize; tb.len()]; - for (i, (t, bx)) in ta.iter().enumerate() { - let [x, y, z] = span(bx); - for ix in x.0..=x.1 { - for iy in y.0..=y.1 { - for iz in z.0..=z.1 { - for j in cells.get((ix * n + iy) * n + iz).into_iter().flatten() { - let (Some(mark), Some((u, ub))) = (seen.get_mut(*j as usize), tb.get(*j as usize)) else { continue }; - if *mark == i + 1 { - continue; - } - *mark = i + 1; - if boxes_meet(bx, ub) && triangles_cross(t, u, tol) { - return true; - } - } - } - } - } - } - false + let mut unlimited = usize::MAX; + find_close_pair(&ta, &tb, 0.0, &mut unlimited, |(t, _), (u, _)| triangles_cross(t, u, tol)).unwrap_or(false) } /// Closest intersection of a ray with the mesh: (distance along `dir`, triangle index). diff --git a/crates/kernel/src/body_lump_tests.rs b/crates/kernel/src/body_lump_tests.rs index 469480b5..19457af3 100644 --- a/crates/kernel/src/body_lump_tests.rs +++ b/crates/kernel/src/body_lump_tests.rs @@ -170,9 +170,15 @@ fn shell_grouping_preserves_large_disconnected_patterns_and_bounds_hostile_work( let mut oversized = mesh.clone(); oversized.positions.resize(shell_groups::MAX_POINTS + 1, Vec3::ZERO); assert!(shell_groups::groups(&[oversized], &[8.0], 1e-8).unwrap_err().to_string().contains("grouping limit")); - let mut repeated = mesh.clone(); - repeated.triangles = mesh.triangles.iter().copied().cycle().take(7000).collect(); - assert!(shell_groups::groups(&[repeated.clone(), repeated], &[8.0, 8.0], 1e-8).unwrap_err().to_string().contains("boundary work")); + // Two stacks of one slanted triangle, in parallel planes 1 µm apart: every pair's boxes + // overlap but none touch, so all 7000² pairs need testing. + let slanted = |dz: f64| Mesh { + positions: vec![Vec3::new(0.0, 0.0, dz), Vec3::new(1.0, 0.0, 1.0 + dz), Vec3::new(0.0, 1.0, 1.0 + dz)], + triangles: vec![[0, 1, 2]; 7000], + ..Mesh::default() + }; + let (repeated, lifted) = (slanted(0.0), slanted(1e-3)); + assert!(shell_groups::groups(&[repeated, lifted], &[8.0, 8.0], 1e-8).unwrap_err().to_string().contains("boundary work")); let mut outer = cube(2.0).tessellate_with(0.01, true).unwrap(); let mut inner = mesh; outer.triangles = outer.triangles.iter().copied().cycle().take(4000).collect(); diff --git a/crates/kernel/src/shell_groups.rs b/crates/kernel/src/shell_groups.rs index cc2c8703..40702e16 100644 --- a/crates/kernel/src/shell_groups.rs +++ b/crates/kernel/src/shell_groups.rs @@ -1,7 +1,7 @@ //! Group disconnected material regions without moving a cavity to another enclosing shell. //! A valid shell hierarchy alternates material and void at every nesting level. -use solidcraft_geom::{Aabb3, Mesh, Vec3}; +use solidcraft_geom::{Aabb3, BoxedTriangle, Mesh, Vec3, find_close_pair}; use crate::{KernelError, Result}; @@ -90,7 +90,7 @@ pub(super) fn meshes_meet(a: &Mesh, b: &Mesh, gap: f64) -> Result { return Ok(false); } // Only triangles near the other body can touch it. - let near = |m: &Mesh, other: Aabb3| -> Result> { + let near = |m: &Mesh, other: Aabb3| -> Result> { let mut out = Vec::new(); for t in &m.triangles { let tri = m.tri(t).ok_or_else(|| failed("invalid contact triangle"))?; @@ -102,16 +102,11 @@ pub(super) fn meshes_meet(a: &Mesh, b: &Mesh, gap: f64) -> Result { Ok(out) }; let (ta, tb) = (near(a, bb)?, near(b, ba)?); - let work = ta.len().checked_mul(tb.len()).ok_or_else(|| failed("contact work exceeds the limit"))?; - if work > MAX_PAIR_WORK { - return Err(failed("contact work exceeds the limit")); - } - for (x, xb) in &ta { - for (y, yb) in &tb { - if !boxes_separate(*xb, *yb, gap) && triangles_meet(x, y, gap) { - return Ok(true); - } - } + let mut budget = MAX_PAIR_WORK; + if find_close_pair(&ta, &tb, gap, &mut budget, |(x, _), (y, _)| triangles_meet(x, y, gap)) + .ok_or_else(|| failed("contact work exceeds the limit"))? + { + return Ok(true); } // No boundary contact: they meet only if one lies wholly inside the other. Ok(b.positions.first().is_some_and(|p| a.contains(*p)) || a.positions.first().is_some_and(|p| b.contains(*p))) @@ -127,7 +122,6 @@ pub(super) fn groups(meshes: &[Mesh], volumes: &[f64], gap: f64) -> Result Result>()); - triangles.push(tris); + triangles.push(tris.into_iter().map(|t| (t, box_of(&t))).collect::>()); bounds.push(mesh.bounds()); } let mut inside = vec![vec![false; meshes.len()]; meshes.len()]; - let mut pair_work = 0usize; + let mut pair_budget = MAX_PAIR_WORK; let mut probes = 0usize; let mut containment_work = 0usize; for (i, mesh) in meshes.iter().enumerate() { @@ -162,21 +155,10 @@ pub(super) fn groups(meshes: &[Mesh], volumes: &[f64], gap: f64) -> Result MAX_PAIR_WORK { - return Err(failed("shell boundary work exceeds the grouping limit")); - } - for (ai, ta) in a.iter().enumerate() { - for (bi, tb) in b.iter().enumerate() { - let (Some(ba), Some(bb)) = (triangle_bounds.get(i).and_then(|x| x.get(ai)), triangle_bounds.get(j).and_then(|x| x.get(bi))) - else { - return Err(failed("missing triangle bounds")); - }; - if !boxes_separate(*ba, *bb, gap) && triangles_meet(ta, tb, gap) { - return Err(failed("shell boundaries intersect or touch; cavity ownership is ambiguous")); - } - } + if find_close_pair(a, b, gap, &mut pair_budget, |(x, _), (y, _)| triangles_meet(x, y, gap)) + .ok_or_else(|| failed("shell boundary work exceeds the grouping limit"))? + { + return Err(failed("shell boundaries intersect or touch; cavity ownership is ambiguous")); } for (inner, outer) in [(i, j), (j, i)] { let ib = bounds.get(inner).ok_or_else(|| failed("missing shell bounds"))?; diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index a8988099..abd2a89a 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -834,6 +834,17 @@ fn sphere_booleans() { assert_eq!((m.merged.faces, m.merged.edges, m.merged.vertices), (7, 13, 9)); } +#[test] +fn similar_spheres_meet() { + // bherbruck/solidcraft#232: the contact test compared every pair of nearby triangles of the + // two fine sphere meshes and gave up past its work limit, so their join was refused. + let a = sphere(Vec3::ZERO, 10.0).unwrap(); + for (x, meet) in [(12.0, true), (4.0, true), (20.0, true), (20.5, false)] { + let b = sphere(Vec3::new(x, 0.0, 0.0), 10.0).unwrap(); + assert_eq!(a.meets(&b).unwrap(), meet, "centres {x} mm apart"); + } +} + #[test] #[ignore] fn debug_convergence() {