carapace_core/model/continuum/
validate.rs

1//! Geometry validation with a scale-aware tolerance: every threshold is
2//! relative to the element's characteristic length squared (its longest
3//! node-to-node distance squared), never an absolute number, so a unit
4//! element and a 1e6-scaled copy are judged identically.
5
6/// A distance below `COINCIDENT * L` is "the same node".
7const COINCIDENT: f64 = 1e-8;
8/// A determinant (or twice the triangle area) below `DEGENERATE * L^2` is zero.
9const DEGENERATE: f64 = 1e-10;
10
11fn length_squared<const N: usize>(p: &[[f64; 2]; N]) -> f64 {
12    let mut longest = 0.0_f64;
13    for i in 0..N {
14        for j in i + 1..N {
15            let (dx, dy) = (p[j][0] - p[i][0], p[j][1] - p[i][1]);
16            longest = longest.max(dx * dx + dy * dy);
17        }
18    }
19    longest
20}
21
22fn check_distinct<const N: usize>(p: &[[f64; 2]; N], l2: f64) -> Result<(), &'static str> {
23    for i in 0..N {
24        for j in i + 1..N {
25            let (dx, dy) = (p[j][0] - p[i][0], p[j][1] - p[i][1]);
26            if dx * dx + dy * dy <= COINCIDENT * COINCIDENT * l2 {
27                return Err("repeated or coincident nodes");
28            }
29        }
30    }
31    Ok(())
32}
33
34/// A triangle needs distinct nodes, area above tolerance, and counter-clockwise order.
35pub fn validate_triangle(p: &[[f64; 2]; 3]) -> Result<(), &'static str> {
36    let l2 = length_squared(p);
37    if !l2.is_finite() || l2 == 0.0 {
38        return Err("repeated or coincident nodes");
39    }
40    check_distinct(p, l2)?;
41    let twice_area =
42        (p[1][0] - p[0][0]) * (p[2][1] - p[0][1]) - (p[2][0] - p[0][0]) * (p[1][1] - p[0][1]);
43    if twice_area.abs() <= DEGENERATE * l2 {
44        return Err("triangle area is zero (collinear nodes)");
45    }
46    if twice_area < 0.0 {
47        return Err("triangle nodes are ordered clockwise");
48    }
49    Ok(())
50}
51
52/// `det J` at reference corner `k` of a bilinear quad, from the two edge
53/// vectors leaving the corner (`J = [e1; e2] / 2` there).
54fn corner_det(p: &[[f64; 2]; 4], k: usize) -> f64 {
55    let (next, prev) = (p[(k + 1) % 4], p[(k + 3) % 4]);
56    let (e1, e2) = (
57        [next[0] - p[k][0], next[1] - p[k][1]],
58        [prev[0] - p[k][0], prev[1] - p[k][1]],
59    );
60    // Counter-clockwise order: the edge to the next node cross the edge to the previous is positive.
61    0.25 * (e1[0] * e2[1] - e1[1] * e2[0])
62}
63
64/// A quad needs distinct nodes and a positive `det J` at all four reference corners.
65///
66/// `det J` of a bilinear map is linear in `(xi, eta)`, so positivity at the
67/// corners is necessary and sufficient for positivity everywhere; positivity
68/// at the Gauss points alone is not (a quad with one reflex corner can have
69/// all Gauss-point determinants positive).
70pub fn validate_quad(p: &[[f64; 2]; 4]) -> Result<(), &'static str> {
71    let l2 = length_squared(p);
72    if !l2.is_finite() || l2 == 0.0 {
73        return Err("repeated or coincident nodes");
74    }
75    check_distinct(p, l2)?;
76    let tol = DEGENERATE * l2;
77    let dets: [f64; 4] = std::array::from_fn(|k| corner_det(p, k));
78    if dets.iter().all(|&d| d < -tol) {
79        return Err("quad nodes are ordered clockwise");
80    }
81    const MESSAGES: [&str; 4] = [
82        "quad corner 0 has a non-positive Jacobian determinant (reflex or degenerate corner)",
83        "quad corner 1 has a non-positive Jacobian determinant (reflex or degenerate corner)",
84        "quad corner 2 has a non-positive Jacobian determinant (reflex or degenerate corner)",
85        "quad corner 3 has a non-positive Jacobian determinant (reflex or degenerate corner)",
86    ];
87    for (k, &d) in dets.iter().enumerate() {
88        if d <= tol {
89            return Err(MESSAGES[k]);
90        }
91    }
92    Ok(())
93}
94
95#[cfg(test)]
96mod tests {
97    use super::*;
98    use crate::model::continuum::{quad4_dshape, Jacobian, QUAD4_GAUSS_2X2};
99
100    fn scaled<const N: usize>(p: [[f64; 2]; N], s: f64) -> [[f64; 2]; N] {
101        p.map(|[x, y]| [x * s, y * s])
102    }
103
104    #[test]
105    fn the_reflex_corner_quad_is_rejected_though_its_gauss_point_determinants_are_positive() {
106        let reflex = [[0.0, 0.0], [1.0, 0.0], [0.4, 0.4], [0.0, 1.0]];
107        let mut gauss: Vec<f64> = QUAD4_GAUSS_2X2
108            .iter()
109            .map(|&(xi, eta, _)| Jacobian::new(&reflex, &quad4_dshape(xi, eta)).det)
110            .collect();
111        gauss.sort_by(f64::total_cmp);
112        assert!(gauss.iter().all(|&d| d > 0.0), "{gauss:?}");
113        for (got, expected) in gauss.iter().zip([0.013, 0.1, 0.1, 0.187]) {
114            assert!((got - expected).abs() < 1e-3, "{gauss:?}");
115        }
116        let corner = Jacobian::new(&reflex, &quad4_dshape(1.0, 1.0)).det;
117        assert!((corner + 0.05).abs() < 1e-12);
118        let error = validate_quad(&reflex).unwrap_err();
119        assert!(error.contains("corner 2"), "{error}");
120    }
121
122    #[test]
123    fn convex_trapezoid_and_parallelogram_are_accepted() {
124        assert!(validate_quad(&[[0.0, 0.0], [2.0, 0.0], [1.5, 1.0], [0.5, 1.0]]).is_ok());
125        assert!(validate_quad(&[[0.0, 0.0], [2.0, 0.0], [3.0, 1.5], [1.0, 1.5]]).is_ok());
126        assert!(validate_triangle(&[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]).is_ok());
127    }
128
129    #[test]
130    fn repeated_clockwise_and_degenerate_inputs_are_rejected() {
131        let quad = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]];
132        let clockwise = [quad[0], quad[3], quad[2], quad[1]];
133        assert!(validate_quad(&clockwise).unwrap_err().contains("clockwise"));
134        let repeated = [quad[0], quad[1], quad[1], quad[3]];
135        assert!(validate_quad(&repeated).unwrap_err().contains("coincident"));
136        assert!(validate_quad(&[[0.0, 0.0]; 4]).is_err());
137        // A bow-tie (self-intersecting) quad.
138        assert!(validate_quad(&[[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]).is_err());
139
140        let tri = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]];
141        assert!(validate_triangle(&[tri[0], tri[2], tri[1]])
142            .unwrap_err()
143            .contains("clockwise"));
144        assert!(validate_triangle(&[tri[0], tri[1], tri[1]])
145            .unwrap_err()
146            .contains("coincident"));
147        assert!(validate_triangle(&[[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]])
148            .unwrap_err()
149            .contains("collinear"));
150    }
151
152    #[test]
153    fn the_tolerance_is_scale_invariant() {
154        let near_degenerate = [[0.0, 0.0], [1.0, 0.0], [0.5, 1e-11]];
155        let sliver_quad = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]];
156        for s in [1.0, 1e-6, 1e6] {
157            assert!(
158                validate_triangle(&scaled(near_degenerate, s)).is_err(),
159                "scale {s}"
160            );
161            assert!(validate_quad(&scaled(sliver_quad, s)).is_ok(), "scale {s}");
162            assert!(
163                validate_triangle(&scaled([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], s)).is_ok(),
164                "scale {s}"
165            );
166        }
167    }
168}
carapace_core/model/continuum/validate.rs

carapace_core/model/continuum/
validate.rs
