diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index 8f7b527..d86c49c 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -35,7 +35,8 @@ use fieldcraft_model::{ElementKind, admissible_jacobian, check_coordinates, tet4 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, + GateOptions, GateReport, Marked, MetricStats, QualityReport, enclosed_volume, mark_elements, quality_report, validity_gate, validity_gate_checks, + validity_gate_elements, }; pub(crate) use shape::{planar_corner, quad_axes}; @@ -77,6 +78,26 @@ pub enum Metric { } impl Metric { + pub const ALL: [Metric; 12] = [ + Metric::ScaledJacobian, + Metric::JacobianRatio, + Metric::Shape, + Metric::MinDihedral, + Metric::MaxDihedral, + Metric::RadiusRatio, + Metric::AspectRatio, + Metric::VolumeLength, + Metric::EdgeRatio, + Metric::MinAngle, + Metric::MaxAngle, + Metric::Skew, + ]; + + /// The metric named `name` ([`Metric::name`]). + pub fn parse(name: &str) -> Option { + Metric::ALL.into_iter().find(|m| m.name() == name) + } + pub fn name(self) -> &'static str { match self { Metric::ScaledJacobian => "scaled_jacobian", diff --git a/crates/mesh/src/quality/report.rs b/crates/mesh/src/quality/report.rs index 1800843..095964f 100644 --- a/crates/mesh/src/quality/report.rs +++ b/crates/mesh/src/quality/report.rs @@ -144,6 +144,71 @@ pub fn quality_report(model: &Model, labels: Option<&[u32]>, bins: usize) -> Res Ok(report) } +/// An element that fails a quality limit, or is not proven valid. +#[derive(Debug, Clone, PartialEq)] +pub struct Marked { + pub label: u32, + pub kind: ElementKind, + /// The mean of its nodes. + pub centroid: [f64; 3], + /// The certified validity, for solids. + pub validity: Option, + /// Each limit it fails: the metric, its value and the limit. + pub failed: Vec<(Metric, f64, f64)>, + /// Its scaled Jacobian, when it has one (the worst-first order). + pub scaled_jacobian: Option, +} + +impl Marked { + /// Not proven valid (inverted or undetermined). + pub fn is_invalid(&self) -> bool { + matches!(self.validity, Some(Validity::Invalid | Validity::Undetermined)) + } +} + +/// The elements (all, or the given labels) that fail a limit: a value below the limit of a +/// metric where higher is better ([`Metric::higher_is_better`]), above it otherwise, or no value +/// (NaN). Solids not proven valid are always marked. A metric an element does not have does not +/// apply to it. Returns how many elements were measured and the marked ones, worst first: not +/// valid, then by scaled Jacobian (lowest first; none counts as 1), then by label. +pub fn mark_elements(model: &Model, labels: Option<&[u32]>, limits: &[(Metric, f64)]) -> Result<(usize, Vec)> { + if let Some((m, v)) = limits.iter().find(|(_, v)| !v.is_finite()) { + return Err(Error::Invalid(format!("the limit of {} must be a finite number, got {v}", m.name()))); + } + let mut measured = 0; + let mut marked = Vec::new(); + for (label, element) in selection(model, labels)? { + if !is_measured(element.kind) { + continue; + } + let xyz = coordinates(model, label, element)?; + let q = element_quality(element.kind, &xyz)?; + measured += 1; + let failed: Vec<(Metric, f64, f64)> = limits + .iter() + .filter_map(|&(metric, limit)| { + let v = q.get(metric)?; + let fails = v.is_nan() || if metric.higher_is_better() { v < limit } else { v > limit }; + fails.then_some((metric, v, limit)) + }) + .collect(); + let validity = q.validity(); + let invalid = matches!(validity, Some(Validity::Invalid | Validity::Undetermined)); + if failed.is_empty() && !invalid { + continue; + } + let n = xyz.len().max(1) as f64; + let centroid = [0, 1, 2].map(|k| xyz.iter().map(|p| p[k]).sum::() / n); + marked.push(Marked { label, kind: element.kind, centroid, validity, failed, scaled_jacobian: q.get(Metric::ScaledJacobian) }); + } + let key = |m: &Marked| (!m.is_invalid(), m.scaled_jacobian.filter(|v| !v.is_nan()).unwrap_or(1.0)); + marked.sort_by(|a, b| { + let (ka, kb) = (key(a), key(b)); + ka.0.cmp(&kb.0).then(ka.1.total_cmp(&kb.1)).then(a.label.cmp(&b.label)) + }); + Ok((measured, marked)) +} + /// What the validity gate checks besides element validity. #[derive(Debug, Clone, Copy, PartialEq)] pub struct GateOptions { diff --git a/crates/mesh/src/quality/tests.rs b/crates/mesh/src/quality/tests.rs index dee4977..0bd50b8 100644 --- a/crates/mesh/src/quality/tests.rs +++ b/crates/mesh/src/quality/tests.rs @@ -309,3 +309,29 @@ fn report_and_gate() { assert!(validity_gate(&model, &GateOptions { expected_volume: Some(f64::NAN), ..Default::default() }).is_err()); assert!(enclosed_volume(&c, &[[0, 1, 99]]).is_err()); } + +/// A cube passes; a sheared hex fails a scaled-Jacobian floor above its value and passes one +/// below; a mirrored hex is marked whatever the limits; worst first. +#[test] +fn marks_elements_below_their_limits() { + let sheared: Vec<[f64; 3]> = UNIT_CUBE.iter().map(|c| [c[0] + 0.75 * c[2], c[1], c[2]]).collect(); + let mirrored: Vec<[f64; 3]> = UNIT_CUBE.iter().map(|c| [-c[0], c[1], c[2]]).collect(); + let m = model_of(&[(ElementKind::Hex8, UNIT_CUBE.to_vec()), (ElementKind::Hex8, sheared), (ElementKind::Hex8, mirrored)]); + let (measured, marked) = mark_elements(&m, None, &[(Metric::ScaledJacobian, 0.9)]).unwrap(); + assert_eq!(measured, 3); + assert_eq!(marked.iter().map(|k| k.label).collect::>(), [3, 2]); + assert!(marked[0].is_invalid()); + let (f, v, limit) = marked[1].failed[0]; + assert_eq!((f, limit), (Metric::ScaledJacobian, 0.9)); + assert!(close(v, 0.8)); + assert!(close(marked[1].centroid[0], 0.5 + 0.375) && close(marked[1].centroid[2], 0.5)); + let (_, loose) = mark_elements(&m, None, &[(Metric::ScaledJacobian, 0.5), (Metric::EdgeRatio, 10.0)]).unwrap(); + assert_eq!(loose.iter().map(|k| k.label).collect::>(), [3]); + // An edge ratio ceiling below the sheared hex's (its slanted edges are 1.25 long). + let (_, edges) = mark_elements(&m, Some(&[2]), &[(Metric::EdgeRatio, 1.2)]).unwrap(); + assert_eq!(edges.len(), 1); + assert!(mark_elements(&m, None, &[(Metric::ScaledJacobian, f64::NAN)]).is_err()); + assert!(mark_elements(&m, Some(&[9]), &[]).is_err()); + assert_eq!(Metric::parse("min_dihedral"), Some(Metric::MinDihedral)); + assert!(Metric::ALL.iter().all(|m| Metric::parse(m.name()) == Some(*m))); +}