conspire/geometry/mesh/quality/metrics/
mod.rs

1mod hexahedron;
2mod quadrilateral;
3mod tetrahedron;
4mod triangle;
5
6use crate::{
7    geometry::{
8        Coordinate, Coordinates,
9        mesh::{Connectivity, Mesh},
10    },
11    math::{Scalar, Tensor, TensorRank2},
12};
13use std::array::from_fn;
14
15pub trait Verdict {
16    fn maximum_edge_ratios(&self) -> Vec<Vec<Scalar>>;
17    fn maximum_skews(&self) -> Vec<Vec<Scalar>>;
18    fn minimum_jacobians(&self) -> Vec<Vec<Scalar>>;
19    fn minimum_scaled_jacobians(&self) -> Vec<Vec<Scalar>>;
20    fn volumes(&self) -> Vec<Vec<Scalar>>;
21}
22
23impl<const D: usize> Verdict for Mesh<D> {
24    fn maximum_edge_ratios(&self) -> Vec<Vec<Scalar>> {
25        let coordinates = self.coordinates();
26        self.iter()
27            .map(|block| match block {
28                Connectivity::Triangular(elements) => elements
29                    .iter()
30                    .map(|element| triangle::maximum_edge_ratio(element, coordinates))
31                    .collect(),
32                Connectivity::Quadrilateral(elements) => elements
33                    .iter()
34                    .map(|element| quadrilateral::maximum_edge_ratio(element, coordinates))
35                    .collect(),
36                Connectivity::Tetrahedral(elements) => elements
37                    .iter()
38                    .map(|element| tetrahedron::maximum_edge_ratio(element, coordinates))
39                    .collect(),
40                Connectivity::Hexahedral(elements) => elements
41                    .iter()
42                    .map(|element| hexahedron::maximum_edge_ratio(element, coordinates))
43                    .collect(),
44                _ => todo!(),
45            })
46            .collect()
47    }
48    fn minimum_jacobians(&self) -> Vec<Vec<Scalar>> {
49        let coordinates = self.coordinates();
50        self.iter()
51            .map(|block| match block {
52                Connectivity::Triangular(elements) => elements
53                    .iter()
54                    .map(|element| triangle::minimum_jacobian(element, coordinates))
55                    .collect(),
56                Connectivity::Quadrilateral(elements) => elements
57                    .iter()
58                    .map(|element| quadrilateral::minimum_jacobian(element, coordinates))
59                    .collect(),
60                Connectivity::Tetrahedral(elements) => elements
61                    .iter()
62                    .map(|element| tetrahedron::minimum_jacobian(element, coordinates))
63                    .collect(),
64                Connectivity::Hexahedral(elements) => elements
65                    .iter()
66                    .map(|element| hexahedron::minimum_jacobian(element, coordinates))
67                    .collect(),
68                _ => todo!(),
69            })
70            .collect()
71    }
72    fn minimum_scaled_jacobians(&self) -> Vec<Vec<Scalar>> {
73        let coordinates = self.coordinates();
74        self.iter()
75            .map(|block| match block {
76                Connectivity::Triangular(elements) => elements
77                    .iter()
78                    .map(|element| triangle::minimum_scaled_jacobian(element, coordinates))
79                    .collect(),
80                Connectivity::Quadrilateral(elements) => elements
81                    .iter()
82                    .map(|element| quadrilateral::minimum_scaled_jacobian(element, coordinates))
83                    .collect(),
84                Connectivity::Tetrahedral(elements) => elements
85                    .iter()
86                    .map(|element| tetrahedron::minimum_scaled_jacobian(element, coordinates))
87                    .collect(),
88                Connectivity::Hexahedral(elements) => elements
89                    .iter()
90                    .map(|element| hexahedron::minimum_scaled_jacobian(element, coordinates))
91                    .collect(),
92                _ => todo!(),
93            })
94            .collect()
95    }
96    fn maximum_skews(&self) -> Vec<Vec<Scalar>> {
97        let coordinates = self.coordinates();
98        self.iter()
99            .map(|block| match block {
100                Connectivity::Triangular(elements) => elements
101                    .iter()
102                    .map(|element| triangle::maximum_skew(element, coordinates))
103                    .collect(),
104                Connectivity::Quadrilateral(elements) => elements
105                    .iter()
106                    .map(|element| quadrilateral::maximum_skew(element, coordinates))
107                    .collect(),
108                Connectivity::Tetrahedral(elements) => elements
109                    .iter()
110                    .map(|element| tetrahedron::maximum_skew(element, coordinates))
111                    .collect(),
112                Connectivity::Hexahedral(elements) => elements
113                    .iter()
114                    .map(|element| hexahedron::maximum_skew(element, coordinates))
115                    .collect(),
116                _ => todo!(),
117            })
118            .collect()
119    }
120    fn volumes(&self) -> Vec<Vec<Scalar>> {
121        let coordinates = self.coordinates();
122        self.iter()
123            .map(|block| match block {
124                Connectivity::Triangular(elements) => elements
125                    .iter()
126                    .map(|element| triangle::volume(element, coordinates))
127                    .collect(),
128                Connectivity::Quadrilateral(elements) => elements
129                    .iter()
130                    .map(|element| quadrilateral::volume(element, coordinates))
131                    .collect(),
132                Connectivity::Tetrahedral(elements) => elements
133                    .iter()
134                    .map(|element| tetrahedron::volume(element, coordinates))
135                    .collect(),
136                Connectivity::Hexahedral(elements) => elements
137                    .iter()
138                    .map(|element| hexahedron::volume(element, coordinates))
139                    .collect(),
140                _ => todo!(),
141            })
142            .collect()
143    }
144}
145
146const EQUIANGLE: Scalar = std::f64::consts::FRAC_PI_3;
147
148fn cross<const D: usize>(a: &Coordinate<D>, b: &Coordinate<D>) -> [Scalar; 3] {
149    let az = if D > 2 { a[2] } else { 0.0 };
150    let bz = if D > 2 { b[2] } else { 0.0 };
151    [
152        a[1] * bz - az * b[1],
153        az * b[0] - a[0] * bz,
154        a[0] * b[1] - a[1] * b[0],
155    ]
156}
157
158fn triple_product<const D: usize>(
159    a: &Coordinate<D>,
160    b: &Coordinate<D>,
161    c: &Coordinate<D>,
162) -> Scalar {
163    let bc = cross(b, c);
164    a[0] * bc[0] + a[1] * bc[1] + a[2] * bc[2]
165}
166
167fn cross_norm<const D: usize>(a: &Coordinate<D>, b: &Coordinate<D>) -> Scalar {
168    let n = cross(a, b);
169    (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt()
170}
171
172fn triangle_area<const D: usize>(triangle: &[usize; 3], coordinates: &Coordinates<D>) -> Scalar {
173    let a = &coordinates[triangle[1]] - &coordinates[triangle[0]];
174    let b = &coordinates[triangle[2]] - &coordinates[triangle[0]];
175    0.5 * cross_norm(&a, &b)
176}
177
178fn tet_volume<const D: usize>(tetrahedron: &[usize; 4], coordinates: &Coordinates<D>) -> Scalar {
179    let a = &coordinates[tetrahedron[1]] - &coordinates[tetrahedron[0]];
180    let b = &coordinates[tetrahedron[2]] - &coordinates[tetrahedron[0]];
181    let c = &coordinates[tetrahedron[3]] - &coordinates[tetrahedron[0]];
182    triple_product(&a, &b, &c) / 6.0
183}
184
185fn triangle_skew<const D: usize>(
186    a: &Coordinate<D>,
187    b: &Coordinate<D>,
188    c: &Coordinate<D>,
189) -> Scalar {
190    let l0 = (c - b).normalized();
191    let l1 = (a - c).normalized();
192    let l2 = (b - a).normalized();
193    let minimum_angle = [
194        (-(&l0 * &l1)).acos(),
195        (-(&l1 * &l2)).acos(),
196        (-(&l2 * &l0)).acos(),
197    ]
198    .into_iter()
199    .fold(Scalar::INFINITY, Scalar::min);
200    (EQUIANGLE - minimum_angle) / EQUIANGLE
201}
202
203fn maximum_edge_ratio<const D: usize, const E: usize>(
204    edges: &[[usize; 2]; E],
205    element: &[usize],
206    coordinates: &Coordinates<D>,
207) -> Scalar {
208    let mut shortest = Scalar::INFINITY;
209    let mut longest: Scalar = 0.0;
210    for [a, b] in edges {
211        let length = (&coordinates[element[*b]] - &coordinates[element[*a]]).norm();
212        shortest = shortest.min(length);
213        longest = longest.max(length);
214    }
215    if shortest > 0.0 {
216        longest / shortest
217    } else {
218        Scalar::INFINITY
219    }
220}
221
222fn min_jacobian<const D: usize, const K: usize, const C: usize>(
223    table: &[[usize; K]; C],
224    element: &[usize],
225    coordinates: &Coordinates<D>,
226) -> Scalar {
227    corners(table, element, coordinates)
228        .into_iter()
229        .map(|(measure, _)| measure)
230        .fold(Scalar::INFINITY, Scalar::min)
231}
232
233fn min_scaled_jacobian<const D: usize, const K: usize, const C: usize>(
234    table: &[[usize; K]; C],
235    element: &[usize],
236    coordinates: &Coordinates<D>,
237    scale: Scalar,
238) -> Scalar {
239    corners(table, element, coordinates)
240        .into_iter()
241        .map(|(measure, normalizer)| {
242            if normalizer > 0.0 {
243                (scale * measure / normalizer).clamp(-1.0, 1.0)
244            } else {
245                0.0
246            }
247        })
248        .fold(Scalar::INFINITY, Scalar::min)
249}
250
251fn corners<const D: usize, const K: usize, const C: usize>(
252    table: &[[usize; K]; C],
253    element: &[usize],
254    coordinates: &Coordinates<D>,
255) -> [(Scalar, Scalar); C] {
256    from_fn(|corner| {
257        let origin = &coordinates[element[corner]];
258        let edges: [Coordinate<D>; K] =
259            from_fn(|edge| &coordinates[element[table[corner][edge]]] - origin);
260        let normalizer: Scalar = edges.iter().map(|edge| edge.norm()).product();
261        (corner_measure(&edges), normalizer)
262    })
263}
264
265fn corner_measure<const D: usize, const K: usize>(edges: &[Coordinate<D>; K]) -> Scalar {
266    if K == D {
267        let matrix: [[Scalar; K]; K] = from_fn(|row| from_fn(|column| edges[row][column]));
268        TensorRank2::<K, 0, 0>::from(matrix).determinant()
269    } else {
270        let gram: [[Scalar; K]; K] = from_fn(|i| from_fn(|j| &edges[i] * &edges[j]));
271        TensorRank2::<K, 0, 0>::from(gram)
272            .determinant()
273            .max(0.0)
274            .sqrt()
275    }
276}
277
278#[derive(Clone, Copy)]
279pub(super) enum Kind {
280    Triangle,
281    Quadrilateral,
282    Tetrahedron,
283    Hexahedron,
284}
285
286impl Kind {
287    pub(super) fn of(connectivity: &Connectivity) -> Option<Self> {
288        match connectivity {
289            Connectivity::Triangular(_) => Some(Self::Triangle),
290            Connectivity::Quadrilateral(_) => Some(Self::Quadrilateral),
291            Connectivity::Tetrahedral(_) => Some(Self::Tetrahedron),
292            Connectivity::Hexahedral(_) => Some(Self::Hexahedron),
293            _ => None,
294        }
295    }
296}
297
298pub(super) fn minimum_jacobian<const D: usize>(
299    kind: Kind,
300    element: &[usize],
301    coordinates: &Coordinates<D>,
302) -> Scalar {
303    match kind {
304        Kind::Triangle => triangle::minimum_jacobian(element, coordinates),
305        Kind::Quadrilateral => quadrilateral::minimum_jacobian(element, coordinates),
306        Kind::Tetrahedron => tetrahedron::minimum_jacobian(element, coordinates),
307        Kind::Hexahedron => hexahedron::minimum_jacobian(element, coordinates),
308    }
309}
310
311pub(super) fn minimum_scaled_jacobian<const D: usize>(
312    kind: Kind,
313    element: &[usize],
314    coordinates: &Coordinates<D>,
315) -> Scalar {
316    match kind {
317        Kind::Triangle => triangle::minimum_scaled_jacobian(element, coordinates),
318        Kind::Quadrilateral => quadrilateral::minimum_scaled_jacobian(element, coordinates),
319        Kind::Tetrahedron => tetrahedron::minimum_scaled_jacobian(element, coordinates),
320        Kind::Hexahedron => hexahedron::minimum_scaled_jacobian(element, coordinates),
321    }
322}
conspire/geometry/mesh/quality/metrics/mod.rs

conspire/geometry/mesh/quality/metrics/
mod.rs
