forked from 1b-agent/fieldcraft-upstream
Element quality metrics and the certified validity gate (mesh.quality) (#12)
* mesh: element quality metrics and a certified Jacobian validity check Bernstein-form Jacobian determinants for tet4/tet10/wedge6/hex8/hex20, certified by bounds and subdivision (Johnen-Remacle-Geuzaine 2013); Shewchuk 2002 tet metrics, Knupp 2001 hex/wedge/quad metrics; quality report and the validity gate meshers run on their output. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * engine: mesh.quality command; model documents its node order mesh.quality reports certified validity, volume and per-metric statistics with fixed-range histograms. The model's ElementKind doc states the Abaqus/VTK node order and that deck code reorders formats like Nastran CHEXA20. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * quality: simplify after review; exact orient3d for tet4; gate on raw mesher arrays - tet4 validity from the exact orient3d of the input coordinates - validity_gate_elements takes (label, kind, coordinates) straight from a mesher - one Jacobian per solid; certify caches piece bounds; product weights precomputed - model::check_coordinates shared by set_node and quality; ElementKind is Ord - engine: shared labels() parameter parser Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
15 files changed
+1981
-17
No files matched your search
@@ -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}")]
|
||||
|
||||
@@ -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<Factor>,
|
||||
coef: Vec<f64>,
|
||||
}
|
||||
|
||||
fn factorial(n: usize) -> f64 {
|
||||
(2..=n).map(|k| k as f64).product()
|
||||
}
|
||||
|
||||
impl Bernstein {
|
||||
pub fn zero(factors: &[Factor]) -> Result<Self> {
|
||||
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<Multi> {
|
||||
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<f64> {
|
||||
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<Self> {
|
||||
if !self.same_domain(o) {
|
||||
return Err(Error::Invalid("Bernstein product of different domains".into()));
|
||||
}
|
||||
let factors: Vec<Factor> = 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<Self> {
|
||||
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<Self> {
|
||||
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<Multi> = (0..=s).map(|r| shift(&start, r)).collect();
|
||||
let mut tri: Vec<f64> = 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<Vec<usize>> {
|
||||
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<f64> {
|
||||
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::<f64>() * 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<Vec<f64>> {
|
||||
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<Factor>> {
|
||||
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<Vec<Vec<f64>>> {
|
||||
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>]) -> 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::<Vec<_>>(), 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<Vec<f64>> = 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());
|
||||
}
|
||||
}
|
||||
@@ -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<Factor>, [(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<Jacobian> {
|
||||
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<Bernstein> { 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<Vec<[f64; 3]>>,
|
||||
depth: u32,
|
||||
}
|
||||
|
||||
impl Piece {
|
||||
fn new(det: Bernstein, verts: Vec<Vec<[f64; 3]>>, 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<JacobianCheck> {
|
||||
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 })
|
||||
}
|
||||
@@ -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<JacobianCheck>,
|
||||
pub metrics: Vec<(Metric, f64)>,
|
||||
}
|
||||
|
||||
impl ElementQuality {
|
||||
pub fn get(&self, metric: Metric) -> Option<f64> {
|
||||
self.metrics.iter().find(|(m, _)| *m == metric).map(|(_, v)| *v)
|
||||
}
|
||||
|
||||
pub fn validity(&self) -> Option<Validity> {
|
||||
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<Vec<[f64; 3]>> {
|
||||
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<const N: usize>(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<Columns>, 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<Sign> {
|
||||
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<JacobianCheck> {
|
||||
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<ElementQuality> {
|
||||
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;
|
||||
@@ -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<f64>,
|
||||
/// 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<usize>,
|
||||
}
|
||||
|
||||
/// 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<u32>,
|
||||
pub undetermined: Vec<u32>,
|
||||
/// 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<MetricStats>,
|
||||
}
|
||||
|
||||
fn coordinates(model: &Model, label: u32, element: &Element) -> Result<Vec<[f64; 3]>> {
|
||||
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<Vec<(u32, &'a Element)>> {
|
||||
match labels {
|
||||
None => Ok(model.elements().iter().map(|(l, e)| (*l, e)).collect()),
|
||||
Some(labels) => labels
|
||||
.iter()
|
||||
.collect::<BTreeSet<_>>()
|
||||
.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<QualityReport> {
|
||||
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<ElementKind, usize> = BTreeMap::new();
|
||||
let mut stats: BTreeMap<Metric, (MetricStats, f64, usize)> = 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<f64>,
|
||||
/// 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<u32>,
|
||||
pub undetermined: Vec<u32>,
|
||||
/// 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<Self> {
|
||||
if self.passed() {
|
||||
return Ok(self);
|
||||
}
|
||||
let mut why = Vec::new();
|
||||
let list = |v: &[u32]| v.iter().take(10).map(u32::to_string).collect::<Vec<_>>().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<Self> {
|
||||
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<I, E>(elements: I, options: &GateOptions) -> Result<GateReport>
|
||||
where
|
||||
I: IntoIterator<Item = (u32, ElementKind, E)>,
|
||||
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<GateReport> {
|
||||
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<f64> {
|
||||
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)
|
||||
}
|
||||
@@ -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::<f64>() / 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)
|
||||
}
|
||||
@@ -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::<usize>()), (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());
|
||||
}
|
||||
Reference in new issue
Block a user