diff --git a/Cargo.lock b/Cargo.lock index 4ceb38a..eee9dea 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -455,6 +455,7 @@ version = "0.1.0" dependencies = [ "log", "modelcraft-geometry", + "modelcraft-mesh", "modelcraft-model", "serde", "serde_json", diff --git a/ROADMAP.md b/ROADMAP.md index 5215267..1263830 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -6,7 +6,7 @@ LS-PrePost, with hexahedral meshing as the centre of gravity, plus its own linea | Milestone | Scope | Status | |---|---|---| | M0 Foundation | Workspace, quality gates (`cargo xtask ci`), FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft's kernel and tessellator | done | -| M1 Surface preparation | Geometry in: open `.solidcraft` designs, rebuilt by SolidCraft, with persistent face names (done); face-tagged 3MF files (`geometry.import.mesh`, docs/formats/tagged-mesh.md, done); exact predicates (Shewchuk adaptive, symbolic perturbation, indirect LPI/TPI, done); weld, closed/manifold/orientation checks and feature edges on the tessellated boundary; mesh quality metrics | in progress | +| M1 Surface preparation | Geometry in: open `.solidcraft` designs, rebuilt by SolidCraft, with persistent face names (done); face-tagged 3MF files (`geometry.import.mesh`, docs/formats/tagged-mesh.md, done); exact predicates (Shewchuk adaptive, symbolic perturbation, indirect LPI/TPI, done); weld, closed/manifold/orientation checks and feature edges on the tessellated boundary; mesh quality metrics and the certified validity gate (`modelcraft-mesh::quality`, `mesh.quality`, done) | in progress | | M2 Tetrahedral meshing | Bowyer–Watson Delaunay, robust constrained Delaunay boundary recovery (Diazzi et al. 2023), Delaunay refinement to a size field, ODT smoothing, Klingner–Shewchuk sliver removal; linear and quadratic tets | planned | | M3 Hexahedral meshing (core) | Octree/grid dual all-hex (valid on any input), then hex-dominant by tet recombination, then polycube for blocky parts; see docs/meshing/literature.md | planned | | M4 Model definition | Materials, sections, sets and surfaces, boundary conditions, loads, steps, output requests | planned | diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index 54ab18d..0d1ba38 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -11,6 +11,7 @@ workspace = true [dependencies] modelcraft-model = { workspace = true } modelcraft-geometry = { workspace = true } +modelcraft-mesh = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index c05b07d..da60020 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -1,7 +1,8 @@ -use modelcraft_model::{ElementKind, MAX_LABEL}; +use modelcraft_mesh::quality::{self, QualityReport}; +use modelcraft_model::{ElementKind, MAX_ELEMENTS}; use serde_json::{Value, json}; -use super::{label, xyz}; +use super::{label, labels, xyz}; use crate::{CommandSpec, EngineError, Result, Session}; pub const NODE_ADD: CommandSpec = CommandSpec { @@ -24,6 +25,13 @@ pub const ELEMENT_ADD: CommandSpec = CommandSpec { pub const ELEMENT_REMOVE: CommandSpec = CommandSpec { id: "mesh.element.remove", label: "Remove Element", params: r#"{"label": 1}"#, run: element_remove }; +pub const QUALITY: CommandSpec = CommandSpec { + id: "mesh.quality", + label: "Mesh Quality", + params: r#"{"elements": [labels] (optional, default all), "bins": 10 (1..=1000)} — certified validity (det J > 0 everywhere), volume, and min/max/mean/worst/histogram per metric (angles in degrees; histograms over each metric's fixed range, out-of-range values in the end bins)"#, + run: quality_report, +}; + fn node_add(s: &mut Session, p: &Value) -> Result { let l = label(p, "label")?; s.model.set_node(l, xyz(p, "xyz")?)?; @@ -39,15 +47,10 @@ fn element_add(s: &mut Session, p: &Value) -> Result { let l = label(p, "label")?; let name = p.get("kind").and_then(Value::as_str).unwrap_or_default(); let kind = ElementKind::from_name(name).ok_or_else(|| EngineError::Params(format!("unknown element kind `{name}`")))?; - let raw = p.get("nodes").and_then(Value::as_array).ok_or_else(|| EngineError::Params("`nodes` must be a list of node labels".into()))?; - if raw.len() != kind.node_count() { - return Err(EngineError::Params(format!("{name} needs {} nodes, got {}", kind.node_count(), raw.len()))); + let nodes = labels(p, "nodes", ElementKind::Hex20.node_count())?; + if nodes.len() != kind.node_count() { + return Err(EngineError::Params(format!("{name} needs {} nodes, got {}", kind.node_count(), nodes.len()))); } - let nodes = raw - .iter() - .map(|v| v.as_u64().and_then(|n| u32::try_from(n).ok()).filter(|n| (1..=MAX_LABEL).contains(n))) - .collect::>>() - .ok_or_else(|| EngineError::Params(format!("`nodes` must be node labels in 1..={MAX_LABEL}")))?; s.model.set_element(l, kind, nodes)?; Ok(json!({"label": l})) } @@ -56,3 +59,96 @@ fn element_remove(s: &mut Session, p: &Value) -> Result { s.model.remove_element(label(p, "label")?)?; Ok(Value::Null) } + +/// Labels listed in a report: the first ones, and how many there are. +const LISTED: usize = 100; + +/// A JSON number, or a string for ±∞ (JSON has no infinity). +fn number(v: f64) -> Value { + if v.is_finite() { json!(v) } else { json!(if v > 0.0 { "inf" } else { "-inf" }) } +} + +fn listed(v: &[u32]) -> Value { + json!({"count": v.len(), "labels": v.iter().take(LISTED).collect::>()}) +} + +fn quality_report(s: &mut Session, p: &Value) -> Result { + let bins = match p.get("bins") { + None => 10, + Some(v) => v.as_u64().and_then(|b| usize::try_from(b).ok()).ok_or_else(|| EngineError::Params("`bins` must be a whole number".into()))?, + }; + let selection = match p.get("elements") { + None => None, + Some(_) => Some(labels(p, "elements", MAX_ELEMENTS)?), + }; + let r: QualityReport = quality::quality_report(&s.model, selection.as_deref(), bins)?; + let skipped: serde_json::Map = r.skipped.iter().map(|(k, n)| (k.name().to_string(), json!(n))).collect(); + let metrics: serde_json::Map = r + .metrics + .iter() + .map(|m| { + let (lo, hi) = m.metric.range(); + let stats = json!({ + "count": m.count, + "min": number(m.min), + "max": number(m.max), + "mean": m.mean.map(number), + "worst": m.worst, + "histogram": {"range": [lo, hi], "counts": m.histogram}, + }); + (m.metric.name().to_string(), stats) + }) + .collect(); + Ok(json!({ + "measured": r.measured, + "skipped": skipped, + "valid": r.valid, + "invalid": listed(&r.invalid), + "undetermined": listed(&r.undetermined), + "volume": r.volume, + "area": r.area, + "metrics": metrics, + })) +} + +#[cfg(test)] +mod tests { + use serde_json::json; + + use crate::Session; + + fn hex(s: &mut Session, label: u32, first_node: u32, mirror: bool) { + let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]; + for (i, c) in corners.iter().enumerate() { + let x = if mirror { -c[0] - 2 } else { c[0] }; + s.execute("mesh.node.add", &json!({"label": first_node + i as u32, "xyz": [x, c[1], c[2]]})).unwrap(); + } + let nodes: Vec = (first_node..first_node + 8).collect(); + s.execute("mesh.element.add", &json!({"label": label, "kind": "hex8", "nodes": nodes})).unwrap(); + } + + #[test] + fn reports_quality_and_inverted_elements() { + let mut s = Session::default(); + hex(&mut s, 1, 1, false); + hex(&mut s, 2, 11, true); + let journal = s.journal.len(); + let r = s.execute("mesh.quality", &json!({"bins": 2})).unwrap(); + assert_eq!(r["measured"], 2); + assert_eq!(r["valid"], 1); + assert_eq!(r["invalid"], json!({"count": 1, "labels": [2]})); + assert_eq!(r["volume"], 0.0); + let sj = &r["metrics"]["scaled_jacobian"]; + assert_eq!((sj["min"].clone(), sj["max"].clone(), sj["worst"].clone()), (json!(-1.0), json!(1.0), json!(2))); + assert_eq!(sj["histogram"]["counts"], json!([1, 1])); + assert_eq!(s.journal.len(), journal, "a report does not change the model"); + + let one = s.execute("mesh.quality", &json!({"elements": [1, 1]})).unwrap(); + assert_eq!((one["measured"].clone(), one["volume"].clone()), (json!(1), json!(1.0))); + assert!(s.execute("mesh.quality", &json!({"elements": [3]})).is_err()); + assert!(s.execute("mesh.quality", &json!({"elements": [0]})).is_err()); + assert!(s.execute("mesh.quality", &json!({"bins": 0})).is_err()); + assert!(s.execute("mesh.quality", &json!({"bins": 1e9})).is_err()); + assert!(s.execute("mesh.quality", &json!({"elements": "all"})).is_err()); + } +} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 6f4300f..06aa4f0 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -4,6 +4,7 @@ mod geometry; mod mesh; mod model; +use modelcraft_model::MAX_LABEL; use serde_json::Value; use crate::{CommandSpec, EngineError, Result}; @@ -18,6 +19,7 @@ pub static ALL: &[CommandSpec] = &[ mesh::NODE_REMOVE, mesh::ELEMENT_ADD, mesh::ELEMENT_REMOVE, + mesh::QUALITY, ]; /// A required label parameter (1..=u32::MAX; the model checks its own range). @@ -38,3 +40,16 @@ pub(crate) fn xyz(p: &Value, key: &str) -> Result<[f64; 3]> { } Ok(out) } + +/// A required list of labels (each in 1..=MAX_LABEL), at most `max_len` long. +pub(crate) fn labels(p: &Value, key: &str, max_len: usize) -> Result> { + let raw = p + .get(key) + .and_then(Value::as_array) + .filter(|a| a.len() <= max_len) + .ok_or_else(|| EngineError::Params(format!("`{key}` must be a list of at most {max_len} labels")))?; + raw.iter() + .map(|v| v.as_u64().and_then(|n| u32::try_from(n).ok()).filter(|n| (1..=MAX_LABEL).contains(n))) + .collect::>>() + .ok_or_else(|| EngineError::Params(format!("`{key}` must be labels in 1..={MAX_LABEL}"))) +} diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 9239335..230e05d 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -26,6 +26,8 @@ pub enum EngineError { Model(#[from] modelcraft_model::Error), #[error(transparent)] Geometry(#[from] modelcraft_geometry::Error), + #[error(transparent)] + Mesh(#[from] modelcraft_mesh::Error), /// A command panicked; the model is unchanged. #[error("internal error in {command}: {message}")] Internal { command: String, message: String }, @@ -173,6 +175,9 @@ mod tests { json!({"path": 7}), json!({"body": "Body1", "tolerance": 0.05}), json!({"body": null, "tolerance": "fine"}), + json!({"elements": [u32::MAX, -1, 1.5]}), + json!({"elements": {}, "bins": -5}), + json!({"elements": [1], "bins": 1e300}), ]; for spec in commands() { for p in &hostile { diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index d802175..0185565 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -7,6 +7,8 @@ #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] +pub mod quality; + #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { #[error("invalid input: {0}")] diff --git a/crates/mesh/src/quality/bernstein.rs b/crates/mesh/src/quality/bernstein.rs new file mode 100644 index 0000000..50dc4dc --- /dev/null +++ b/crates/mesh/src/quality/bernstein.rs @@ -0,0 +1,505 @@ +//! Bernstein–Bézier polynomials on the reference domains of solid elements. +//! +//! A reference domain is a product of simplices: three segments for a hexahedron, a triangle +//! times a segment for a wedge, one tetrahedron for a tetrahedron. On a simplex of dimension `d` +//! a polynomial of degree `n` is written in the Bernstein basis +//! `B_α(λ) = n! / (α_0! … α_d!) · λ_0^α_0 … λ_d^α_d` (`|α| = n`, `λ` barycentric coordinates); +//! on a product the basis is the product of each factor's basis. +//! +//! Three properties make the basis the right tool for certifying an element's Jacobian +//! (Johnen, Remacle, Geuzaine, *Geometrical validity of curvilinear finite elements*, +//! J. Comput. Phys. 233:359–372, 2013): +//! 1. The basis is a partition of unity with non-negative terms, so the polynomial lies between +//! its smallest and its largest coefficient everywhere on the domain. +//! 2. A coefficient whose multi-index sits at a vertex of the domain is the polynomial's value +//! there. +//! 3. Subdividing the domain (de Casteljau) gives the coefficients on each piece, and they converge +//! to the values as the pieces shrink. +//! +//! So a positive smallest coefficient proves positivity, a non-positive vertex coefficient proves +//! the opposite, and subdivision settles the cases in between. The product rule +//! `B^m_α · B^n_β = C(m;α) C(n;β) / C(m+n;α+β) · B^{m+n}_{α+β}` (multinomials `C`) and the +//! derivative rule `∂p/∂λ_k − ∂p/∂λ_0 = n Σ_β (c_{β+e_k} − c_{β+e_0}) B^{n−1}_β` are the standard +//! ones (Farin, *Curves and Surfaces for CAGD*, 5th ed., ch. 17–18). + +use crate::{Error, Result}; + +/// Highest degree per factor. The Jacobian determinant of a 20-node hexahedron has degree 5. +pub(crate) const MAX_DEGREE: usize = 8; + +/// Free barycentric components of a multi-index: three for every 3-D reference domain (the +/// component `α_0` of each factor is implied by its degree). +const FREE: usize = 3; +type Multi = [usize; FREE]; + +/// One simplex factor of a reference domain (`dim` 1 is a segment). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) struct Factor { + pub dim: usize, + pub deg: usize, +} + +/// A polynomial in Bernstein form on a product of simplices. +/// +/// Coefficients are stored densely over `(deg + 1)^dim` slots per factor; slots whose free +/// components exceed the factor's degree are unused and stay zero. +#[derive(Debug, Clone, PartialEq)] +pub(crate) struct Bernstein { + factors: Vec, + coef: Vec, +} + +fn factorial(n: usize) -> f64 { + (2..=n).map(|k| k as f64).product() +} + +impl Bernstein { + pub fn zero(factors: &[Factor]) -> Result { + let free: usize = factors.iter().map(|f| f.dim).sum(); + if free > FREE || factors.iter().any(|f| f.dim == 0 || f.deg > MAX_DEGREE) { + return Err(Error::Invalid(format!("unsupported Bernstein domain {factors:?}"))); + } + let len = factors.iter().map(|f| (f.deg + 1).pow(f.dim as u32)).product(); + Ok(Self { factors: factors.to_vec(), coef: vec![0.0; len] }) + } + + pub fn factors(&self) -> &[Factor] { + &self.factors + } + + /// Position of factor `f`'s first free component in a multi-index. + fn offset(&self, f: usize) -> usize { + self.factors.iter().take(f).map(|x| x.dim).sum() + } + + fn slot(&self, m: &Multi) -> usize { + let mut s = 0; + let mut p = 0; + for f in &self.factors { + for _ in 0..f.dim { + s = s * (f.deg + 1) + m.get(p).copied().unwrap_or(0); + p += 1; + } + } + s + } + + fn valid(&self, m: &Multi) -> bool { + let mut p = 0; + self.factors.iter().all(|f| { + let sum: usize = m.iter().skip(p).take(f.dim).sum(); + p += f.dim; + sum <= f.deg + }) + } + + /// Every valid multi-index, in storage order. + fn multis(&self) -> Vec { + let mut out = vec![[0; FREE]]; + let mut p = 0; + for f in &self.factors { + for _ in 0..f.dim { + out = out + .into_iter() + .flat_map(|m| { + (0..=f.deg).map(move |c| { + let mut m = m; + if let Some(x) = m.get_mut(p) { + *x = c; + } + m + }) + }) + .collect(); + p += 1; + } + } + out.retain(|m| self.valid(m)); + out + } + + fn get(&self, m: &Multi) -> f64 { + self.coef.get(self.slot(m)).copied().unwrap_or(0.0) + } + + fn add_at(&mut self, m: &Multi, v: f64) { + let s = self.slot(m); + if let Some(c) = self.coef.get_mut(s) { + *c += v; + } + } + + /// Set the coefficient of `m` (free components, factor after factor). + pub fn set(&mut self, m: Multi, v: f64) { + let s = self.slot(&m); + if self.valid(&m) + && let Some(c) = self.coef.get_mut(s) + { + *c = v; + } + } + + /// Barycentric component `c` (0..=dim) of factor `f` in `m`. + fn component(&self, m: &Multi, f: usize, c: usize) -> usize { + let (Some(fac), off) = (self.factors.get(f), self.offset(f)) else { return 0 }; + if c == 0 { fac.deg.saturating_sub(m.iter().skip(off).take(fac.dim).sum()) } else { m.get(off + c - 1).copied().unwrap_or(0) } + } + + /// The multinomial weight `n! / α!` of `m`, over all factors. + fn weight(&self, m: &Multi) -> f64 { + let mut w = 1.0; + for (f, fac) in self.factors.iter().enumerate() { + w *= factorial(fac.deg); + for c in 0..=fac.dim { + w /= factorial(self.component(m, f, c)); + } + } + w + } + + fn same_domain(&self, o: &Self) -> bool { + self.factors.len() == o.factors.len() && self.factors.iter().zip(&o.factors).all(|(a, b)| a.dim == b.dim) + } + + /// The multinomial weight of every slot (zero for unused slots). + fn weights(&self) -> Vec { + let mut w = vec![0.0; self.coef.len()]; + for m in self.multis() { + let s = self.slot(&m); + if let Some(x) = w.get_mut(s) { + *x = self.weight(&m); + } + } + w + } + + /// The product of two polynomials on the same domain. + pub fn mul(&self, o: &Self) -> Result { + if !self.same_domain(o) { + return Err(Error::Invalid("Bernstein product of different domains".into())); + } + let factors: Vec = self.factors.iter().zip(&o.factors).map(|(a, b)| Factor { dim: a.dim, deg: a.deg + b.deg }).collect(); + let mut out = Self::zero(&factors)?; + let out_w = out.weights(); + let terms = + |p: &Self| -> Vec<(Multi, f64)> { p.multis().into_iter().map(|m| (m, p.get(&m) * p.weight(&m))).filter(|(_, c)| *c != 0.0).collect() }; + let right = terms(o); + for (a, ca) in terms(self) { + for (b, cb) in &right { + let mut k = a; + for (x, y) in k.iter_mut().zip(b) { + *x += y; + } + let s = out.slot(&k); + if let (Some(c), Some(w)) = (out.coef.get_mut(s), out_w.get(s)) { + *c += ca * cb / w; + } + } + } + Ok(out) + } + + /// `self + s · o`, both of the same domain and degrees. + pub fn axpy(&self, s: f64, o: &Self) -> Result { + if self.factors != o.factors { + return Err(Error::Invalid("Bernstein sum of different spaces".into())); + } + let coef = self.coef.iter().zip(&o.coef).map(|(a, b)| a + s * b).collect(); + Ok(Self { factors: self.factors.clone(), coef }) + } + + /// The derivative along reference direction `k` (1..=dim) of factor `f`: the direction from + /// that factor's vertex 0 to its vertex `k`, with the other coordinates held fixed. + pub fn derivative(&self, f: usize, k: usize) -> Result { + let fac = self.factors.get(f).copied().ok_or_else(|| Error::Invalid(format!("no factor {f}")))?; + if k == 0 || k > fac.dim || fac.deg == 0 { + return Err(Error::Invalid(format!("no derivative {k} of factor {fac:?}"))); + } + let mut factors = self.factors.clone(); + if let Some(x) = factors.get_mut(f) { + x.deg -= 1; + } + let mut out = Self::zero(&factors)?; + let pos = self.offset(f) + k - 1; + for b in out.multis() { + let mut hi = b; + if let Some(x) = hi.get_mut(pos) { + *x += 1; + } + let v = fac.deg as f64 * (self.get(&hi) - self.get(&b)); + out.add_at(&b, v); + } + Ok(out) + } + + /// Split factor `f` at the midpoint of its edge `(i, j)` (vertex numbers 0..=dim) with de + /// Casteljau's algorithm. The first piece keeps vertex `i` (vertex `j` becomes the midpoint), + /// the second keeps vertex `j` (vertex `i` becomes the midpoint). + pub fn bisect(&self, f: usize, i: usize, j: usize) -> Result<(Self, Self)> { + let fac = self.factors.get(f).copied().ok_or_else(|| Error::Invalid(format!("no factor {f}")))?; + if i == j || i > fac.dim || j > fac.dim { + return Err(Error::Invalid(format!("no edge ({i}, {j}) in {fac:?}"))); + } + let off = self.offset(f); + let shift = |start: &Multi, r: usize| { + let mut m = *start; + if i > 0 + && let Some(x) = m.get_mut(off + i - 1) + { + *x -= r; + } + if j > 0 + && let Some(x) = m.get_mut(off + j - 1) + { + *x += r; + } + m + }; + let (mut a, mut b) = (self.clone(), self.clone()); + // Each line starts where component j is zero; walking it moves units from i to j. + for start in self.multis() { + if self.component(&start, f, j) != 0 { + continue; + } + let s = self.component(&start, f, i); + let line: Vec = (0..=s).map(|r| shift(&start, r)).collect(); + let mut tri: Vec = line.iter().map(|m| self.get(m)).collect(); + let n = tri.len(); + let mut left = vec![0.0; n]; + let mut right = vec![0.0; n]; + for level in 0..n { + let last = n - 1 - level; + if let (Some(l), Some(&first)) = (left.get_mut(level), tri.first()) { + *l = first; + } + if let (Some(r), Some(&end)) = (right.get_mut(last), tri.get(last)) { + *r = end; + } + for r in 0..last { + if let (Some(&x), Some(&y)) = (tri.get(r), tri.get(r + 1)) + && let Some(t) = tri.get_mut(r) + { + *t = 0.5 * (x + y); + } + } + } + for ((m, l), r) in line.iter().zip(&left).zip(&right) { + let s = self.slot(m); + if let Some(c) = a.coef.get_mut(s) { + *c = *l; + } + if let Some(c) = b.coef.get_mut(s) { + *c = *r; + } + } + } + Ok((a, b)) + } + + /// The smallest and largest coefficient: bounds of the polynomial on its domain. + pub fn bounds(&self) -> (f64, f64) { + self.multis().iter().map(|m| self.get(m)).fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), c| (lo.min(c), hi.max(c))) + } + + /// The value at a vertex of the domain: one vertex number (0..=dim) per factor. + pub fn at_vertex(&self, verts: &[usize]) -> f64 { + let mut m = [0; FREE]; + for (f, (fac, &v)) in self.factors.iter().zip(verts).enumerate() { + if v > 0 + && v <= fac.dim + && let Some(x) = m.get_mut(self.offset(f) + v - 1) + { + *x = fac.deg; + } + } + self.get(&m) + } + + /// Every vertex of the domain, as one vertex number per factor. + pub fn vertices(&self) -> Vec> { + self.factors.iter().fold(vec![Vec::new()], |acc, f| { + acc.into_iter() + .flat_map(|v| { + (0..=f.dim).map(move |k| { + let mut v = v.clone(); + v.push(k); + v + }) + }) + .collect() + }) + } + + /// The values at every vertex of the domain. + pub fn vertex_values(&self) -> Vec { + self.vertices().iter().map(|v| self.at_vertex(v)).collect() + } + + /// The integral over the reference domain (unit segment, unit right triangle and unit right + /// tetrahedron factors). Every basis function of a factor has the same integral, + /// `1 / (dim! · C(deg + dim, dim))`. + pub fn integral(&self) -> f64 { + let per: f64 = self.factors.iter().map(|f| factorial(f.deg) / factorial(f.deg + f.dim)).product(); + self.multis().iter().map(|m| self.get(m)).sum::() * per + } + + /// The value at a point given by barycentric coordinates per factor (`dim + 1` each). + pub fn eval(&self, bary: &[&[f64]]) -> f64 { + self.multis() + .iter() + .map(|m| { + let mut t = self.get(m) * self.weight(m); + for (f, fac) in self.factors.iter().enumerate() { + for c in 0..=fac.dim { + let l = bary.get(f).and_then(|b| b.get(c)).copied().unwrap_or(0.0); + t *= l.powi(self.component(m, f, c) as i32); + } + } + t + }) + .sum() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn seg(deg: usize) -> Factor { + Factor { dim: 1, deg } + } + + /// A pseudo-random polynomial (deterministic). + fn sample(factors: &[Factor], seed: u64) -> Bernstein { + let mut p = Bernstein::zero(factors).unwrap(); + let mut s = seed; + for m in p.multis() { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + p.set(m, ((s >> 33) as f64 / (1u64 << 31) as f64) - 0.5); + } + p + } + + fn points(dim: usize) -> Vec> { + match dim { + 1 => vec![vec![0.7, 0.3], vec![0.1, 0.9], vec![1.0, 0.0]], + 2 => vec![vec![0.2, 0.3, 0.5], vec![0.6, 0.1, 0.3], vec![0.0, 0.0, 1.0]], + _ => vec![vec![0.1, 0.2, 0.3, 0.4], vec![0.25; 4], vec![0.0, 1.0, 0.0, 0.0]], + } + } + + fn domains() -> Vec> { + vec![vec![seg(2), seg(1), seg(3)], vec![Factor { dim: 2, deg: 2 }, seg(1)], vec![Factor { dim: 3, deg: 2 }]] + } + + /// Every combination of sample points, one per factor. + fn grid(factors: &[Factor]) -> Vec>> { + factors.iter().fold(vec![Vec::new()], |acc, f| { + acc.into_iter().flat_map(|v| points(f.dim).into_iter().map(move |p| [v.clone(), vec![p]].concat())).collect() + }) + } + + fn eval(p: &Bernstein, x: &[Vec]) -> f64 { + let refs: Vec<&[f64]> = x.iter().map(|v| v.as_slice()).collect(); + p.eval(&refs) + } + + #[test] + fn product_matches_pointwise_product() { + for d in domains() { + let a = sample(&d, 1); + let b = sample(&d.iter().map(|f| Factor { dim: f.dim, deg: 1 }).collect::>(), 2); + let ab = a.mul(&b).unwrap(); + for x in grid(&d) { + assert!((eval(&ab, &x) - eval(&a, &x) * eval(&b, &x)).abs() < 1e-12); + } + } + } + + #[test] + fn derivative_matches_finite_difference() { + for d in domains() { + let p = sample(&d, 3); + for (f, fac) in d.iter().enumerate() { + for k in 1..=fac.dim { + let dp = p.derivative(f, k).unwrap(); + for x in grid(&d) { + // Move along direction k: λ_k += h, λ_0 -= h. + let h = 1e-6; + let mut plus = x.clone(); + let mut minus = x.clone(); + plus[f][k] += h; + plus[f][0] -= h; + minus[f][k] -= h; + minus[f][0] += h; + let fd = (eval(&p, &plus) - eval(&p, &minus)) / (2.0 * h); + assert!((eval(&dp, &x) - fd).abs() < 1e-6, "{fd} vs {}", eval(&dp, &x)); + } + } + } + } + } + + #[test] + fn bisection_reproduces_the_polynomial() { + for d in domains() { + let p = sample(&d, 4); + for (f, fac) in d.iter().enumerate() { + for (i, j) in [(0, 1), (1, 0), (fac.dim, 0)] { + if i == j { + continue; + } + let (a, b) = p.bisect(f, i, j).unwrap(); + for x in grid(&d) { + // Piece a: vertex j is the midpoint, so λ_i += λ_j / 2, λ_j /= 2. + let mut xa = x.clone(); + let lj = xa[f][j]; + xa[f][i] += lj / 2.0; + xa[f][j] = lj / 2.0; + assert!((eval(&a, &x) - eval(&p, &xa)).abs() < 1e-12); + let mut xb = x.clone(); + let li = xb[f][i]; + xb[f][j] += li / 2.0; + xb[f][i] = li / 2.0; + assert!((eval(&b, &x) - eval(&p, &xb)).abs() < 1e-12); + } + } + } + } + } + + #[test] + fn vertices_bounds_and_integral() { + for d in domains() { + let p = sample(&d, 5); + let (lo, hi) = p.bounds(); + for v in p.vertices() { + let x: Vec> = d.iter().zip(&v).map(|(f, &k)| (0..=f.dim).map(|c| if c == k { 1.0 } else { 0.0 }).collect()).collect(); + let value = eval(&p, &x); + assert!((p.at_vertex(&v) - value).abs() < 1e-12); + assert!(lo <= value && value <= hi); + } + } + // ∫ over the unit tetrahedron of 1 is 1/6; over the triangle × segment of 1 is 1/2. + let mut one = Bernstein::zero(&[Factor { dim: 3, deg: 2 }]).unwrap(); + for m in one.multis() { + one.set(m, 1.0); + } + assert!((one.integral() - 1.0 / 6.0).abs() < 1e-15); + let mut one = Bernstein::zero(&[Factor { dim: 2, deg: 1 }, seg(2)]).unwrap(); + for m in one.multis() { + one.set(m, 1.0); + } + assert!((one.integral() - 0.5).abs() < 1e-15); + } + + #[test] + fn rejects_unsupported_domains() { + assert!(Bernstein::zero(&[Factor { dim: 3, deg: 1 }, seg(1)]).is_err()); + assert!(Bernstein::zero(&[seg(MAX_DEGREE + 1)]).is_err()); + let p = Bernstein::zero(&[Factor { dim: 3, deg: 1 }]).unwrap(); + assert!(p.derivative(0, 4).is_err()); + assert!(p.bisect(0, 2, 2).is_err()); + assert!(p.mul(&Bernstein::zero(&[seg(1), seg(1), seg(1)]).unwrap()).is_err()); + } +} diff --git a/crates/mesh/src/quality/jacobian.rs b/crates/mesh/src/quality/jacobian.rs new file mode 100644 index 0000000..d249b06 --- /dev/null +++ b/crates/mesh/src/quality/jacobian.rs @@ -0,0 +1,325 @@ +//! The Jacobian determinant of a solid element as an exact polynomial, its certified sign and +//! the element's exact volume. +//! +//! The geometry map of every solid element ModelCraft knows is a polynomial on its reference +//! domain, so its Jacobian determinant is one too: constant for tet4, cubic for tet10, degree +//! (2, 2) on triangle × segment for wedge6, degree (2, 2, 2) for hex8 and (5, 5, 5) for hex20. It +//! is built here in Bernstein form ([`super::bernstein`]) and its sign is settled by bounds and +//! subdivision, the method of Johnen, Remacle and Geuzaine (J. Comput. Phys. 233, 2013). +//! +//! Positive corner Jacobians do not make a trilinear hexahedron valid: its determinant has +//! degree 2 in each variable and can be negative inside while positive at all eight corners +//! (Marschner, Palmer, Zhang, Solomon, *Hexahedral mesh repair via sum-of-squares relaxation*, +//! CGF 39(5), 2020; Knupp, *Achieving finite element mesh quality via optimization of the +//! Jacobian matrix norm and associated quantities*, IJNME 48, 2000). Sampling at Gauss points is +//! no proof either. The check here is a proof: `Valid` means `det J > 0` on the whole element. +//! +//! Node orders (corners first, then mid-side nodes): +//! - tet4/tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; positive when node 4 +//! lies on the side of face 1-2-3 that its right-hand normal points to. +//! - wedge6: triangle 1-2-3, then 4-5-6 above it; positive when 4-5-6 lies on the side of 1-2-3's +//! right-hand normal. +//! - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; mid-side nodes 9–12 on the edges of 1-2-3-4, +//! 13–16 on those of 5-6-7-8, 17–20 on edges 1-5, 2-6, 3-7, 4-8. + +use modelcraft_model::ElementKind; + +use super::bernstein::{Bernstein, Factor}; +use crate::{Error, Result}; + +/// Reference coordinates of the hex corners on the unit cube, in node order. +pub(crate) const HEX_CORNERS: [[usize; 3]; 8] = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]; +/// Hex edges in the order of hex20's mid-side nodes. +pub(crate) const HEX_EDGES: [(usize, usize); 12] = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)]; +/// Tet edges in the order of tet10's mid-side nodes. +pub(crate) const TET_EDGES: [(usize, usize); 6] = [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)]; + +/// A Jacobian determinant is positive when it exceeds this fraction of its largest magnitude. +/// Anything smaller is degenerate to rounding and counts as invalid. +pub const RELATIVE_TOLERANCE: f64 = 1e-12; +/// Subdivision caps for the certificate: depth of one piece and pieces per element. +pub const MAX_DEPTH: u32 = 30; +pub const MAX_PIECES: usize = 4096; + +/// Whether an element's Jacobian determinant is positive everywhere. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Validity { + /// Proven positive on the whole element. + Valid, + /// Zero or negative somewhere (degenerate or inverted): a point where it is was found. + Invalid, + /// The subdivision caps were reached without a proof either way. Treat as invalid. + Undetermined, +} + +/// The certified Jacobian check of one solid element. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct JacobianCheck { + pub validity: Validity, + /// A lower bound of `min det J` divided by an upper bound of `max |det J|`, in `[-1, 1]`; + /// 1 for an affine element (straight-sided tet, parallelepiped). For a valid element it comes + /// from the subdivision that proved it, so it is conservative. + pub ratio: f64, + /// The exact signed volume (`∫ det J` over the reference domain). + pub volume: f64, +} + +/// The geometry map's Jacobian of one element: its columns (each three coordinate polynomials) +/// and its determinant. +pub(crate) struct Jacobian { + pub columns: [[Bernstein; 3]; 3], + pub det: Bernstein, +} + +fn domain(kind: ElementKind) -> Result<(Vec, [(usize, usize); 3])> { + let seg = |deg| Factor { dim: 1, deg }; + Ok(match kind { + ElementKind::Tet4 => (vec![Factor { dim: 3, deg: 1 }], [(0, 1), (0, 2), (0, 3)]), + ElementKind::Tet10 => (vec![Factor { dim: 3, deg: 2 }], [(0, 1), (0, 2), (0, 3)]), + ElementKind::Wedge6 => (vec![Factor { dim: 2, deg: 1 }, seg(1)], [(0, 1), (0, 2), (1, 1)]), + ElementKind::Hex8 => (vec![seg(1); 3], [(0, 1), (1, 1), (2, 1)]), + ElementKind::Hex20 => (vec![seg(2); 3], [(0, 1), (1, 1), (2, 1)]), + other => return Err(Error::Invalid(format!("{} is not a solid element", other.name()))), + }) +} + +/// Free components of a simplex vertex `v` scaled by `deg` (vertex 0 is all zeros). +fn simplex_vertex(v: usize, deg: usize) -> [usize; 3] { + let mut m = [0; 3]; + if v > 0 + && let Some(x) = m.get_mut(v - 1) + { + *x = deg; + } + m +} + +fn mid(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1]), 0.5 * (a[2] + b[2])] +} + +/// The quadratic Bernstein control point between two end points for a curve through `m` at the +/// parameter midpoint: `2m − (a + b)/2`. +fn quadratic_control(a: [f64; 3], m: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [2.0 * m[0] - 0.5 * (a[0] + b[0]), 2.0 * m[1] - 0.5 * (a[1] + b[1]), 2.0 * m[2] - 0.5 * (a[2] + b[2])] +} + +/// The 20-node serendipity hexahedron's shape functions at reference point `x ∈ [-1, 1]³`. +fn serendipity(x: [f64; 3], xyz: &[[f64; 3]]) -> [f64; 3] { + let mut out = [0.0; 3]; + let mut add = |n: f64, p: Option<&[f64; 3]>| { + if let Some(p) = p { + for (o, c) in out.iter_mut().zip(p) { + *o += n * c; + } + } + }; + let corner = |k: usize| HEX_CORNERS.get(k).map(|c| c.map(|v| 2.0 * v as f64 - 1.0)).unwrap_or([0.0; 3]); + for k in 0..8 { + let c = corner(k); + let n = 0.125 * (1.0 + x[0] * c[0]) * (1.0 + x[1] * c[1]) * (1.0 + x[2] * c[2]) * (x[0] * c[0] + x[1] * c[1] + x[2] * c[2] - 2.0); + add(n, xyz.get(k)); + } + for (e, &(a, b)) in HEX_EDGES.iter().enumerate() { + let r = mid(corner(a), corner(b)); + let mut n = 0.25; + for (xi, ri) in x.iter().zip(r) { + n *= if ri == 0.0 { 1.0 - xi * xi } else { 1.0 + xi * ri }; + } + add(n, xyz.get(8 + e)); + } + out +} + +/// The geometry map in Bernstein form: one polynomial per coordinate. +fn geometry(kind: ElementKind, xyz: &[[f64; 3]]) -> Result<[Bernstein; 3]> { + let (factors, _) = domain(kind)?; + let zero = Bernstein::zero(&factors)?; + let mut g = [zero.clone(), zero.clone(), zero]; + let mut put = |m: [usize; 3], p: [f64; 3]| { + for (poly, c) in g.iter_mut().zip(p) { + poly.set(m, c); + } + }; + let node = |i: usize| xyz.get(i).copied().ok_or_else(|| Error::Invalid(format!("{} needs {} nodes", kind.name(), kind.node_count()))); + match kind { + ElementKind::Tet4 => { + for v in 0..4 { + put(simplex_vertex(v, 1), node(v)?); + } + } + ElementKind::Tet10 => { + for v in 0..4 { + put(simplex_vertex(v, 2), node(v)?); + } + for (e, &(a, b)) in TET_EDGES.iter().enumerate() { + let (ma, mb) = (simplex_vertex(a, 1), simplex_vertex(b, 1)); + put([ma[0] + mb[0], ma[1] + mb[1], ma[2] + mb[2]], quadratic_control(node(a)?, node(4 + e)?, node(b)?)); + } + } + ElementKind::Wedge6 => { + for (i, base) in [[0, 0], [1, 0], [0, 1]].into_iter().enumerate() { + put([base[0], base[1], 0], node(i)?); + put([base[0], base[1], 1], node(i + 3)?); + } + } + ElementKind::Hex8 => { + for (i, c) in HEX_CORNERS.iter().enumerate() { + put(*c, node(i)?); + } + } + ElementKind::Hex20 => { + if xyz.len() < 20 { + return Err(Error::Invalid("hex20 needs 20 nodes".into())); + } + // The serendipity map lies in the triquadratic space, so interpolating it at the 27 + // triquadratic nodes reproduces it; then convert each axis from Lagrange to Bernstein. + let mut v = [[[[0.0; 3]; 3]; 3]; 3]; + for (a, plane) in v.iter_mut().enumerate() { + for (b, row) in plane.iter_mut().enumerate() { + for (c, p) in row.iter_mut().enumerate() { + *p = serendipity([a as f64 - 1.0, b as f64 - 1.0, c as f64 - 1.0], xyz); + } + } + } + // One full pass per axis (the conversion is a tensor product). + for axis in 0..3 { + for i in 0..3 { + for j in 0..3 { + let at = |k: usize| match axis { + 0 => (k, i, j), + 1 => (i, k, j), + _ => (i, j, k), + }; + let r = [0, 1, 2].map(|k| { + let (a, b, c) = at(k); + v[a][b][c] + }); + let (a, b, c) = at(1); + v[a][b][c] = quadratic_control(r[0], r[1], r[2]); + } + } + } + for (a, plane) in v.iter().enumerate() { + for (b, row) in plane.iter().enumerate() { + for (c, p) in row.iter().enumerate() { + put([a, b, c], *p); + } + } + } + } + other => return Err(Error::Invalid(format!("{} is not a solid element", other.name()))), + } + Ok(g) +} + +/// The Jacobian of a solid element whose node coordinates are already checked. +pub(crate) fn jacobian(kind: ElementKind, xyz: &[[f64; 3]]) -> Result { + let (_, dirs) = domain(kind)?; + let g = geometry(kind, xyz)?; + let column = |(f, k): (usize, usize)| -> Result<[Bernstein; 3]> { Ok([g[0].derivative(f, k)?, g[1].derivative(f, k)?, g[2].derivative(f, k)?]) }; + let columns = [column(dirs[0])?, column(dirs[1])?, column(dirs[2])?]; + let [a, b, c] = &columns; + // det [a b c] = a · (b × c) + let cross = |i: usize, j: usize| -> Result { b[i].mul(&c[j])?.axpy(-1.0, &b[j].mul(&c[i])?) }; + let det = a[0].mul(&cross(1, 2)?)?.axpy(1.0, &a[1].mul(&cross(2, 0)?)?)?.axpy(1.0, &a[2].mul(&cross(0, 1)?)?)?; + Ok(Jacobian { columns, det }) +} + +/// A piece of the reference domain during certification: the determinant on it and its vertices +/// in reference coordinates, per factor. +struct Piece { + det: Bernstein, + /// The determinant's smallest coefficient on this piece. + min: f64, + verts: Vec>, + depth: u32, +} + +impl Piece { + fn new(det: Bernstein, verts: Vec>, depth: u32) -> Self { + let min = det.bounds().0; + Piece { det, min, verts, depth } + } + + /// The longest edge over all factors (first on ties, so the choice is deterministic). + fn longest_edge(&self) -> (usize, usize, usize) { + let mut best = (0, 0, 1); + let mut len = -1.0; + for (f, vs) in self.verts.iter().enumerate() { + for (i, a) in vs.iter().enumerate() { + for (j, b) in vs.iter().enumerate().skip(i + 1) { + let l = (a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2); + if l > len { + len = l; + best = (f, i, j); + } + } + } + } + best + } + + fn split(&self) -> Result<(Piece, Piece)> { + let (f, i, j) = self.longest_edge(); + let (a, b) = self.det.bisect(f, i, j)?; + let vs = self.verts.get(f).ok_or_else(|| Error::Failed("certificate: lost a factor".into()))?; + let (pi, pj) = match (vs.get(i), vs.get(j)) { + (Some(p), Some(q)) => (*p, *q), + _ => return Err(Error::Failed("certificate: lost a vertex".into())), + }; + let m = mid(pi, pj); + let child = |det: Bernstein, replace: usize| { + let mut verts = self.verts.clone(); + if let Some(v) = verts.get_mut(f).and_then(|vs| vs.get_mut(replace)) { + *v = m; + } + Piece::new(det, verts, self.depth + 1) + }; + Ok((child(a, j), child(b, i))) + } +} + +/// Decide the sign of a Jacobian determinant on the whole reference domain. +pub(crate) fn certify(det: Bernstein) -> Result { + let (lo, hi) = det.bounds(); + let volume = det.integral(); + let scale = lo.abs().max(hi.abs()); + if !scale.is_finite() || !volume.is_finite() { + return Err(Error::Invalid("the element's Jacobian is not a finite number".into())); + } + let tol = RELATIVE_TOLERANCE * scale; + let ratio_of = |v: f64| if scale > 0.0 { v / scale } else { 0.0 }; + let verts = det.factors().iter().map(|f| (0..=f.dim).map(|k| simplex_vertex(k, 1).map(|x| x as f64)).collect()).collect(); + let mut stack = vec![Piece::new(det, verts, 0)]; + let mut lower = f64::INFINITY; + let mut undetermined = false; + let mut pieces = 0; + while let Some(piece) = stack.pop() { + pieces += 1; + if pieces > MAX_PIECES { + undetermined = true; + break; + } + if piece.det.vertex_values().iter().any(|v| !(*v > tol)) { + return Ok(JacobianCheck { validity: Validity::Invalid, ratio: ratio_of(lo), volume }); + } + if piece.min > tol { + lower = lower.min(piece.min); + continue; + } + if piece.depth >= MAX_DEPTH { + undetermined = true; + continue; + } + let (a, b) = piece.split()?; + // Look at the more doubtful piece first: a negative value is found sooner. + let (first, second) = if a.min <= b.min { (a, b) } else { (b, a) }; + stack.push(second); + stack.push(first); + } + if undetermined { + return Ok(JacobianCheck { validity: Validity::Undetermined, ratio: ratio_of(lo), volume }); + } + Ok(JacobianCheck { validity: Validity::Valid, ratio: ratio_of(lower), volume }) +} diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs new file mode 100644 index 0000000..e69751f --- /dev/null +++ b/crates/mesh/src/quality/mod.rs @@ -0,0 +1,228 @@ +//! Element quality: certified validity, exact volume and shape metrics, and the validity gate +//! every mesher runs before its mesh reaches the model. +//! +//! Sources (see `docs/meshing/literature.md`): +//! - Validity is a proof that `det J > 0` on the whole element, from the Jacobian determinant in +//! Bernstein form (Johnen, Remacle, Geuzaine, J. Comput. Phys. 233:359–372, 2013). Corner +//! Jacobians are not enough for a trilinear hexahedron (Marschner, Palmer, Zhang, Solomon, +//! CGF 39(5):133–147, 2020); see [`jacobian`]. +//! - Tetrahedron shape: Shewchuk, *What is a good linear finite element?*, 2002; Klingner and +//! Shewchuk, *Aggressive tetrahedral mesh improvement*, IMR 2007. +//! - Hexahedron, wedge and quadrilateral shape: Knupp, *Algebraic mesh quality metrics*, SIAM J. +//! Sci. Comput. 23(1):193–218, 2001. +//! +//! Inputs are hostile: coordinates must be finite and within [`modelcraft_model::MAX_COORD`], +//! and every metric has a defined value on degenerate and inverted elements. +//! +//! Linear tetrahedra take their validity from the exact `orient3d` ([`crate::predicates`]); the +//! other kinds from the floating-point Bernstein certificate, which counts a Jacobian within +//! [`jacobian::RELATIVE_TOLERANCE`] of zero as degenerate. + +mod bernstein; +pub mod jacobian; +mod report; +mod shape; + +use modelcraft_model::{ElementKind, check_coordinates}; + +pub use jacobian::{JacobianCheck, Validity}; +pub use report::{GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_elements}; + +use crate::predicates::{Point3, Sign, orient3d}; +use crate::{Error, Result}; + +/// A quality metric. Angles are in degrees. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum Metric { + /// Smallest normalized corner Jacobian, in `[-1, 1]`; 1 for a cube, a regular tetrahedron + /// and a square. Hexahedra also measure the centre. + ScaledJacobian, + /// Certified lower bound of `min det J` over an upper bound of `max |det J|`, in `[-1, 1]`; + /// 1 for affine elements. See [`JacobianCheck::ratio`]. + JacobianRatio, + /// Knupp's shape metric of the hexahedron corners (and centre), in `[0, 1]`; 1 for a cube. + Shape, + /// Smallest dihedral angle of a tetrahedron (70.53° for the regular one). + MinDihedral, + /// Largest dihedral angle of a tetrahedron. + MaxDihedral, + /// Three times the inradius over the circumradius, signed by orientation; 1 for the regular + /// tetrahedron. + RadiusRatio, + /// Longest edge over smallest altitude, normalized so the regular tetrahedron and the + /// equilateral triangle score 1; `+∞` when degenerate. + AspectRatio, + /// `6√2 V / ℓ_rms³`, signed; 1 for the regular tetrahedron. + VolumeLength, + /// Longest over shortest edge; `+∞` when an edge has zero length. + EdgeRatio, + /// Smallest corner angle of a triangle or quadrilateral. + MinAngle, + /// Largest corner angle of a triangle or quadrilateral. + MaxAngle, +} + +impl Metric { + pub fn name(self) -> &'static str { + match self { + Metric::ScaledJacobian => "scaled_jacobian", + Metric::JacobianRatio => "jacobian_ratio", + Metric::Shape => "shape", + Metric::MinDihedral => "min_dihedral", + Metric::MaxDihedral => "max_dihedral", + Metric::RadiusRatio => "radius_ratio", + Metric::AspectRatio => "aspect_ratio", + Metric::VolumeLength => "volume_length", + Metric::EdgeRatio => "edge_ratio", + Metric::MinAngle => "min_angle", + Metric::MaxAngle => "max_angle", + } + } + + /// Whether a larger value is a better element. + pub fn higher_is_better(self) -> bool { + !matches!(self, Metric::MaxDihedral | Metric::AspectRatio | Metric::EdgeRatio | Metric::MaxAngle) + } + + /// The fixed histogram range; values outside it fall in the first or last bin. + pub fn range(self) -> (f64, f64) { + match self { + Metric::ScaledJacobian | Metric::JacobianRatio | Metric::RadiusRatio | Metric::VolumeLength => (-1.0, 1.0), + Metric::Shape => (0.0, 1.0), + Metric::MinDihedral | Metric::MaxDihedral | Metric::MinAngle | Metric::MaxAngle => (0.0, 180.0), + Metric::AspectRatio | Metric::EdgeRatio => (1.0, 11.0), + } + } +} + +/// The quality of one element. +#[derive(Debug, Clone, PartialEq)] +pub struct ElementQuality { + pub kind: ElementKind, + /// Signed volume of a solid (exact for curved elements), area of a shell. + pub size: f64, + /// The certified Jacobian check, for solids. + pub jacobian: Option, + pub metrics: Vec<(Metric, f64)>, +} + +impl ElementQuality { + pub fn get(&self, metric: Metric) -> Option { + self.metrics.iter().find(|(m, _)| *m == metric).map(|(_, v)| *v) + } + + pub fn validity(&self) -> Option { + self.jacobian.map(|j| j.validity) + } +} + +/// Whether [`element_quality`] measures this kind (solids and shells; not rods or beams). +pub fn is_measured(kind: ElementKind) -> bool { + kind.dimension() >= 2 +} + +/// Check the node count and coordinates, and translate the element so its first node is at the +/// origin (less cancellation far from the origin). +fn local(kind: ElementKind, xyz: &[[f64; 3]]) -> Result> { + if xyz.len() != kind.node_count() { + return Err(Error::Invalid(format!("{} needs {} nodes, got {}", kind.name(), kind.node_count(), xyz.len()))); + } + check_coordinates(xyz.as_flattened())?; + let o = xyz.first().copied().unwrap_or([0.0; 3]); + Ok(xyz.iter().map(|p| shape::sub(*p, o)).collect()) +} + +/// The first `N` points (corner nodes). +fn corners(p: &[[f64; 3]]) -> Result<[[f64; 3]; N]> { + p.get(..N).and_then(|s| <[[f64; 3]; N]>::try_from(s).ok()).ok_or_else(|| Error::Invalid(format!("need {N} corner nodes"))) +} + +/// The Jacobian matrix's columns as polynomials (each three coordinates). +type Columns = [[bernstein::Bernstein; 3]; 3]; + +/// The certified check of a solid, and its Jacobian columns when they were built. A linear +/// tetrahedron's Jacobian is constant: its sign is the exact `orient3d` of the input +/// coordinates, its volume the determinant. +fn checked(kind: ElementKind, xyz: &[[f64; 3]], p: &[[f64; 3]]) -> Result<(Option, JacobianCheck)> { + if kind == ElementKind::Tet4 + && let Some(sign) = exact_tet_sign(xyz) + { + let [a, b, c, d]: [[f64; 3]; 4] = corners(p)?; + let volume = shape::det(shape::sub(b, a), shape::sub(c, a), shape::sub(d, a)) / 6.0; + let (validity, ratio) = match sign { + Sign::Positive => (Validity::Valid, 1.0), + Sign::Negative => (Validity::Invalid, -1.0), + Sign::Zero => (Validity::Invalid, 0.0), + }; + return Ok((None, JacobianCheck { validity, ratio, volume })); + } + let jacobian::Jacobian { columns, det } = jacobian::jacobian(kind, p)?; + Ok((Some(columns), jacobian::certify(det)?)) +} + +/// The exact orientation of a tetrahedron (positive for node 4 on the side of face 1-2-3's +/// right-hand normal), or `None` when a coordinate is outside the predicates' domain. +fn exact_tet_sign(xyz: &[[f64; 3]]) -> Option { + let p = |i: usize| xyz.get(i).and_then(|c| Point3::try_from(*c).ok()); + // orient3d(a, b, c, d) is positive when d lies below the plane of a, b, c seen + // counterclockwise, i.e. opposite to our convention. + Some(match orient3d(&p(0)?, &p(1)?, &p(2)?, &p(3)?) { + Sign::Positive => Sign::Negative, + Sign::Negative => Sign::Positive, + Sign::Zero => Sign::Zero, + }) +} + +/// The certified Jacobian check and exact volume of a solid element (tet4, tet10, wedge6, hex8, +/// hex20) from its node coordinates in node order. +pub fn element_validity(kind: ElementKind, xyz: &[[f64; 3]]) -> Result { + let p = local(kind, xyz)?; + checked(kind, xyz, &p).map(|(_, check)| check) +} + +/// Knupp's scaled Jacobian and shape of a hexahedron: the Jacobian matrix at the eight corners +/// and the centre (hex20 uses its curved map). +fn hex_metrics(columns: &Columns) -> (f64, f64) { + let at = |value: &dyn Fn(&bernstein::Bernstein) -> f64| { + let [a, b, c] = columns.each_ref().map(|col| col.each_ref().map(value)); + shape::jacobian_metrics(a, b, c) + }; + let half: &[f64] = &[0.5, 0.5]; + let centre = at(&|q| q.eval(&[half, half, half])); + jacobian::HEX_CORNERS.iter().map(|corner| at(&|q| q.at_vertex(corner))).fold(centre, |(s, k), (s2, k2)| (s.min(s2), k.min(k2))) +} + +/// Every metric of one element from its node coordinates in node order. Quadratic elements are +/// certified on their curved geometry; tet10 and wedge shape metrics are those of their corners. +/// Rods and beams are an error ([`is_measured`]). +pub fn element_quality(kind: ElementKind, xyz: &[[f64; 3]]) -> Result { + let p = local(kind, xyz)?; + let mut metrics = match kind { + ElementKind::Tri3 | ElementKind::Tri6 => { + let (area, metrics) = shape::tri(corners(&p)?); + return Ok(ElementQuality { kind, size: area, jacobian: None, metrics }); + } + ElementKind::Quad4 | ElementKind::Quad8 => { + let (area, metrics) = shape::quad(corners(&p)?); + return Ok(ElementQuality { kind, size: area, jacobian: None, metrics }); + } + ElementKind::Rod2 | ElementKind::Beam2 => return Err(Error::Invalid(format!("quality metrics are not defined for {}", kind.name()))), + ElementKind::Tet4 | ElementKind::Tet10 => shape::tet(corners(&p)?), + ElementKind::Wedge6 => { + let c: [[f64; 3]; 6] = corners(&p)?; + let edges = [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (0, 3), (1, 4), (2, 5)]; + vec![(Metric::ScaledJacobian, shape::wedge_scaled_jacobian(c)), (Metric::EdgeRatio, shape::edge_ratio(&c, &edges))] + } + ElementKind::Hex8 | ElementKind::Hex20 => vec![(Metric::EdgeRatio, shape::edge_ratio(&corners::<8>(&p)?, &jacobian::HEX_EDGES))], + }; + let (columns, check) = checked(kind, xyz, &p)?; + if let Some(columns) = columns.as_ref().filter(|_| matches!(kind, ElementKind::Hex8 | ElementKind::Hex20)) { + let (scaled, knupp) = hex_metrics(columns); + metrics.extend([(Metric::ScaledJacobian, scaled), (Metric::Shape, knupp)]); + } + metrics.push((Metric::JacobianRatio, check.ratio)); + Ok(ElementQuality { kind, size: check.volume, jacobian: Some(check), metrics }) +} + +#[cfg(test)] +mod tests; diff --git a/crates/mesh/src/quality/report.rs b/crates/mesh/src/quality/report.rs new file mode 100644 index 0000000..3835aa0 --- /dev/null +++ b/crates/mesh/src/quality/report.rs @@ -0,0 +1,271 @@ +//! Mesh-wide quality: statistics per metric, and the validity gate meshers run on their output. + +use std::collections::{BTreeMap, BTreeSet}; + +use modelcraft_model::{Element, ElementKind, Model, check_coordinates}; + +use super::{ElementQuality, Metric, Validity, element_quality, element_validity, is_measured}; +use crate::{Error, Result}; + +/// The most histogram bins a report may ask for. +pub const MAX_BINS: usize = 1000; + +/// Statistics of one metric over the measured elements. +#[derive(Debug, Clone, PartialEq)] +pub struct MetricStats { + pub metric: Metric, + /// Elements that have this metric. + pub count: usize, + pub min: f64, + pub max: f64, + /// Mean of the finite values (`None` when there are none). + pub mean: Option, + /// The element with the worst value (the first one on ties). + pub worst: u32, + /// Counts over [`Metric::range`] in equal bins; values outside the range fall in the end bins. + pub histogram: Vec, +} + +/// Quality of a set of elements. +#[derive(Debug, Clone, PartialEq)] +pub struct QualityReport { + /// Elements measured (solids and shells). + pub measured: usize, + /// Elements of kinds with no quality metrics (rods, beams), by kind. + pub skipped: Vec<(ElementKind, usize)>, + /// Solids proven valid. + pub valid: usize, + pub invalid: Vec, + pub undetermined: Vec, + /// Total signed volume of the solids and total area of the shells. + pub volume: f64, + pub area: f64, + /// One entry per metric present, in [`Metric`] declaration order. + pub metrics: Vec, +} + +fn coordinates(model: &Model, label: u32, element: &Element) -> Result> { + element + .nodes + .iter() + .map(|n| { + model.node(*n).map(|node| node.xyz).ok_or_else(|| Error::Model(modelcraft_model::Error::NotFound(format!("element {label}: node {n}")))) + }) + .collect() +} + +/// The elements to look at: all of them, or the given labels (sorted, without repeats). +fn selection<'a>(model: &'a Model, labels: Option<&[u32]>) -> Result> { + match labels { + None => Ok(model.elements().iter().map(|(l, e)| (*l, e)).collect()), + Some(labels) => labels + .iter() + .collect::>() + .into_iter() + .map(|l| model.element(*l).map(|e| (*l, e)).ok_or_else(|| Error::Model(modelcraft_model::Error::NotFound(format!("element {l}"))))) + .collect(), + } +} + +/// Quality statistics of the model's elements, or of the given ones, with `bins` histogram bins +/// (1..=[`MAX_BINS`]) per metric. Deterministic: elements are visited in label order. +pub fn quality_report(model: &Model, labels: Option<&[u32]>, bins: usize) -> Result { + if !(1..=MAX_BINS).contains(&bins) { + return Err(Error::Invalid(format!("bins must be in 1..={MAX_BINS}, got {bins}"))); + } + let mut report = QualityReport { + measured: 0, + skipped: Vec::new(), + valid: 0, + invalid: Vec::new(), + undetermined: Vec::new(), + volume: 0.0, + area: 0.0, + metrics: Vec::new(), + }; + let mut skipped: BTreeMap = BTreeMap::new(); + let mut stats: BTreeMap = BTreeMap::new(); + for (label, element) in selection(model, labels)? { + if !is_measured(element.kind) { + *skipped.entry(element.kind).or_insert(0) += 1; + continue; + } + let q: ElementQuality = element_quality(element.kind, &coordinates(model, label, element)?)?; + report.measured += 1; + match q.validity() { + Some(Validity::Valid) => report.valid += 1, + Some(Validity::Invalid) => report.invalid.push(label), + Some(Validity::Undetermined) => report.undetermined.push(label), + None => {} + } + if element.kind.dimension() == 3 { + report.volume += q.size; + } else { + report.area += q.size; + } + for (metric, v) in q.metrics { + let (s, sum, finite) = stats.entry(metric).or_insert_with(|| { + ( + MetricStats { metric, count: 0, min: f64::INFINITY, max: f64::NEG_INFINITY, mean: None, worst: label, histogram: vec![0; bins] }, + 0.0, + 0, + ) + }); + if v.is_nan() { + continue; + } + let worse = if metric.higher_is_better() { v < s.min } else { v > s.max }; + if worse || s.count == 0 { + s.worst = label; + } + s.count += 1; + s.min = s.min.min(v); + s.max = s.max.max(v); + if v.is_finite() { + *sum += v; + *finite += 1; + } + let (lo, hi) = metric.range(); + // `as` saturates: −∞ lands in bin 0, +∞ in the last bin. + let bin = (((v - lo) / (hi - lo)) * bins as f64).floor().max(0.0) as usize; + if let Some(b) = s.histogram.get_mut(bin.min(bins - 1)) { + *b += 1; + } + } + } + report.skipped = skipped.into_iter().collect(); + report.metrics = stats + .into_values() + .map(|(mut s, sum, finite)| { + s.mean = (finite > 0).then(|| sum / finite as f64); + s + }) + .collect(); + Ok(report) +} + +/// What the validity gate checks besides element validity. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct GateOptions { + /// The volume the mesh must fill: the boundary's enclosed volume ([`enclosed_volume`]). + pub expected_volume: Option, + /// Relative tolerance on the volume. Curved (quadratic) meshes of a faceted boundary need a + /// looser one than linear meshes. + pub volume_tolerance: f64, +} + +impl Default for GateOptions { + fn default() -> Self { + Self { expected_volume: None, volume_tolerance: 1e-9 } + } +} + +/// The result of the validity gate. +#[derive(Debug, Clone, PartialEq)] +pub struct GateReport { + /// Solid elements checked (other kinds are ignored). + pub checked: usize, + pub invalid: Vec, + pub undetermined: Vec, + /// Total signed volume of the solids. + pub volume: f64, + pub options: GateOptions, +} + +impl GateReport { + /// Whether the volume matches the expected one (true when none is expected). + pub fn volume_matches(&self) -> bool { + self.options.expected_volume.is_none_or(|v| (self.volume - v).abs() <= self.options.volume_tolerance * v.abs()) + } + + /// Every solid is proven valid and the volume matches. + pub fn passed(&self) -> bool { + self.invalid.is_empty() && self.undetermined.is_empty() && self.volume_matches() + } + + /// The report, or [`Error::Failed`] describing why the gate failed. + pub fn into_result(self) -> Result { + if self.passed() { + return Ok(self); + } + let mut why = Vec::new(); + let list = |v: &[u32]| v.iter().take(10).map(u32::to_string).collect::>().join(", "); + if !self.invalid.is_empty() { + why.push(format!("{} inverted or degenerate elements ({})", self.invalid.len(), list(&self.invalid))); + } + if !self.undetermined.is_empty() { + why.push(format!("{} elements whose validity could not be proven ({})", self.undetermined.len(), list(&self.undetermined))); + } + if let (false, Some(v)) = (self.volume_matches(), self.options.expected_volume) { + why.push(format!("mesh volume {} differs from the boundary's {v}", self.volume)); + } + Err(Error::Failed(format!("mesh validity gate: {}", why.join("; ")))) + } +} + +/// Accumulates the gate's report one element at a time. +struct Gate(GateReport); + +impl Gate { + fn new(options: &GateOptions) -> Result { + if !(options.volume_tolerance >= 0.0 && options.volume_tolerance.is_finite()) || options.expected_volume.is_some_and(|v| !v.is_finite()) { + return Err(Error::Invalid("gate options must be finite, with a non-negative tolerance".into())); + } + Ok(Gate(GateReport { checked: 0, invalid: Vec::new(), undetermined: Vec::new(), volume: 0.0, options: *options })) + } + + fn add(&mut self, label: u32, kind: ElementKind, xyz: &[[f64; 3]]) -> Result<()> { + if kind.dimension() != 3 { + return Ok(()); + } + let check = element_validity(kind, xyz)?; + let r = &mut self.0; + r.checked += 1; + r.volume += check.volume; + match check.validity { + Validity::Valid => {} + Validity::Invalid => r.invalid.push(label), + Validity::Undetermined => r.undetermined.push(label), + } + Ok(()) + } +} + +/// The validity gate every mesher runs on its output before it reaches the model: each solid +/// element's Jacobian determinant is proven positive everywhere, and the elements' total volume +/// equals the volume the boundary encloses. Elements come as `(label, kind, node coordinates in +/// node order)`, straight from the mesher's arrays; other than solids are ignored. Call +/// [`GateReport::into_result`] to turn a failure into an error. +pub fn validity_gate_elements(elements: I, options: &GateOptions) -> Result +where + I: IntoIterator, + E: AsRef<[[f64; 3]]>, +{ + let mut gate = Gate::new(options)?; + for (label, kind, xyz) in elements { + gate.add(label, kind, xyz.as_ref())?; + } + Ok(gate.0) +} + +/// [`validity_gate_elements`] over a model's elements. +pub fn validity_gate(model: &Model, options: &GateOptions) -> Result { + let mut gate = Gate::new(options)?; + for (label, element) in model.elements() { + gate.add(*label, element.kind, &coordinates(model, *label, element)?)?; + } + Ok(gate.0) +} + +/// The signed volume enclosed by a closed triangle surface (positive when the triangles' right- +/// hand normals point outwards), by the divergence theorem. +pub fn enclosed_volume(points: &[[f64; 3]], triangles: &[[usize; 3]]) -> Result { + check_coordinates(points.as_flattened())?; + let o = points.first().copied().unwrap_or([0.0; 3]); + let mut six_v = 0.0; + for t in triangles { + let p = |i: usize| points.get(i).map(|p| super::shape::sub(*p, o)).ok_or_else(|| Error::Invalid(format!("triangle {t:?}: no point {i}"))); + six_v += super::shape::det(p(t[0])?, p(t[1])?, p(t[2])?); + } + Ok(six_v / 6.0) +} diff --git a/crates/mesh/src/quality/shape.rs b/crates/mesh/src/quality/shape.rs new file mode 100644 index 0000000..ad610d7 --- /dev/null +++ b/crates/mesh/src/quality/shape.rs @@ -0,0 +1,182 @@ +//! Shape metrics from corner positions and corner Jacobians. +//! +//! - Tetrahedra: Shewchuk, *What is a good linear finite element? Interpolation, conditioning, +//! anisotropy, and quality measures*, 2002 (dihedral angles, radius ratio, aspect ratio), and +//! the volume–length measure Klingner and Shewchuk recommend for mesh improvement (IMR 2007, §2). +//! - Hexahedra, wedges and quadrilaterals: Knupp's algebraic metrics built on the corner Jacobian +//! matrices (Knupp, *Algebraic mesh quality metrics*, SIAM J. Sci. Comput. 23(1), 2001): the +//! scaled Jacobian `det(A) / (|a1| |a2| |a3|)` and the shape metric `3 det(A)^{2/3} / |A|²_F`. +//! +//! Every metric is defined for degenerate input (no division by zero, no NaN): a degenerate +//! element scores 0 on bounded metrics and `+∞` on aspect and edge ratios. + +use std::f64::consts::SQRT_2; + +use super::Metric; + +pub(crate) type P = [f64; 3]; + +pub(crate) fn sub(a: P, b: P) -> P { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} +pub(crate) fn dot(a: P, b: P) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} +pub(crate) fn cross(a: P, b: P) -> P { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} +pub(crate) fn norm(a: P) -> f64 { + dot(a, a).sqrt() +} +pub(crate) fn det(a: P, b: P, c: P) -> f64 { + dot(a, cross(b, c)) +} + +/// The angle between two vectors in degrees, 0 when either is zero. +fn angle(a: P, b: P) -> f64 { + norm(cross(a, b)).atan2(dot(a, b)).to_degrees() +} + +/// `num / den`, or `fallback` when the denominator is zero or the quotient is not finite. +fn ratio(num: f64, den: f64, fallback: f64) -> f64 { + let r = num / den; + if den != 0.0 && r.is_finite() { r } else { fallback } +} + +/// Normalized scaled Jacobian of one corner: the determinant of its three unit edge vectors. +fn corner_det(a: P, b: P, c: P) -> f64 { + ratio(det(a, b, c), norm(a) * norm(b) * norm(c), 0.0) +} + +/// Largest over smallest edge length. +pub(crate) fn edge_ratio(p: &[P], edges: &[(usize, usize)]) -> f64 { + let (mut lo, mut hi) = (f64::INFINITY, 0.0_f64); + for &(a, b) in edges { + if let (Some(&x), Some(&y)) = (p.get(a), p.get(b)) { + let l = norm(sub(x, y)); + lo = lo.min(l); + hi = hi.max(l); + } + } + ratio(hi, lo, f64::INFINITY).max(1.0) +} + +/// Tetrahedron metrics on corners `p` (positive orientation as in [`super::jacobian`]). +pub(crate) fn tet(p: [P; 4]) -> Vec<(Metric, f64)> { + let [a, b, c, d] = p; + let volume = det(sub(b, a), sub(c, a), sub(d, a)) / 6.0; + let edges = super::jacobian::TET_EDGES.map(|(i, j)| norm(sub(p[j], p[i]))); + let faces = + [cross(sub(c, b), sub(d, b)), cross(sub(d, a), sub(c, a)), cross(sub(b, a), sub(d, a)), cross(sub(c, a), sub(b, a))].map(|n| 0.5 * norm(n)); + let area: f64 = faces.iter().sum(); + let l_max = edges.iter().fold(0.0_f64, |m, &l| m.max(l)); + let a_max = faces.iter().fold(0.0_f64, |m, &l| m.max(l)); + let l_rms = (edges.iter().map(|l| l * l).sum::() / 6.0).sqrt(); + + // Dihedral angle at edge (i, j): the angle between the other two vertices as seen + // perpendicular to the edge. + let mut dihedral = (180.0_f64, 0.0_f64); + // The two vertices off each edge of TET_EDGES, in the same order. + let opposite = [(2, 3), (0, 3), (1, 3), (1, 2), (0, 2), (0, 1)]; + for ((i, j), (k, l)) in super::jacobian::TET_EDGES.into_iter().zip(opposite) { + let e = sub(p[j], p[i]); + let ee = dot(e, e); + let perp = |v: P| { + let w = sub(v, p[i]); + let t = ratio(dot(w, e), ee, 0.0); + sub(w, [t * e[0], t * e[1], t * e[2]]) + }; + let th = angle(perp(p[k]), perp(p[l])); + dihedral = (dihedral.0.min(th), dihedral.1.max(th)); + } + + // Circumradius R = |num| / (12 |V|) with num = |a|²(b×c) + |b|²(c×a) + |c|²(a×b) for edge + // vectors a, b, c from one vertex; inradius r = 3V / area; radius ratio 3r/R. + let (u, v, w) = (sub(b, a), sub(c, a), sub(d, a)); + let (bc, ca, ab) = (cross(v, w), cross(w, u), cross(u, v)); + let (lu, lv, lw) = (dot(u, u), dot(v, v), dot(w, w)); + let num = norm([lu * bc[0] + lv * ca[0] + lw * ab[0], lu * bc[1] + lv * ca[1] + lw * ab[1], lu * bc[2] + lv * ca[2] + lw * ab[2]]); + let radius_ratio = ratio(108.0 * volume * volume.abs(), area * num, 0.0); + + // Smallest altitude 3|V| / largest face; normalized so the regular tetrahedron scores 1. + let aspect = ratio((2.0_f64 / 3.0).sqrt() * l_max * a_max, 3.0 * volume.abs(), f64::INFINITY); + let volume_length = ratio(6.0 * SQRT_2 * volume, l_rms.powi(3), 0.0); + + // Scaled Jacobian: 6V over the largest product of the three edge lengths at a corner, × √2. + let corner = [(0, 1, 2, 3), (1, 2, 0, 3), (2, 0, 1, 3), (3, 0, 2, 1)] + .map(|(o, x, y, z)| norm(sub(p[x], p[o])) * norm(sub(p[y], p[o])) * norm(sub(p[z], p[o]))) + .iter() + .fold(0.0_f64, |m, &x| m.max(x)); + let scaled = ratio(SQRT_2 * 6.0 * volume, corner, 0.0); + + vec![ + (Metric::ScaledJacobian, scaled), + (Metric::MinDihedral, dihedral.0), + (Metric::MaxDihedral, dihedral.1), + (Metric::RadiusRatio, radius_ratio), + (Metric::AspectRatio, aspect.max(1.0)), + (Metric::VolumeLength, volume_length), + ] +} + +/// Triangle area and metrics. +pub(crate) fn tri(p: [P; 3]) -> (f64, Vec<(Metric, f64)>) { + let area = 0.5 * norm(cross(sub(p[1], p[0]), sub(p[2], p[0]))); + let angles = [0, 1, 2].map(|i| angle(sub(p[(i + 1) % 3], p[i]), sub(p[(i + 2) % 3], p[i]))); + let l_max = (0..3).map(|i| norm(sub(p[(i + 1) % 3], p[i]))).fold(0.0_f64, f64::max); + // Longest edge over smallest altitude, normalized so the equilateral triangle scores 1. + let aspect = ratio(3.0_f64.sqrt() * l_max * l_max, 4.0 * area, f64::INFINITY).max(1.0); + ( + area, + vec![ + (Metric::MinAngle, angles.iter().fold(180.0, |m: f64, &a| m.min(a))), + (Metric::MaxAngle, angles.iter().fold(0.0, |m: f64, &a| m.max(a))), + (Metric::AspectRatio, aspect), + ], + ) +} + +/// Quadrilateral area (of its projection on the mean plane) and metrics. The scaled Jacobian of +/// a corner is measured against the normal of the diagonals, so a concave or bow-tie +/// quadrilateral scores ≤ 0. +pub(crate) fn quad(p: [P; 4]) -> (f64, Vec<(Metric, f64)>) { + let n = cross(sub(p[2], p[0]), sub(p[3], p[1])); + let area = 0.5 * norm(n); + let mut scaled = f64::INFINITY; + let mut angles = (180.0_f64, 0.0_f64); + for i in 0..4 { + let next = sub(p[(i + 1) % 4], p[i]); + let prev = sub(p[(i + 3) % 4], p[i]); + scaled = scaled.min(ratio(dot(cross(next, prev), n), norm(next) * norm(prev) * norm(n), 0.0)); + let a = angle(next, prev); + angles = (angles.0.min(a), angles.1.max(a)); + } + let edges = [(0, 1), (1, 2), (2, 3), (3, 0)]; + ( + area, + vec![ + (Metric::ScaledJacobian, scaled), + (Metric::MinAngle, angles.0), + (Metric::MaxAngle, angles.1), + (Metric::EdgeRatio, edge_ratio(&p, &edges)), + ], + ) +} + +/// Wedge scaled Jacobian: the smallest determinant of the three unit edge vectors at a corner, +/// taken in right-handed order (1 for a right prism on a right-angled corner, sin 60° ≈ 0.866 on +/// an equilateral triangle). +pub(crate) fn wedge_scaled_jacobian(p: [P; 6]) -> f64 { + [(0, 1, 2, 3), (1, 2, 0, 4), (2, 0, 1, 5), (3, 5, 4, 0), (4, 3, 5, 1), (5, 4, 3, 2)] + .map(|(o, a, b, c)| corner_det(sub(p[a], p[o]), sub(p[b], p[o]), sub(p[c], p[o]))) + .iter() + .fold(f64::INFINITY, |m, &x| m.min(x)) +} + +/// Scaled Jacobian and Knupp shape of one Jacobian matrix `[a b c]` (hexahedron corners and centre). +pub(crate) fn jacobian_metrics(a: P, b: P, c: P) -> (f64, f64) { + let d = det(a, b, c); + let frob = dot(a, a) + dot(b, b) + dot(c, c); + let shape = if d > 0.0 { ratio(3.0 * d.cbrt().powi(2), frob, 0.0) } else { 0.0 }; + (corner_det(a, b, c), shape) +} diff --git a/crates/mesh/src/quality/tests.rs b/crates/mesh/src/quality/tests.rs new file mode 100644 index 0000000..34ab06a --- /dev/null +++ b/crates/mesh/src/quality/tests.rs @@ -0,0 +1,302 @@ +//! Analytic tests: every expected value is closed-form. + +use modelcraft_model::{ElementKind, Model}; + +use super::*; + +const UNIT_CUBE: [[f64; 3]; 8] = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]]; +const RIGHT_TET: [[f64; 3]; 4] = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]; + +fn close(a: f64, b: f64) -> bool { + (a - b).abs() <= 1e-12 * (1.0 + b.abs()) +} + +fn regular_tet() -> [[f64; 3]; 4] { + let s = 1.0 / 2.0_f64.sqrt(); + // Alternate cube corners: edge √2·(1/√2) = 1. + [[0., 0., 0.], [s, 0., s], [s, s, 0.], [0., s, s]] +} + +fn metric(q: &ElementQuality, m: Metric) -> f64 { + q.get(m).unwrap() +} + +/// Quadratic element nodes from corner nodes with straight edges. +fn midpoints(corners: &[[f64; 3]], edges: &[(usize, usize)]) -> Vec<[f64; 3]> { + let mut out = corners.to_vec(); + for &(a, b) in edges { + out.push([0.5 * (corners[a][0] + corners[b][0]), 0.5 * (corners[a][1] + corners[b][1]), 0.5 * (corners[a][2] + corners[b][2])]); + } + out +} + +#[test] +fn regular_tetrahedron_scores_one() { + let p = regular_tet(); + let q = element_quality(ElementKind::Tet4, &p).unwrap(); + assert!(close(q.size, 1.0 / (6.0 * 2.0_f64.sqrt())), "{}", q.size); + assert_eq!(q.validity(), Some(Validity::Valid)); + let dihedral = (1.0_f64 / 3.0).acos().to_degrees(); + assert!(close(metric(&q, Metric::MinDihedral), dihedral)); + assert!(close(metric(&q, Metric::MaxDihedral), dihedral)); + for m in [Metric::RadiusRatio, Metric::AspectRatio, Metric::VolumeLength, Metric::ScaledJacobian, Metric::JacobianRatio] { + assert!(close(metric(&q, m), 1.0), "{m:?} = {}", metric(&q, m)); + } +} + +#[test] +fn right_corner_tetrahedron() { + let q = element_quality(ElementKind::Tet4, &RIGHT_TET).unwrap(); + assert!(close(q.size, 1.0 / 6.0)); + assert!(close(metric(&q, Metric::MinDihedral), (1.0 / 3.0_f64.sqrt()).acos().to_degrees())); + assert!(close(metric(&q, Metric::MaxDihedral), 90.0)); + // Inradius 1/(3 + √3), circumradius √3/2. + let r = 1.0 / (3.0 + 3.0_f64.sqrt()); + assert!(close(metric(&q, Metric::RadiusRatio), 3.0 * r / (3.0_f64.sqrt() / 2.0))); + // Corner at (1,0,0): edges 1, √2, √2 → 6V/2 · √2. + assert!(close(metric(&q, Metric::ScaledJacobian), 2.0_f64.sqrt() / 2.0)); + // ℓ_rms² = (3·1 + 3·2)/6 = 3/2. + assert!(close(metric(&q, Metric::VolumeLength), 2.0_f64.sqrt() / 1.5_f64.powf(1.5))); +} + +#[test] +fn inverted_and_flat_tetrahedra() { + let mut p = regular_tet(); + p.swap(1, 2); + let q = element_quality(ElementKind::Tet4, &p).unwrap(); + assert!(q.size < 0.0); + assert_eq!(q.validity(), Some(Validity::Invalid)); + assert!(close(metric(&q, Metric::ScaledJacobian), -1.0)); + assert!(close(metric(&q, Metric::RadiusRatio), -1.0)); + assert!(close(metric(&q, Metric::JacobianRatio), -1.0)); + + let flat = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.]]; + let q = element_quality(ElementKind::Tet4, &flat).unwrap(); + assert_eq!(q.size, 0.0); + assert_eq!(q.validity(), Some(Validity::Invalid)); + assert_eq!(metric(&q, Metric::RadiusRatio), 0.0); + assert_eq!(metric(&q, Metric::AspectRatio), f64::INFINITY); + assert!(q.metrics.iter().all(|(_, v)| !v.is_nan())); + + // A sliver with an exactly positive orientation is valid, however thin (exact orient3d). + let sliver = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0.3, 0.3, 1e-30]]; + assert_eq!(element_validity(ElementKind::Tet4, &sliver).unwrap().validity, Validity::Valid); + let below = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0.3, 0.3, -1e-30]]; + assert_eq!(element_validity(ElementKind::Tet4, &below).unwrap().validity, Validity::Invalid); + + // All four nodes at one point. + let q = element_quality(ElementKind::Tet4, &[[1., 2., 3.]; 4]).unwrap(); + assert_eq!(q.validity(), Some(Validity::Invalid)); + assert!(q.metrics.iter().all(|(_, v)| !v.is_nan())); +} + +#[test] +fn unit_cube_scores_one() { + for (kind, nodes) in [(ElementKind::Hex8, UNIT_CUBE.to_vec()), (ElementKind::Hex20, midpoints(&UNIT_CUBE, &jacobian::HEX_EDGES))] { + let q = element_quality(kind, &nodes).unwrap(); + assert!(close(q.size, 1.0), "{kind:?} {}", q.size); + assert_eq!(q.validity(), Some(Validity::Valid)); + for m in [Metric::ScaledJacobian, Metric::JacobianRatio, Metric::Shape, Metric::EdgeRatio] { + assert!(close(metric(&q, m), 1.0), "{kind:?} {m:?} = {}", metric(&q, m)); + } + } +} + +#[test] +fn sheared_hexahedron() { + // x += s·z: a parallelepiped, so det J is constant and the volume stays 1. + let s = 0.75; + let p: Vec<[f64; 3]> = UNIT_CUBE.iter().map(|c| [c[0] + s * c[2], c[1], c[2]]).collect(); + let q = element_quality(ElementKind::Hex8, &p).unwrap(); + assert!(close(q.size, 1.0)); + assert_eq!(q.validity(), Some(Validity::Valid)); + assert!(close(metric(&q, Metric::JacobianRatio), 1.0)); + assert!(close(metric(&q, Metric::ScaledJacobian), 1.0 / (1.0 + s * s).sqrt())); +} + +#[test] +fn mirrored_hexahedron_is_inverted() { + let p: Vec<[f64; 3]> = UNIT_CUBE.iter().map(|c| [-c[0], c[1], c[2]]).collect(); + let q = element_quality(ElementKind::Hex8, &p).unwrap(); + assert!(close(q.size, -1.0)); + assert_eq!(q.validity(), Some(Validity::Invalid)); + assert!(close(metric(&q, Metric::ScaledJacobian), -1.0)); + assert_eq!(metric(&q, Metric::Shape), 0.0); +} + +/// The top face is the bottom face mapped by A = diag(−2, −3) (centred on the axis). The +/// trilinear map is x = ((1−w) I + w A) p, z = w, so det J = (1 − 3w)(1 − 4w): 1 at the bottom +/// corners, 6 at the top corners, 1/2 at the centre, −1/48 at w = 7/24, and the volume is +/// ∫(1 − 7w + 12w²) = 3/2. +fn inside_out_hex() -> Vec<[f64; 3]> { + let bottom = [[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]]; + let mut p: Vec<[f64; 3]> = bottom.iter().map(|b| [b[0], b[1], 0.0]).collect(); + p.extend(bottom.iter().map(|b| [-2.0 * b[0], -3.0 * b[1], 1.0])); + p +} + +#[test] +fn hex_with_positive_corners_inverted_inside_is_invalid() { + let p = inside_out_hex(); + // Every corner Jacobian (three edge vectors in right-handed order) is positive. + let e = |a: usize, b: usize| shape::sub(p[b], p[a]); + for (o, x, y, z) in [(0, 1, 3, 4), (1, 2, 0, 5), (2, 3, 1, 6), (3, 0, 2, 7), (4, 7, 5, 0), (5, 4, 6, 1), (6, 5, 7, 2), (7, 6, 4, 3)] { + assert!(shape::det(e(o, x), e(o, y), e(o, z)) > 0.0, "corner {o}"); + } + let q = element_quality(ElementKind::Hex8, &p).unwrap(); + assert!(metric(&q, Metric::ScaledJacobian) > 0.0, "corners and centre look fine"); + assert!(close(q.size, 1.5)); + assert_eq!(q.validity(), Some(Validity::Invalid)); + assert!(metric(&q, Metric::JacobianRatio) < 0.0); +} + +#[test] +fn quarter_point_hex20_is_degenerate() { + // Move the mid-side node of edge 1-2 (node 9) along the edge to parameter m. On that edge + // x'(0) = 4m − 1, so the element is valid for 1/4 < m < 3/4 and degenerate at m = 1/4. + for (m, valid) in [(0.5, true), (0.3, true), (0.7, true), (0.25, false), (0.2, false), (0.8, false)] { + let mut p = midpoints(&UNIT_CUBE, &jacobian::HEX_EDGES); + p[8] = [m, 0.0, 0.0]; + let check = element_validity(ElementKind::Hex20, &p).unwrap(); + assert_eq!(check.validity == Validity::Valid, valid, "m = {m}: {check:?}"); + assert!(check.validity != Validity::Undetermined); + } +} + +#[test] +fn tet10_with_a_moved_mid_side_node() { + // Moving node 5 (edge 1-2) to parameter m gives det J = det_lin (1 + 4(m − ½)(λ1 − λ2)), + // extreme at the two edge ends: valid iff |m − ½| < ¼, ratio (1 − 4d)/(1 + 4d). + for m in [0.5, 0.3, 0.6, 0.25, 0.1, 0.9] { + let mut p = midpoints(&RIGHT_TET, &jacobian::TET_EDGES); + p[4] = [m, 0.0, 0.0]; + let q = element_quality(ElementKind::Tet10, &p).unwrap(); + let d: f64 = (m - 0.5_f64).abs(); + let valid = d < 0.25; + assert_eq!(q.validity() == Some(Validity::Valid), valid, "m = {m}"); + if valid { + assert!(close(metric(&q, Metric::JacobianRatio), (1.0 - 4.0 * d) / (1.0 + 4.0 * d)), "m = {m}"); + } + // The node moves along a straight edge, so the volume is unchanged. + assert!(close(q.size, 1.0 / 6.0), "m = {m}: {}", q.size); + } +} + +#[test] +fn wedges() { + let right = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]]; + let q = element_quality(ElementKind::Wedge6, &right).unwrap(); + assert!(close(q.size, 0.5)); + assert_eq!(q.validity(), Some(Validity::Valid)); + assert!(close(metric(&q, Metric::ScaledJacobian), 0.5_f64.sqrt())); + assert!(close(metric(&q, Metric::JacobianRatio), 1.0)); + assert!(close(metric(&q, Metric::EdgeRatio), 2.0_f64.sqrt())); + + let h = 3.0_f64.sqrt() / 2.0; + let equilateral = [[0., 0., 0.], [1., 0., 0.], [0.5, h, 0.], [0., 0., 1.], [1., 0., 1.], [0.5, h, 1.]]; + let q = element_quality(ElementKind::Wedge6, &equilateral).unwrap(); + assert!(close(metric(&q, Metric::ScaledJacobian), h)); + + let upside_down: Vec<[f64; 3]> = right.iter().map(|p| [p[0], p[1], -p[2]]).collect(); + let q = element_quality(ElementKind::Wedge6, &upside_down).unwrap(); + assert!(close(q.size, -0.5)); + assert_eq!(q.validity(), Some(Validity::Invalid)); +} + +#[test] +fn shells() { + let h = 3.0_f64.sqrt() / 2.0; + let q = element_quality(ElementKind::Tri3, &[[0., 0., 0.], [1., 0., 0.], [0.5, h, 0.]]).unwrap(); + assert!(close(q.size, h / 2.0)); + assert!(close(metric(&q, Metric::MinAngle), 60.0) && close(metric(&q, Metric::MaxAngle), 60.0)); + assert!(close(metric(&q, Metric::AspectRatio), 1.0)); + assert_eq!(q.validity(), None); + + let square = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]]; + let q = element_quality(ElementKind::Quad4, &square).unwrap(); + assert!(close(q.size, 1.0)); + assert!(close(metric(&q, Metric::ScaledJacobian), 1.0)); + assert!(close(metric(&q, Metric::MinAngle), 90.0)); + // A concave quadrilateral: corner 3 pushed past the diagonal. + let dart = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0.7, 0.3, 0.]]; + assert!(metric(&element_quality(ElementKind::Quad4, &dart).unwrap(), Metric::ScaledJacobian) < 0.0); + + assert!(element_quality(ElementKind::Rod2, &[[0.; 3], [1., 0., 0.]]).is_err()); +} + +#[test] +fn hostile_coordinates_are_errors() { + let mut p = UNIT_CUBE; + p[3][1] = f64::NAN; + assert!(element_quality(ElementKind::Hex8, &p).is_err()); + p[3][1] = 1e300; + assert!(element_validity(ElementKind::Hex8, &p).is_err()); + assert!(element_quality(ElementKind::Hex8, &UNIT_CUBE[..7]).is_err()); + assert!(element_validity(ElementKind::Tri3, &[[0.; 3]; 3]).is_err()); + // Far from the origin and at the coordinate limit: still exact. + let far: Vec<[f64; 3]> = UNIT_CUBE.iter().map(|c| [c[0] * 1e9, c[1] * 1e9, c[2] * 1e9 - 1e9]).collect(); + let q = element_quality(ElementKind::Hex8, &far).unwrap(); + assert_eq!(q.validity(), Some(Validity::Valid)); + assert!(close(q.size, 1e27)); +} + +fn model_of(elements: &[(ElementKind, Vec<[f64; 3]>)]) -> Model { + let mut m = Model::default(); + let mut next = 1; + for (e, (kind, nodes)) in elements.iter().enumerate() { + let mut labels = Vec::new(); + for p in nodes { + m.set_node(next, *p).unwrap(); + labels.push(next); + next += 1; + } + m.set_element(e as u32 + 1, *kind, labels).unwrap(); + } + m +} + +#[test] +fn report_and_gate() { + let mut flipped = UNIT_CUBE.to_vec(); + flipped.iter_mut().for_each(|p| p[0] = -p[0] - 2.0); + let model = model_of(&[ + (ElementKind::Hex8, UNIT_CUBE.to_vec()), + (ElementKind::Hex8, inside_out_hex()), + (ElementKind::Tet4, regular_tet().to_vec()), + (ElementKind::Hex8, flipped), + (ElementKind::Rod2, vec![[0.; 3], [1., 0., 0.]]), + ]); + let r = quality_report(&model, None, 4).unwrap(); + assert_eq!(r.measured, 4); + assert_eq!(r.skipped, vec![(ElementKind::Rod2, 1)]); + assert_eq!((r.valid, r.invalid.clone()), (2, vec![2, 4])); + let sj = r.metrics.iter().find(|s| s.metric == Metric::ScaledJacobian).unwrap(); + assert_eq!((sj.count, sj.worst, sj.min, sj.histogram.iter().sum::()), (4, 4, -1.0, 4)); + assert_eq!(sj.histogram[0], 1); + assert_eq!(r.metrics.iter().find(|s| s.metric == Metric::RadiusRatio).unwrap().count, 1); + assert!(quality_report(&model, None, 0).is_err()); + assert!(quality_report(&model, Some(&[99]), 4).is_err()); + assert_eq!(quality_report(&model, Some(&[1, 1, 3]), 4).unwrap().measured, 2); + + let gate = validity_gate(&model, &GateOptions::default()).unwrap(); + assert_eq!((gate.checked, gate.invalid.clone()), (4, vec![2, 4])); + assert!(gate.into_result().is_err()); + // The same gate straight from a mesher's arrays. + let raw = validity_gate_elements([(7, ElementKind::Hex8, inside_out_hex()), (8, ElementKind::Hex8, UNIT_CUBE.to_vec())], &GateOptions::default()) + .unwrap(); + assert_eq!((raw.checked, raw.invalid), (2, vec![7])); + + // A unit cube split into six tetrahedra around the 1-7 diagonal passes against its boundary. + let c = UNIT_CUBE; + let tets = [[0, 1, 2, 6], [0, 2, 3, 6], [0, 3, 7, 6], [0, 7, 4, 6], [0, 4, 5, 6], [0, 5, 1, 6]]; + let model = model_of(&tets.map(|t| (ElementKind::Tet4, t.iter().map(|&i| c[i]).collect()))); + let faces = [[0, 3, 2], [0, 2, 1], [4, 5, 6], [4, 6, 7], [0, 1, 5], [0, 5, 4], [1, 2, 6], [1, 6, 5], [2, 3, 7], [2, 7, 6], [3, 0, 4], [3, 4, 7]]; + let boundary = enclosed_volume(&c, &faces).unwrap(); + assert!(close(boundary, 1.0)); + let ok = validity_gate(&model, &GateOptions { expected_volume: Some(boundary), ..Default::default() }).unwrap(); + assert!(ok.passed(), "{ok:?}"); + let short = validity_gate(&model, &GateOptions { expected_volume: Some(1.1), ..Default::default() }).unwrap(); + assert!(!short.volume_matches() && short.into_result().is_err()); + assert!(validity_gate(&model, &GateOptions { expected_volume: Some(f64::NAN), ..Default::default() }).is_err()); + assert!(enclosed_volume(&c, &[[0, 1, 99]]).is_err()); +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 8dc6612..e4454b2 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -47,9 +47,17 @@ pub struct Node { pub xyz: [f64; 3], } -/// Element topology. Node order follows each family's usual corner-first convention (corners, -/// then mid-side nodes), which all three target formats share for these shapes. -#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)] +/// Element topology. Node order is corners first, then mid-side nodes, with positive volume for +/// a right-handed ordering: +/// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; +/// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to; +/// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the +/// edges of 1-2-3-4, 13–16 on those of 5-6-7-8, and 17–20 on edges 1-5, 2-6, 3-7, 4-8. +/// +/// This is the Abaqus (and VTK) order. Formats that order nodes differently, such as Nastran's +/// CHEXA with 20 nodes (vertical edges before the top face's edges), are reordered by their deck +/// readers and writers; the model never stores another order. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] #[serde(rename_all = "snake_case")] pub enum ElementKind { /// Two-node axial member. @@ -148,6 +156,15 @@ pub struct Model { elements: BTreeMap, } +/// Check that every coordinate is a finite number within ±[`MAX_COORD`] (the rule for nodes, +/// shared by everything that takes coordinates from outside). +pub fn check_coordinates(xyz: &[f64]) -> Result<()> { + match xyz.iter().find(|c| !(c.abs() <= MAX_COORD)) { + Some(c) => Err(Error::Invalid(format!("coordinate {c} is not a finite number within ±{MAX_COORD}"))), + None => Ok(()), + } +} + fn check_label(what: &str, label: u32) -> Result<()> { if label == 0 || label > MAX_LABEL { return Err(Error::Invalid(format!("{what} label {label} is outside 1..={MAX_LABEL}"))); @@ -175,9 +192,10 @@ impl Model { /// Add a node, or move an existing one (decks redefine nodes the same way). pub fn set_node(&mut self, label: u32, xyz: [f64; 3]) -> Result<()> { check_label("node", label)?; - if let Some(c) = xyz.iter().find(|c| !(c.abs() <= MAX_COORD)) { - return Err(Error::Invalid(format!("node {label}: coordinate {c} is not a finite number within ±{MAX_COORD}"))); - } + check_coordinates(&xyz).map_err(|e| match e { + Error::Invalid(m) => Error::Invalid(format!("node {label}: {m}")), + e => e, + })?; if !self.nodes.contains_key(&label) && self.nodes.len() >= MAX_NODES { return Err(Error::Limit(format!("a model holds at most {MAX_NODES} nodes"))); } diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index ffbb5f4..89821a7 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -746,6 +746,19 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant - **HexaLab** (MIT, https://github.com/mlivesu/HexaLab, https://www.hexalab.net) is READ OK. It is a hex quality viewer whose metrics match the literature's. +**B6.6 Element validity and quality metrics (implemented in `modelcraft-mesh::quality`).** **Verdict: BUILD ON.** +- **Johnen, Remacle, Geuzaine 2013.** *Geometrical validity of curvilinear finite elements.* J. Comput. Phys. + 233:359–372, 2013. The Jacobian + determinant in Bernstein form is bounded by its coefficients, equals them at the vertices and converges under + subdivision: a proof of `det J > 0` on the whole element. We use it for every solid kind, including the trilinear + hex whose corner Jacobians do not certify it (B6.3). (The paper's authors also wrote Gmsh code, which stays + forbidden; the paper is enough.) +- **Shewchuk 2002.** *What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality + measures.* Unpublished preprint, https://people.eecs.berkeley.edu/~jrs/papers/elemj.pdf. Tet dihedral angles, + radius ratio and aspect ratio; with Klingner–Shewchuk's volume–length measure (A3.6). +- **Knupp 2001.** *Algebraic mesh quality metrics.* SIAM J. Sci. Comput. 23(1):193–218. Corner-Jacobian scaled Jacobian and + shape metrics for hexes, wedges and quads. + --- ## Part C: Datasets and benchmarks