From 5cee94f947052af1b9beb86c4b0ca0b7158f7189 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 23:01:46 -0400 Subject: [PATCH] geometry: check bodies for short edges, small and sliver faces, near faces Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/geometry/src/check.rs | 404 +++++++++++++++++++++++++++++++++++ crates/geometry/src/edges.rs | 4 +- crates/geometry/src/lib.rs | 1 + 3 files changed, 407 insertions(+), 2 deletions(-) create mode 100644 crates/geometry/src/check.rs diff --git a/crates/geometry/src/check.rs b/crates/geometry/src/check.rs new file mode 100644 index 0000000..5f330b8 --- /dev/null +++ b/crates/geometry/src/check.rs @@ -0,0 +1,404 @@ +//! Geometry check before meshing: the features of a body that are small against the element +//! size an analyst intends to mesh it with, and so will force tiny or distorted elements. +//! +//! Four kinds, each measured on a fine tessellation of the B-rep (its edge polylines and its +//! triangles tagged by face) and compared with a fraction of the target size: +//! - **short edges**: an edge shorter than `short_edge × size`; +//! - **small faces**: a face whose area is below `small_face × size²`; +//! - **sliver faces**: a face that is not small but narrow, its mean width `2 × area / perimeter` +//! (the width of a long strip, the radius of a disc) below `sliver × size`; +//! - **near faces**: two nearly parallel faces closer than `gap × size` without touching: a gap +//! between bodies or cells that almost meet, faces that overlap, or a wall thinner than the gap. +//! +//! Every finding names the body, its persistent face and edge names, the measured value, the +//! limit it fell below and a point to look at. The check reads the geometry and changes nothing. + +use std::collections::{BTreeMap, BTreeSet, HashMap}; + +use crate::edges::EDGE_TOLERANCE; +use crate::{Error, Result, Shape}; + +/// Most triangles the check looks at, over every body. +pub const MAX_TRIANGLES: usize = 1_000_000; + +/// Most grid-cell insertions of the near-face search (a triangle goes into every cell its box +/// covers). +const MAX_CELL_ENTRIES: usize = 20_000_000; + +/// The thresholds of a check, as fractions of the target element size. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct CheckOptions { + /// The target element size (mm). + pub size: f64, + pub short_edge: f64, + pub small_face: f64, + pub sliver: f64, + pub gap: f64, + /// Faces whose normals differ from parallel (or opposite) by more than this (radians) are + /// never near faces. + pub parallel_angle: f64, +} + +impl CheckOptions { + /// The default thresholds for `size`. + pub fn new(size: f64) -> Self { + CheckOptions { size, short_edge: 0.3, small_face: 0.1, sliver: 0.3, gap: 0.3, parallel_angle: 10f64.to_radians() } + } +} + +/// What a finding is about. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum Kind { + ShortEdge, + SmallFace, + SliverFace, + NearFaces, +} + +impl Kind { + pub const ALL: [Kind; 4] = [Kind::ShortEdge, Kind::SmallFace, Kind::SliverFace, Kind::NearFaces]; + + pub fn name(self) -> &'static str { + match self { + Kind::ShortEdge => "short_edge", + Kind::SmallFace => "small_face", + Kind::SliverFace => "sliver_face", + Kind::NearFaces => "near_faces", + } + } +} + +/// One feature below its threshold. +#[derive(Debug, Clone, PartialEq)] +pub struct Finding { + pub kind: Kind, + pub body: String, + /// The edge, for a short edge. + pub edge: Option, + /// The face (small and sliver faces), the faces of the edge (short edges), or the two faces + /// (near faces: the first on `body`, the second on `other_body`). + pub faces: Vec, + /// The body of the second face, for near faces. + pub other_body: Option, + /// The measure: length (mm), area (mm²), width (mm) or distance (mm). + pub value: f64, + /// The limit `value` fell below, in the same unit. + pub limit: f64, + /// A point on the feature to look at. + pub at: [f64; 3], +} + +/// Checks `bodies` against `options`. Findings come grouped by kind, worst first within a kind. +pub fn check(bodies: &[&Shape], options: &CheckOptions) -> Result> { + let o = options; + let fractions = [o.short_edge, o.small_face, o.sliver, o.gap]; + if !(o.size > 0.0 && o.size.is_finite()) || fractions.iter().any(|f| !(*f >= 0.0 && f.is_finite())) { + return Err(Error::Invalid("the check needs a positive size and thresholds that are finite and not negative".into())); + } + if !(o.parallel_angle >= 0.0 && o.parallel_angle < std::f64::consts::FRAC_PI_2) { + return Err(Error::Invalid("the parallel angle must be at least 0 and below 90°".into())); + } + let mut findings = Vec::new(); + let mut triangles = Vec::new(); + let mut adjacent = BTreeSet::new(); + let mut touch = 0.0f64; + for (b, shape) in bodies.iter().enumerate() { + let tolerance = (shape.size()? * EDGE_TOLERANCE).max(1e-4); + touch = touch.max(2.0 * tolerance); + let mesh = shape.tessellate(tolerance)?; + if triangles.len() + mesh.triangles.len() > MAX_TRIANGLES { + return Err(Error::Invalid(format!("the bodies have more than {MAX_TRIANGLES} triangles to check"))); + } + let edges = shape.edges_of(&mesh); + // Short edges, and each face's perimeter (a seam, both sides on one face, bounds nothing). + let mut perimeter: BTreeMap = BTreeMap::new(); + for e in &edges { + let length = e.length(); + let distinct: BTreeSet = e.faces.iter().copied().collect(); + if !(e.faces.len() == 2 && distinct.len() == 1) { + for &f in &distinct { + *perimeter.entry(f).or_default() += length; + } + } + if let [f, g] = e.faces[..] { + adjacent.insert((b, f.min(g), f.max(g))); + } + let limit = o.short_edge * o.size; + if length < limit { + findings.push(Finding { + kind: Kind::ShortEdge, + body: shape.name.clone(), + edge: Some(e.name.clone()), + faces: distinct.iter().map(|&f| shape.face_label(f)).collect(), + other_body: None, + value: length, + limit, + at: e.point_at(0.5).unwrap_or_default(), + }); + } + } + // Face areas, area centroids and the triangle centres to place them by. + let mut faces: BTreeMap = BTreeMap::new(); + for (t, &f) in mesh.triangles.iter().zip(&mesh.tri_face) { + let corner = |i: u32| mesh.positions.get(i as usize).copied(); + let (Some(p), Some(q), Some(r)) = (corner(t[0]), corner(t[1]), corner(t[2])) else { continue }; + let n = cross(sub(q, p), sub(r, p)); + let area = 0.5 * norm(n); + let centre = [0, 1, 2].map(|k| (p[k] + q[k] + r[k]) / 3.0); + let sum = faces.entry(f).or_default(); + sum.area += area; + sum.moment = add(sum.moment, scale(centre, area)); + sum.centres.push(centre); + if area > 0.0 && area.is_finite() { + triangles.push(Tri { body: b, face: f, corners: [p, q, r], normal: scale(n, 0.5 / area) }); + } + } + for (&f, sum) in &faces { + let at = sum.location(); + let small = o.small_face * o.size * o.size; + if sum.area < small { + findings.push(face_finding(Kind::SmallFace, shape, f, sum.area, small, at)); + continue; + } + let Some(p) = perimeter.get(&f).copied().filter(|p| *p > 0.0) else { continue }; + let width = 2.0 * sum.area / p; + let limit = o.sliver * o.size; + if width < limit { + findings.push(face_finding(Kind::SliverFace, shape, f, width, limit, at)); + } + } + } + findings.extend(near_faces(bodies, &triangles, &adjacent, touch, o)?); + findings.sort_by(|a, b| a.kind.cmp(&b.kind).then(a.value.total_cmp(&b.value))); + Ok(findings) +} + +fn face_finding(kind: Kind, shape: &Shape, f: u32, value: f64, limit: f64, at: [f64; 3]) -> Finding { + Finding { kind, body: shape.name.clone(), edge: None, faces: vec![shape.face_label(f)], other_body: None, value, limit, at } +} + +#[derive(Default)] +struct FaceSum { + area: f64, + moment: [f64; 3], + centres: Vec<[f64; 3]>, +} + +impl FaceSum { + /// The triangle centre nearest the area centroid: a point on the face (the centroid of a + /// curved or ring-shaped face lies off it). + fn location(&self) -> [f64; 3] { + let c = if self.area > 0.0 { scale(self.moment, 1.0 / self.area) } else { self.centres.first().copied().unwrap_or_default() }; + self.centres.iter().copied().min_by(|a, b| norm(sub(*a, c)).total_cmp(&norm(sub(*b, c)))).unwrap_or(c) + } +} + +struct Tri { + body: usize, + face: u32, + corners: [[f64; 3]; 3], + /// Unit normal. + normal: [f64; 3], +} + +/// Pairs of nearly parallel faces closer than the gap but farther apart than `touch` (faces +/// within it are the same surface tessellated twice, such as the shared face of two cells), +/// measured from points on each triangle to their perpendicular foot on the other face. +fn near_faces(bodies: &[&Shape], tris: &[Tri], adjacent: &BTreeSet<(usize, u32, u32)>, touch: f64, o: &CheckOptions) -> Result> { + let gap = o.gap * o.size; + if gap <= touch || tris.is_empty() { + return Ok(Vec::new()); + } + // A uniform grid over the triangles' boxes, cells no smaller than the gap. + let mut extents: Vec = tris.iter().map(|t| norm(sub(bbox(t).1, bbox(t).0))).collect(); + extents.sort_by(f64::total_cmp); + let cell = extents.get(extents.len() / 2).copied().unwrap_or(gap).max(gap); + let key = |p: [f64; 3]| p.map(|x| (x / cell).floor() as i64); + let mut grid: HashMap<[i64; 3], Vec> = HashMap::new(); + let mut entries = 0usize; + for (i, t) in tris.iter().enumerate() { + let (lo, hi) = bbox(t); + let (lo, hi) = (key(lo), key(hi)); + let count = (0..3).map(|k| (hi[k] - lo[k] + 1).max(1) as usize).product::(); + entries = entries.saturating_add(count); + if entries > MAX_CELL_ENTRIES { + return Err(Error::Invalid("the near-face search is too large for these bodies; check them one at a time".into())); + } + for x in lo[0]..=hi[0] { + for y in lo[1]..=hi[1] { + for z in lo[2]..=hi[2] { + grid.entry([x, y, z]).or_default().push(i); + } + } + } + } + let parallel = o.parallel_angle.cos(); + // Per face pair (the lower first): the smallest distance and where. + let mut pairs: BTreeMap<((usize, u32), (usize, u32)), (f64, [f64; 3])> = BTreeMap::new(); + for a in tris { + let [p, q, r] = a.corners; + let c = [0, 1, 2].map(|k| (p[k] + q[k] + r[k]) / 3.0); + for s in [c, mid(c, p), mid(c, q), mid(c, r)] { + let (lo, hi) = (key(s.map(|x| x - gap)), key(s.map(|x| x + gap))); + let mut seen = BTreeSet::new(); + for x in lo[0]..=hi[0] { + for y in lo[1]..=hi[1] { + for z in lo[2]..=hi[2] { + for &j in grid.get(&[x, y, z]).into_iter().flatten() { + let Some(b) = tris.get(j) else { continue }; + if !seen.insert(j) || (a.body, a.face) >= (b.body, b.face) { + continue; + } + if a.body == b.body && adjacent.contains(&(a.body, a.face.min(b.face), a.face.max(b.face))) { + continue; + } + if dot(a.normal, b.normal).abs() < parallel { + continue; + } + let foot = closest_on_triangle(s, b.corners); + let d = norm(sub(s, foot)); + // Only a perpendicular foot: a point beside the other face is not over it. + if d <= touch || d >= gap || dot(sub(s, foot), b.normal).abs() < 0.999 * d { + continue; + } + let slot = pairs.entry(((a.body, a.face), (b.body, b.face))).or_insert((d, s)); + if d < slot.0 { + *slot = (d, s); + } + } + } + } + } + } + } + let mut out = Vec::new(); + for (((ba, fa), (bb, fb)), (d, at)) in pairs { + let (Some(sa), Some(sb)) = (bodies.get(ba), bodies.get(bb)) else { continue }; + out.push(Finding { + kind: Kind::NearFaces, + body: sa.name.clone(), + edge: None, + faces: vec![sa.face_label(fa), sb.face_label(fb)], + other_body: Some(sb.name.clone()), + value: d, + limit: gap, + at, + }); + } + Ok(out) +} + +fn bbox(t: &Tri) -> ([f64; 3], [f64; 3]) { + let [p, q, r] = t.corners; + ([0, 1, 2].map(|k| p[k].min(q[k]).min(r[k])), [0, 1, 2].map(|k| p[k].max(q[k]).max(r[k]))) +} + +/// The point of triangle `abc` closest to `p` (by its Voronoi regions: Ericson, *Real-Time +/// Collision Detection*, 2005, §5.1.5). +fn closest_on_triangle(p: [f64; 3], [a, b, c]: [[f64; 3]; 3]) -> [f64; 3] { + let (ab, ac, ap) = (sub(b, a), sub(c, a), sub(p, a)); + let (d1, d2) = (dot(ab, ap), dot(ac, ap)); + if d1 <= 0.0 && d2 <= 0.0 { + return a; + } + let bp = sub(p, b); + let (d3, d4) = (dot(ab, bp), dot(ac, bp)); + if d3 >= 0.0 && d4 <= d3 { + return b; + } + let vc = d1 * d4 - d3 * d2; + if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 { + return add(a, scale(ab, d1 / (d1 - d3))); + } + let cp = sub(p, c); + let (d5, d6) = (dot(ab, cp), dot(ac, cp)); + if d6 >= 0.0 && d5 <= d6 { + return c; + } + let vb = d5 * d2 - d1 * d6; + if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 { + return add(a, scale(ac, d2 / (d2 - d6))); + } + let va = d3 * d6 - d5 * d4; + if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 { + return add(b, scale(sub(c, b), (d4 - d3) / ((d4 - d3) + (d5 - d6)))); + } + let denom = va + vb + vc; + if !(denom.abs() > 0.0) { + return a; + } + add(a, add(scale(ab, vb / denom), scale(ac, vc / denom))) +} + +fn add(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] + b[0], a[1] + b[1], a[2] + b[2]] +} +fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} +fn scale(a: [f64; 3], k: f64) -> [f64; 3] { + a.map(|x| x * k) +} +fn mid(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + scale(add(a, b), 0.5) +} +fn dot(a: [f64; 3], b: [f64; 3]) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} +fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} +fn norm(a: [f64; 3]) -> f64 { + dot(a, a).sqrt() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn closest_point_regions() { + let t = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]; + let near = |p, q: [f64; 3]| norm(sub(closest_on_triangle(p, t), q)) < 1e-12; + assert!(near([0.2, 0.2, 3.0], [0.2, 0.2, 0.0])); + assert!(near([-1.0, -1.0, 0.0], [0.0, 0.0, 0.0])); + assert!(near([0.5, -2.0, 0.0], [0.5, 0.0, 0.0])); + assert!(near([2.0, 2.0, 0.0], [0.5, 0.5, 0.0])); + } + + #[test] + fn hostile_options_are_errors() { + let body = crate::block([0.0; 3], [1.0; 3]).unwrap(); + for o in [ + CheckOptions::new(0.0), + CheckOptions::new(f64::NAN), + CheckOptions { gap: -1.0, ..CheckOptions::new(1.0) }, + CheckOptions { parallel_angle: 2.0, ..CheckOptions::new(1.0) }, + ] { + assert!(check(&[&body], &o).is_err(), "{o:?}"); + } + } + + /// A unit cube at size 1 has nothing small; at size 10 every edge is short and every face small. + #[test] + fn a_cube_is_clean_at_its_own_size() { + let body = crate::block([0.0; 3], [1.0; 3]).unwrap(); + assert!(check(&[&body], &CheckOptions::new(1.0)).unwrap().is_empty()); + let f = check(&[&body], &CheckOptions::new(10.0)).unwrap(); + assert_eq!(f.iter().filter(|f| f.kind == Kind::ShortEdge).count(), 12); + assert_eq!(f.iter().filter(|f| f.kind == Kind::SmallFace).count(), 6); + } + + /// Two blocks 0.1 apart face each other across the gap; a block on the other's face touches it. + #[test] + fn near_faces_across_a_gap_but_not_a_contact() { + let a = crate::block([0.0; 3], [10.0; 3]).unwrap(); + let b = crate::block([0.0, 0.0, 10.1], [10.0, 10.0, 20.0]).unwrap(); + let f = check(&[&a, &b], &CheckOptions::new(2.0)).unwrap(); + let near: Vec<&Finding> = f.iter().filter(|f| f.kind == Kind::NearFaces).collect(); + assert_eq!(near.len(), 1, "{f:?}"); + assert!((near[0].value - 0.1).abs() < 1e-9 && (near[0].at[2] - 10.0).abs() < 1e-9, "{:?}", near[0]); + let c = crate::block([0.0, 0.0, 10.0], [10.0, 10.0, 20.0]).unwrap(); + assert!(check(&[&a, &c], &CheckOptions::new(2.0)).unwrap().is_empty()); + } +} diff --git a/crates/geometry/src/edges.rs b/crates/geometry/src/edges.rs index c7fe92b..a2f43cd 100644 --- a/crates/geometry/src/edges.rs +++ b/crates/geometry/src/edges.rs @@ -128,7 +128,7 @@ fn number_repeats(names: &mut [String], centres: &[[f64; 3]]) { impl Shape { /// The name of face `f`: its persistent name, or its index for a shape without names. - fn face_label(&self, f: u32) -> String { + pub(crate) fn face_label(&self, f: u32) -> String { match self.faces.get(f as usize) { Some(face) if !face.name.is_empty() => face.name.clone(), _ => f.to_string(), @@ -201,7 +201,7 @@ impl Shape { } /// The edges of `mesh` (a tessellation of this shape) with their persistent names. - fn edges_of(&self, mesh: &SurfaceMesh) -> Vec { + pub(crate) fn edges_of(&self, mesh: &SurfaceMesh) -> Vec { let names = self.edge_names(mesh); names .into_iter() diff --git a/crates/geometry/src/lib.rs b/crates/geometry/src/lib.rs index 3a41e65..3dd6e3e 100644 --- a/crates/geometry/src/lib.rs +++ b/crates/geometry/src/lib.rs @@ -17,6 +17,7 @@ #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] +pub mod check; mod design; mod edges; pub mod features;