diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index 82db91e..e1029df 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -250,6 +250,7 @@ pub(super) fn quality_json(q: &TetQuality) -> Value { "max_dihedral": q.max_dihedral, "min_radius_ratio": q.min_radius_ratio, "mean_radius_ratio": q.mean_radius_ratio, + "min_scaled_jacobian": q.min_scaled_jacobian, "slivers": q.slivers, "inadmissible": q.inadmissible, }) diff --git a/crates/mesh/src/quality/tests.rs b/crates/mesh/src/quality/tests.rs index 113444a..b002884 100644 --- a/crates/mesh/src/quality/tests.rs +++ b/crates/mesh/src/quality/tests.rs @@ -44,6 +44,43 @@ fn regular_tetrahedron_scores_one() { } } +/// The dihedral angles and the scaled Jacobian measure different things on the same corners +/// (#547): a flat tetrahedron is low on both; a needle has every dihedral angle between 60° and +/// 90° but a nearly flat apex corner, so only the scaled Jacobian sees it. A tet10 with straight +/// edges scores what its corners score, and the tet mesher's report reads the same numbers. +#[test] +fn dihedral_angles_and_scaled_jacobian_on_the_same_tetrahedra() { + let both = |p: &[[f64; 3]]| { + let q = element_quality(ElementKind::Tet4, p).unwrap(); + (metric(&q, Metric::MinDihedral), metric(&q, Metric::MaxDihedral), metric(&q, Metric::ScaledJacobian)) + }; + let (lo, hi, sj) = both(®ular_tet()); + assert!((lo - 70.528_779).abs() < 1e-5 && (hi - lo).abs() < 1e-9 && close(sj, 1.0), "{lo} {hi} {sj}"); + + let flat = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.01]]; + let (lo, hi, sj) = both(&flat); + assert!(lo < 2.0 && hi > 178.0 && sj > 0.0 && sj < 0.02, "{lo} {hi} {sj}"); + + // Equilateral base of side 1, apex 12 above its centre. + let h = 3.0_f64.sqrt() / 2.0; + let needle = [[0., 0., 0.], [1., 0., 0.], [0.5, h, 0.], [0.5, h / 3.0, 12.0]]; + let (lo, hi, sj) = both(&needle); + assert!(lo > 59.0 && hi < 91.0, "{lo} {hi}"); + // Apex corner: 6V = (√3/2)·12, edges √(1/3 + 144) each. + let edge = (1.0 / 3.0 + 144.0_f64).sqrt(); + assert!(close(sj, 2.0_f64.sqrt() * h * 12.0 / edge.powi(3)), "{sj}"); + assert!(sj < 0.01); + + for p in [regular_tet(), flat, needle] { + let tet10 = element_quality(ElementKind::Tet10, &midpoints(&p, &jacobian::TET_EDGES)).unwrap(); + let (lo, hi, sj) = both(&p); + assert!(close(metric(&tet10, Metric::MinDihedral), lo) && close(metric(&tet10, Metric::MaxDihedral), hi)); + assert!(close(metric(&tet10, Metric::ScaledJacobian), sj)); + let r = crate::tet::tet_quality_stats(&p, &[[0, 1, 2, 3]]).unwrap(); + assert!(close(r.min_dihedral, lo) && close(r.max_dihedral, hi) && close(r.min_scaled_jacobian, sj), "{r:?}"); + } +} + #[test] fn right_corner_tetrahedron() { let q = element_quality(ElementKind::Tet4, &RIGHT_TET).unwrap(); diff --git a/crates/mesh/src/tet/mod.rs b/crates/mesh/src/tet/mod.rs index a877b6a..27e69d3 100644 --- a/crates/mesh/src/tet/mod.rs +++ b/crates/mesh/src/tet/mod.rs @@ -140,6 +140,10 @@ pub struct TetQuality { pub max_dihedral: f64, pub min_radius_ratio: f64, pub mean_radius_ratio: f64, + /// The smallest scaled Jacobian at the worst corner (as `mesh.quality` reports it). The + /// dihedral range does not bound it: a needle has every dihedral angle near 60° or 90° and a + /// nearly flat corner (#547). + pub min_scaled_jacobian: f64, /// Tetrahedra with a dihedral angle below 10° or above 170°. pub slivers: usize, /// Tetrahedra not admissible by FieldCraft's validity criterion. @@ -149,7 +153,7 @@ pub struct TetQuality { /// The quality of tetrahedra given by corner coordinates, with the metrics `mesh.quality` /// reports ([`crate::quality::element_quality`]). pub fn tet_quality_stats(points: &[[f64; 3]], tets: &[[u32; 4]]) -> Result { - let mut s = TetQuality { min_dihedral: 180.0, min_radius_ratio: 1.0, ..TetQuality::default() }; + let mut s = TetQuality { min_dihedral: 180.0, min_radius_ratio: 1.0, min_scaled_jacobian: 1.0, ..TetQuality::default() }; let mut sum = 0.0; for t in tets { let p = |v: u32| points.get(v as usize).copied().ok_or_else(|| Error::Failed("quality: vertex id".into())); @@ -159,6 +163,7 @@ pub fn tet_quality_stats(points: &[[f64; 3]], tets: &[[u32; 4]]) -> Result= 10.0 && hi <= 170.0)); s.inadmissible += usize::from(q.validity() != Some(Validity::Valid)); @@ -167,6 +172,7 @@ pub fn tet_quality_stats(points: &[[f64; 3]], tets: &[[u32; 4]]) -> Result 30.0 && f(&mesh["max_dihedral"]) < 145.0, "{order}: {mesh}"); + assert!(f(&mesh["min_scaled_jacobian"]) >= 0.3, "{order}: {mesh}"); + } + let _ = std::fs::remove_dir_all(&dir); +}