From 9e38fe866d0b48a14ce51202bd68e65f4af90256 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 21:57:24 -0400 Subject: [PATCH] mesh: conformity check finds cracks (free faces, open edges, duplicate nodes, mid-side mismatches) `quality::conformity` checks a solid mesh for cracks: the free faces must close up (every corner edge on exactly two), no two distinct nodes may coincide, and every element edge has one mid-side node. The hex tests assert it for every method, hex20 included, and the crank arm sample's 7274 hex20 is conforming (#260: the disjoint lines are a display issue). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/engine/src/sample.rs | 17 ++ crates/mesh/src/hex/tests.rs | 36 +++++ crates/mesh/src/quality/conformity.rs | 219 ++++++++++++++++++++++++++ crates/mesh/src/quality/mod.rs | 2 + 4 files changed, 274 insertions(+) create mode 100644 crates/mesh/src/quality/conformity.rs diff --git a/crates/engine/src/sample.rs b/crates/engine/src/sample.rs index ca1f72b..b6c5004 100644 --- a/crates/engine/src/sample.rs +++ b/crates/engine/src/sample.rs @@ -817,6 +817,23 @@ mod tests { assert!((reaction["sets"][0]["force"][1].as_f64().unwrap() - PEDAL_LOAD).abs() < 1e-6, "{reaction}"); } + /// The conformity of `s`'s mesh ([`fieldcraft_mesh::quality::conformity`]). + fn conformity(s: &Session) -> fieldcraft_mesh::quality::Conformity { + let nodes: Vec<(u32, [f64; 3])> = s.model.nodes().iter().map(|(l, n)| (*l, n.xyz)).collect(); + let elements = s.model.elements().iter().map(|(l, e)| (*l, e.kind, e.nodes.as_slice())); + fieldcraft_mesh::quality::conformity(&nodes, elements, 1e-6).unwrap() + } + + /// The crank arm swept all-hex20 has no crack (#260): its free faces close up, no two nodes + /// coincide and every element edge has one mid-side node. + #[test] + fn the_crank_arm_hex20_mesh_is_conforming() { + let s = open("crank-arm"); + let c = conformity(&s); + assert_eq!(s.model.elements().len(), 7274); + assert!(c.is_conforming(), "{:?} {:?} {:?}", c.open_edges.first(), c.coincident.first(), c.midside_mismatches.first()); + } + /// The mirrored crank opens as one conforming all-hex20 mesh: its two halves alike, joined /// on the symmetry plane, every element valid, and the spindle and pedal faces found on both /// halves (the image's faces carry the names of the half they are the image of). diff --git a/crates/mesh/src/hex/tests.rs b/crates/mesh/src/hex/tests.rs index 28d34ff..71aa3cb 100644 --- a/crates/mesh/src/hex/tests.rs +++ b/crates/mesh/src/hex/tests.rs @@ -61,8 +61,19 @@ pub(crate) fn prism(poly: &[[f64; 2]], z0: f64, z1: f64) -> BoundarySurface { polyhedron(&pts, &faces) } +/// The mesh has no crack ([`crate::quality::conformity`]): its free faces close up, no two nodes +/// coincide, one mid-side node per element edge. +fn assert_conforming(m: &HexMesh) { + let nodes: Vec<(u32, [f64; 3])> = m.nodes.iter().enumerate().map(|(i, p)| (i as u32, *p)).collect(); + let elements = (0..m.element_count()).filter_map(|e| Some((e as u32, m.kind, m.element(e)?))); + let c = crate::quality::conformity(&nodes, elements, 1e-9).unwrap(); + assert!(c.is_conforming(), "{:?} {:?} {:?}", c.open_edges.first(), c.coincident.first(), c.midside_mismatches.first()); + assert_eq!(c.free_faces.len(), m.boundary.len()); +} + fn check(m: &HexMesh, volume: f64) { assert!(m.report.gate.passed()); + assert_conforming(m); assert!((m.report.volume - volume).abs() < 1e-9 * volume, "{} vs {volume}", m.report.volume); assert!(m.report.min_scaled_jacobian > 0.0); // Every boundary face is tagged and the boundary faces are exactly the hex faces used once. @@ -146,6 +157,7 @@ fn mapped_block_with_a_curved_top() { // Top nodes lie on the tessellated top. let quadratic = hex_mesh(&s, &[], &UniformSize(0.5), &HexOptions { quadratic: true, ..opts }).unwrap(); assert!(quadratic.report.gate.passed()); + assert_conforming(&quadratic); assert!((quadratic.report.volume - s.volume()).abs() < 1e-3 * s.volume(), "{} vs {}", quadratic.report.volume, s.volume()); } @@ -264,6 +276,29 @@ fn an_l_prism_sweeps_one_to_one() { check(&m, 40.0); } +/// Hex20 meshes are conforming whichever method makes them (#260), on a curved body too: the +/// mid-side nodes put on the curved top are shared by the elements on both sides of an edge. +#[test] +fn hex20_meshes_are_conforming_by_every_method() { + let l = [[0.0, 0.0], [6.0, 0.0], [6.0, 2.0], [2.0, 2.0], [2.0, 6.0], [0.0, 6.0]]; + let s = prism(&l, 0.0, 2.0); + for method in [HexMethod::Extrude, HexMethod::Sweep] { + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { method, quadratic: true, ..HexOptions::default() }).unwrap(); + assert_eq!(m.kind, ElementKind::Hex20); + assert_conforming(&m); + } + let curved = curved_block(); + for method in [HexMethod::Extrude, HexMethod::Mapped] { + let m = hex_mesh(&curved, &[], &UniformSize(0.5), &HexOptions { method, quadratic: true, ..HexOptions::default() }).unwrap(); + assert_conforming(&m); + } + let (pipe, faces) = pipe(1.0, 2.0, 2.0, 16); + for method in [HexMethod::Extrude, HexMethod::Revolve] { + let m = hex_mesh(&pipe, &faces, &UniformSize(0.5), &HexOptions { method, quadratic: true, ..HexOptions::default() }).unwrap(); + assert_conforming(&m); + } +} + #[test] fn unsweepable_bodies_say_why_and_hostile_input_is_an_error() { // A tetrahedron: every face is slanted to every other. @@ -369,6 +404,7 @@ fn boxes_side_by_side_conform() { // Corner nodes plus one node per lattice edge: 9·9·5 + (8·9·5)·2 + 9·9·4. assert_eq!(q.mesh.nodes.len(), 405 + 720 + 324); assert!(distinct_positions(&q.mesh)); + assert_conforming(&q.mesh); } #[test] diff --git a/crates/mesh/src/quality/conformity.rs b/crates/mesh/src/quality/conformity.rs new file mode 100644 index 0000000..a825e15 --- /dev/null +++ b/crates/mesh/src/quality/conformity.rs @@ -0,0 +1,219 @@ +//! Whether a solid mesh is conforming: the checks that find cracks. +//! +//! A mesh of solids without cracks has +//! - free faces (element faces used by one element) that form a closed 2-manifold: every corner +//! edge of a free face lies on exactly two free faces ([`Conformity::open_edges`]), +//! - no two distinct nodes at the same place ([`Conformity::coincident`]): a duplicate node +//! splits the elements around it, and the faces between them become internal free faces, +//! - one mid-side node per element edge for quadratic kinds ([`Conformity::midside_mismatches`]): +//! two elements on the same corner edge with different mid-side nodes are joined at the +//! corners and open along the edge. +//! +//! The report is deterministic: every list is sorted. + +use fieldcraft_model::ElementKind; + +use crate::{Error, Result}; + +/// What [`conformity`] found. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Conformity { + /// Free faces, as (element, face number from 1 as [`ElementKind::faces`] numbers them). + pub free_faces: Vec<(u32, u8)>, + /// Corner edges (sorted node pair) of free faces that are not on exactly two free faces. + pub open_edges: Vec<[u32; 2]>, + /// Pairs of distinct nodes (lower first) closer than the tolerance. + pub coincident: Vec<[u32; 2]>, + /// Corner edges (sorted node pair) that elements give different mid-side nodes. + pub midside_mismatches: Vec<[u32; 2]>, +} + +impl Conformity { + /// No open edges, no coincident nodes and no mismatched mid-side nodes. + pub fn is_conforming(&self) -> bool { + self.open_edges.is_empty() && self.coincident.is_empty() && self.midside_mismatches.is_empty() + } +} + +/// Check the solid elements (`(label, kind, nodes)`; shells and line elements are skipped) over +/// the nodes `(label, xyz)`: nodes closer than `tolerance` count as coincident. Elements that +/// name a node missing from `nodes` are an error, as are non-finite coordinates or tolerance. +pub fn conformity<'a>( + nodes: &[(u32, [f64; 3])], + elements: impl IntoIterator, + tolerance: f64, +) -> Result { + if !(tolerance >= 0.0 && tolerance.is_finite()) { + return Err(Error::Invalid(format!("the tolerance must be a non-negative number, got {tolerance}"))); + } + let mut labels: Vec = nodes.iter().map(|n| n.0).collect(); + labels.sort_unstable(); + // (sorted corners, corners in order, element, face), triangles padded with `u32::MAX`. + let mut faces: Vec<([u32; 4], [u32; 4], u32, u8)> = Vec::new(); + // (corner pair, mid-side node) of every edge of a quadratic element. + let mut mids: Vec<([u32; 2], u32)> = Vec::new(); + for (label, kind, conn) in elements { + if kind.dimension() != 3 { + continue; + } + if conn.len() != kind.node_count() { + return Err(Error::Invalid(format!("element {label} has {} nodes; a {} has {}", conn.len(), kind.name(), kind.node_count()))); + } + if let Some(n) = conn.iter().find(|n| labels.binary_search(n).is_err()) { + return Err(Error::Invalid(format!("element {label} names node {n}, which does not exist"))); + } + let node = |i: usize| conn.get(i).copied().unwrap_or(u32::MAX); + for (k, f) in kind.faces().iter().enumerate() { + let mut ordered = [u32::MAX; 4]; + for (slot, &i) in ordered.iter_mut().zip(f.iter().take(face_corners(f.len()))) { + *slot = node(i); + } + let mut key = ordered; + key.sort_unstable(); + faces.push((key, ordered, label, u8::try_from(k + 1).unwrap_or(u8::MAX))); + } + if kind.is_quadratic() { + for e in kind.edges() { + if let (Some(&a), Some(&b), Some(&m)) = (e.first(), e.get(1), e.get(2)) { + let (a, b) = (node(a), node(b)); + mids.push(([a.min(b), a.max(b)], node(m))); + } + } + } + } + faces.sort_unstable(); + let mut out = Conformity::default(); + let mut rim: Vec<[u32; 2]> = Vec::new(); + for run in faces.chunk_by(|a, b| a.0 == b.0) { + let [(_, ordered, element, face)] = run else { continue }; + out.free_faces.push((*element, *face)); + let ordered: Vec = ordered.iter().copied().filter(|&n| n != u32::MAX).collect(); + for (j, &a) in ordered.iter().enumerate() { + let b = ordered.get((j + 1) % ordered.len()).copied().unwrap_or(a); + rim.push([a.min(b), a.max(b)]); + } + } + out.free_faces.sort_unstable(); + rim.sort_unstable(); + out.open_edges = rim.chunk_by(|a, b| a == b).filter(|r| r.len() != 2).filter_map(|r| r.first().copied()).collect(); + mids.sort_unstable(); + mids.dedup(); + out.midside_mismatches = mids.chunk_by(|a, b| a.0 == b.0).filter(|r| r.len() > 1).filter_map(|r| r.first().map(|m| m.0)).collect(); + out.coincident = coincident(nodes, tolerance)?; + Ok(out) +} + +/// Corners of a face with `n` node positions: a quadratic face lists its mid-side nodes after +/// its corners, as many of each. +fn face_corners(n: usize) -> usize { + match n { + 6 | 8 => n / 2, + n => n, + } +} + +/// Pairs of distinct nodes closer than `tol`: nodes are bucketed in a grid of cells `tol` wide +/// (sorted, so the work stays O(n log n) on any input) and compared with the 27 cells around. +fn coincident(nodes: &[(u32, [f64; 3])], tol: f64) -> Result> { + if let Some((l, _)) = nodes.iter().find(|(_, p)| !p.iter().all(|v| v.is_finite())) { + return Err(Error::Invalid(format!("node {l} has a non-finite coordinate"))); + } + // A zero tolerance still needs a cell size; exact copies share a cell at any size. + let h = if tol > 0.0 { tol } else { 1.0 }; + let cell = |p: &[f64; 3]| p.map(|v| (v / h).floor() as i64); + let mut grid: Vec<([i64; 3], u32, [f64; 3])> = nodes.iter().map(|(l, p)| (cell(p), *l, *p)).collect(); + grid.sort_by(|a, b| a.0.cmp(&b.0).then(a.1.cmp(&b.1))); + let mut out = Vec::new(); + for (c, la, pa) in &grid { + for d in 0..27i64 { + let near = [c[0].saturating_add(d % 3 - 1), c[1].saturating_add((d / 3) % 3 - 1), c[2].saturating_add(d / 9 - 1)]; + let start = grid.partition_point(|g| g.0 < near); + for (_, lb, pb) in grid.get(start..).unwrap_or(&[]).iter().take_while(|g| g.0 == near) { + let d2: f64 = (0..3).map(|k| (pa[k] - pb[k]).powi(2)).sum(); + if la < lb && d2 <= tol * tol { + out.push([*la, *lb]); + } + } + } + } + out.sort_unstable(); + out.dedup(); + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Two unit hex8 cells side by side along x (nodes 1–12), elements 1 and 2. + fn pair() -> (Vec<(u32, [f64; 3])>, Vec>) { + let mut nodes = Vec::new(); + for k in 0..2 { + for j in 0..2 { + for i in 0..3 { + nodes.push((1 + i + 3 * j + 6 * k, [f64::from(i), f64::from(j), f64::from(k)])); + } + } + } + let hex = |i: u32| vec![1 + i, 2 + i, 5 + i, 4 + i, 7 + i, 8 + i, 11 + i, 10 + i]; + (nodes, vec![hex(0), hex(1)]) + } + + fn check(nodes: &[(u32, [f64; 3])], elements: &[Vec], kind: ElementKind) -> Conformity { + conformity(nodes, elements.iter().enumerate().map(|(e, n)| (e as u32 + 1, kind, n.as_slice())), 1e-9).unwrap() + } + + #[test] + fn a_conforming_pair_has_a_closed_skin() { + let (nodes, elements) = pair(); + let c = check(&nodes, &elements, ElementKind::Hex8); + assert!(c.is_conforming(), "{c:?}"); + assert_eq!(c.free_faces.len(), 10); + assert!(!c.free_faces.contains(&(1, 4)) && !c.free_faces.contains(&(2, 6)), "the shared face is inside"); + } + + #[test] + fn a_duplicated_node_is_a_crack() { + let (mut nodes, mut elements) = pair(); + // Element 2 takes its own copy of node 2 (at x = 1). + nodes.push((99, [1.0, 0.0, 0.0])); + elements[1][0] = 99; + let c = check(&nodes, &elements, ElementKind::Hex8); + assert_eq!(c.coincident, vec![[2, 99]]); + assert!(!c.open_edges.is_empty(), "the faces beside the copy no longer close up"); + assert_eq!(c.free_faces.len(), 12, "the faces between the elements are free"); + } + + #[test] + fn different_mid_side_nodes_on_one_edge_are_a_crack() { + let (mut nodes, elements) = pair(); + let mut next = 100; + let mut quadratic = Vec::new(); + for (e, n) in elements.iter().enumerate() { + let mut q = n.clone(); + for edge in ElementKind::Hex20.edges() { + next += 1; + let (a, b) = (n[edge[0]], n[edge[1]]); + let (pa, pb) = (nodes.iter().find(|x| x.0 == a).unwrap().1, nodes.iter().find(|x| x.0 == b).unwrap().1); + nodes.push((next, [0.5 * (pa[0] + pb[0]), 0.5 * (pa[1] + pb[1]), 0.5 * (pa[2] + pb[2]) + e as f64 * 1e-3])); + q.push(next); + } + quadratic.push(q); + } + let c = check(&nodes, &quadratic, ElementKind::Hex20); + // The four edges of the shared face got a mid-side node from each element. + assert_eq!(c.midside_mismatches.len(), 4); + assert!(c.open_edges.is_empty(), "the corners still close up"); + } + + #[test] + fn hostile_input_is_an_error() { + let (nodes, elements) = pair(); + assert!(conformity(&nodes, std::iter::once((1, ElementKind::Hex8, &elements[0][..4])), 0.0).is_err()); + assert!(conformity(&nodes, std::iter::once((1, ElementKind::Hex8, &[1, 2, 3, 4, 5, 6, 7, 1000][..])), 0.0).is_err()); + assert!(conformity(&nodes, std::iter::empty(), f64::NAN).is_err()); + assert!(conformity(&[(1, [f64::INFINITY, 0.0, 0.0])], std::iter::empty(), 0.0).is_err()); + let far = [(1, [1e300, 0.0, 0.0]), (2, [1e300, 0.0, 0.0]), (3, [-1e300, 0.0, 0.0])]; + assert_eq!(conformity(&far, std::iter::empty(), 1e-300).unwrap().coincident, vec![[1, 2]]); + } +} diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index b81f4ac..8f7b527 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -25,12 +25,14 @@ //! their integration points, so whatever passes the gate passes the solvers. mod bernstein; +pub mod conformity; pub mod jacobian; mod report; pub(crate) mod shape; use fieldcraft_model::{ElementKind, admissible_jacobian, check_coordinates, tet4_jacobian}; +pub use conformity::{Conformity, conformity}; pub use jacobian::{JacobianCheck, Validity}; pub use report::{ GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_checks, validity_gate_elements,