diff --git a/crates/mesh/src/tet/complex.rs b/crates/mesh/src/tet/complex.rs index 726c41a..7bdfee6 100644 --- a/crates/mesh/src/tet/complex.rs +++ b/crates/mesh/src/tet/complex.rs @@ -1033,6 +1033,17 @@ impl TetComplex { if self.curve(vs[k], vs[l]).is_some() { return Ok(false); } + // Each neighbour gains a boundary face (both, when they are one tetrahedron); one left + // with every face on the boundary would float free of the mesh, joined to it by its + // vertices alone (a singular stiffness). + for m in [k, l] { + let n = ns[m]; + let gained = if ns[k] == ns[l] { 2 } else { 1 }; + let open = self.adj(n)?.iter().filter(|&&x| x == NONE).count(); + if open + gained >= 4 { + return Ok(false); + } + } let raw = self.journal.recording.then(|| Box::new(self.snapshot(&[t, ns[k], ns[l]], &vs))); for m in [k, l] { let n = ns[m]; diff --git a/crates/mesh/src/tet/mod.rs b/crates/mesh/src/tet/mod.rs index 577bee8..a877b6a 100644 --- a/crates/mesh/src/tet/mod.rs +++ b/crates/mesh/src/tet/mod.rs @@ -559,6 +559,26 @@ pub fn quality_mesh(surface: &BoundarySurface, cdt: &TetMesh, analytic: Vec Vec { + if tets.len() < 2 { + return Vec::new(); + } + let face = |t: &[u32; 4], k: usize| { + let mut f = [t[(k + 1) % 4], t[(k + 2) % 4], t[(k + 3) % 4]]; + f.sort_unstable(); + f + }; + let mut count: std::collections::BTreeMap<[u32; 3], u32> = std::collections::BTreeMap::new(); + for t in tets { + for k in 0..4 { + *count.entry(face(t, k)).or_default() += 1; + } + } + tets.iter().enumerate().filter(|(_, t)| (0..4).all(|k| count.get(&face(t, k)).copied().unwrap_or(0) < 2)).map(|(i, _)| i as u32 + 1).collect() +} + /// The constrained Delaunay tetrahedralization itself as a [`QualityMesh`] (no quality stage), /// through the same validity gate: every tetrahedron admissible, the volume equal to what its /// boundary encloses. @@ -576,9 +596,18 @@ pub fn raw_mesh(surface: &BoundarySurface, cdt: &TetMesh) -> Result Ok(QualityMesh { points: cdt.points.clone(), elements: cdt.tets.iter().map(|t| t.to_vec()).collect(), kind: ElementKind::Tet4, boundary, report }) } -/// The validity gate of a linear mesh: every tetrahedron admissible, and the volume equal to -/// what its boundary encloses within 1e-9 (no overlap, no gap). +/// The validity gate of a linear mesh: every tetrahedron admissible, joined to the others by a +/// face (when there are several), and the volume equal to what its boundary encloses within +/// 1e-9 (no overlap, no gap). pub fn linear_gate(points: &[[f64; 3]], tets: &[[u32; 4]], boundary_volume: f64) -> Result { + let detached = detached_tets(tets); + if !detached.is_empty() { + let list = detached.iter().take(10).map(u32::to_string).collect::>().join(", "); + return Err(Error::Failed(format!( + "mesh validity gate: {} tetrahedra share no face with the rest of the mesh ({list}): held by vertices alone, they would make the stiffness singular", + detached.len() + ))); + } let coords = |t: &[u32; 4]| -> Result<[[f64; 3]; 4]> { let p = |v: u32| points.get(v as usize).copied().ok_or_else(|| Error::Failed("gate: vertex id".into())); Ok([p(t[0])?, p(t[1])?, p(t[2])?, p(t[3])?]) diff --git a/crates/mesh/src/tet/tests.rs b/crates/mesh/src/tet/tests.rs index 7054bba..9ad065c 100644 --- a/crates/mesh/src/tet/tests.rs +++ b/crates/mesh/src/tet/tests.rs @@ -454,3 +454,46 @@ mod determinism { } } } + +/// A cap tetrahedron (two boundary faces of one B-rep face) whose two other neighbours each have +/// three boundary faces: deleting it would leave them joined to nothing by a face, floating free +/// (a singular stiffness: the dome-on-a-plate mesh on Linux, #438), so it stays. +#[test] +fn deleting_a_cap_never_leaves_a_neighbour_floating() { + let p: Vec<[f64; 3]> = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.3, 0.3, -1.0], [1.0, 1.0, 1.0]]; + let pts: Vec = p.iter().map(|c| Point3::try_from(*c).unwrap()).collect(); + let positive = |t: [u32; 4]| if measure::volume6(&t.map(|v| pts[v as usize].coords())) > 0.0 { t } else { [t[1], t[0], t[2], t[3]] }; + // The cap, and a tetrahedron on each of its faces 012 and 123. + let tets = [positive([0, 1, 2, 3]), positive([0, 1, 2, 4]), positive([1, 2, 3, 5])]; + let key = |f: [u32; 3]| { + let mut f = f; + f.sort_unstable(); + f + }; + let mut boundary = Vec::new(); + for (t, vs) in tets.iter().enumerate() { + for k in 0..4 { + let f = key(face_toward(*vs, k)); + let shared = tets.iter().enumerate().any(|(u, ws)| u != t && (0..4).any(|j| key(face_toward(*ws, j)) == f)); + if !shared { + boundary.push((t as u32, k, 0)); + } + } + } + let mut cx = TetComplex::new(&p, vec![VertexKind::Corner; 6], &tets, &boundary, &[], 100, 100).unwrap(); + assert!(!cx.delete_boundary_tet(0).unwrap(), "the cap was deleted, leaving its neighbours floating"); + assert_eq!(cx.live_tets().count(), 3); + cx.validate().unwrap(); +} + +/// The validity gate refuses a tetrahedron that shares no face with the rest of the mesh, even +/// when the volume adds up: it says why instead of leaving the solver to find it singular. +#[test] +fn the_gate_refuses_a_detached_tetrahedron() { + let p = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 1.0, 1.0], [2.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 2.0]]; + let tets = [[0, 1, 2, 3], [4, 5, 6, 7]]; + let v = 2.0 / 6.0; + let e = super::linear_gate(&p, &tets, v).unwrap_err().to_string(); + assert!(e.contains("share no face") && e.contains("(1, 2)"), "{e}"); + assert!(super::linear_gate(&p, &tets[..1], v / 2.0).is_ok()); +}