carapace_core/model/continuum/
quadrature.rs

1//! Gauss rules as `(xi, eta, weight)` triples (`(t, weight)` for the edge
2//! rule). Triangle weights sum to the reference area `1/2`, quadrilateral
3//! weights to `4`, edge weights to `2`.
4
5const G: f64 = 0.577_350_269_189_625_8; // 1 / sqrt(3)
6
7/// One-point triangle rule (centroid), exact for degree 1.
8pub const TRI3_GAUSS_1: [(f64, f64, f64); 1] = [(1.0 / 3.0, 1.0 / 3.0, 0.5)];
9
10/// Three-point triangle rule (interior points), exact for degree 2.
11pub const TRI3_GAUSS_3: [(f64, f64, f64); 3] = [
12    (1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0),
13    (2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0),
14    (1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0),
15];
16
17/// 2x2 rule on `[-1, 1]^2`, exact for degree 3 per axis. Point order is
18/// `(-,-), (+,-), (+,+), (-,+)`, matching the quad's node order.
19pub const QUAD4_GAUSS_2X2: [(f64, f64, f64); 4] =
20    [(-G, -G, 1.0), (G, -G, 1.0), (G, G, 1.0), (-G, G, 1.0)];
21
22/// Two-point rule on `[-1, 1]` for element edges, exact for degree 3.
23pub const EDGE_GAUSS_2: [(f64, f64); 2] = [(-G, 1.0), (G, 1.0)];
24
25#[cfg(test)]
26mod tests {
27    use super::*;
28
29    fn integrate(rule: &[(f64, f64, f64)], f: impl Fn(f64, f64) -> f64) -> f64 {
30        rule.iter().map(|&(x, y, w)| w * f(x, y)).sum()
31    }
32
33    #[test]
34    fn weights_give_the_reference_areas() {
35        assert!((integrate(&TRI3_GAUSS_1, |_, _| 1.0) - 0.5).abs() < 1e-15);
36        assert!((integrate(&TRI3_GAUSS_3, |_, _| 1.0) - 0.5).abs() < 1e-15);
37        assert!((integrate(&QUAD4_GAUSS_2X2, |_, _| 1.0) - 4.0).abs() < 1e-15);
38        assert!((EDGE_GAUSS_2.iter().map(|p| p.1).sum::<f64>() - 2.0).abs() < 1e-15);
39    }
40
41    #[test]
42    fn triangle_rules_are_exact_to_their_degree() {
43        // Exact: int x^a y^b over the unit triangle = a! b! / (a + b + 2)!.
44        assert!((integrate(&TRI3_GAUSS_1, |x, y| x + 2.0 * y) - 3.0 / 6.0).abs() < 1e-15);
45        for (a, b, exact) in [(2, 0, 1.0 / 12.0), (1, 1, 1.0 / 24.0), (0, 2, 1.0 / 12.0)] {
46            let got = integrate(&TRI3_GAUSS_3, |x, y| x.powi(a) * y.powi(b));
47            assert!((got - exact).abs() < 1e-15, "x^{a} y^{b}: {got} vs {exact}");
48        }
49    }
50
51    #[test]
52    fn quad_and_edge_rules_integrate_cubics_exactly() {
53        // Exact over [-1,1]^2: x^2 y^2 -> 4/9; x^3 y -> 0; x^2 + y^3 + 1 -> 4/3 + 0 + 4.
54        assert!((integrate(&QUAD4_GAUSS_2X2, |x, y| x * x * y * y) - 4.0 / 9.0).abs() < 1e-15);
55        assert!(integrate(&QUAD4_GAUSS_2X2, |x, y| x.powi(3) * y).abs() < 1e-15);
56        assert!(
57            (integrate(&QUAD4_GAUSS_2X2, |x, y| x * x + y.powi(3) + 1.0) - (4.0 / 3.0 + 4.0)).abs()
58                < 1e-14
59        );
60        let edge = |f: fn(f64) -> f64| EDGE_GAUSS_2.iter().map(|&(t, w)| w * f(t)).sum::<f64>();
61        assert!((edge(|t| t * t) - 2.0 / 3.0).abs() < 1e-15);
62        assert!(edge(|t| t.powi(3)).abs() < 1e-15);
63    }
64}
carapace_core/model/continuum/quadrature.rs

carapace_core/model/continuum/
quadrature.rs
