carapace_core/model/continuum/
shape.rs

1//! Shape functions and reference derivatives of the 3-node triangle (natural
2//! coordinates `xi, eta` on the unit triangle) and the 4-node quadrilateral
3//! (`xi, eta` in `[-1, 1]`). Node order is counter-clockwise; quad node `k`
4//! sits at the reference corner `QUAD4_CORNERS[k]`.
5
6/// Reference corners of the bilinear quadrilateral, in node order.
7pub const QUAD4_CORNERS: [[f64; 2]; 4] = [[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]];
8
9/// `N_i(xi, eta)` for the triangle: `1 - xi - eta`, `xi`, `eta`.
10pub fn tri3_shape(xi: f64, eta: f64) -> [f64; 3] {
11    [1.0 - xi - eta, xi, eta]
12}
13
14/// `[dN_i/dxi, dN_i/deta]` for the triangle (constant).
15pub fn tri3_dshape() -> [[f64; 2]; 3] {
16    [[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]]
17}
18
19/// `N_i(xi, eta) = (1 + xi_i xi)(1 + eta_i eta) / 4`.
20pub fn quad4_shape(xi: f64, eta: f64) -> [f64; 4] {
21    std::array::from_fn(|i| {
22        let [xi_i, eta_i] = QUAD4_CORNERS[i];
23        0.25 * (1.0 + xi_i * xi) * (1.0 + eta_i * eta)
24    })
25}
26
27/// `[dN_i/dxi, dN_i/deta]` for the quadrilateral.
28pub fn quad4_dshape(xi: f64, eta: f64) -> [[f64; 2]; 4] {
29    std::array::from_fn(|i| {
30        let [xi_i, eta_i] = QUAD4_CORNERS[i];
31        [
32            0.25 * xi_i * (1.0 + eta_i * eta),
33            0.25 * eta_i * (1.0 + xi_i * xi),
34        ]
35    })
36}
37
38#[cfg(test)]
39mod tests {
40    use super::*;
41
42    const POINTS: [[f64; 2]; 4] = [[0.1, 0.2], [0.3, 0.3], [-0.7, 0.4], [0.5, -0.9]];
43
44    #[test]
45    fn shape_functions_partition_unity_and_derivatives_sum_to_zero() {
46        for [xi, eta] in POINTS {
47            let tri = [xi.abs() * 0.4, eta.abs() * 0.4];
48            assert!((tri3_shape(tri[0], tri[1]).iter().sum::<f64>() - 1.0).abs() < 1e-15);
49            assert!((quad4_shape(xi, eta).iter().sum::<f64>() - 1.0).abs() < 1e-15);
50            for axis in 0..2 {
51                assert!(tri3_dshape().iter().map(|d| d[axis]).sum::<f64>().abs() < 1e-15);
52                assert!(
53                    quad4_dshape(xi, eta)
54                        .iter()
55                        .map(|d| d[axis])
56                        .sum::<f64>()
57                        .abs()
58                        < 1e-15
59                );
60            }
61        }
62    }
63
64    #[test]
65    fn shape_functions_are_one_at_their_own_node_and_zero_elsewhere() {
66        let tri_nodes = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]];
67        for (k, [xi, eta]) in tri_nodes.into_iter().enumerate() {
68            for (i, n) in tri3_shape(xi, eta).into_iter().enumerate() {
69                assert_eq!(n, if i == k { 1.0 } else { 0.0 });
70            }
71        }
72        for (k, [xi, eta]) in QUAD4_CORNERS.into_iter().enumerate() {
73            for (i, n) in quad4_shape(xi, eta).into_iter().enumerate() {
74                assert_eq!(n, if i == k { 1.0 } else { 0.0 });
75            }
76        }
77    }
78
79    #[test]
80    fn quad_derivatives_match_central_finite_differences() {
81        let h = 1e-6;
82        for [xi, eta] in POINTS {
83            let d = quad4_dshape(xi, eta);
84            let (dxi_p, dxi_m) = (quad4_shape(xi + h, eta), quad4_shape(xi - h, eta));
85            let (deta_p, deta_m) = (quad4_shape(xi, eta + h), quad4_shape(xi, eta - h));
86            for i in 0..4 {
87                assert!((d[i][0] - (dxi_p[i] - dxi_m[i]) / (2.0 * h)).abs() < 1e-9);
88                assert!((d[i][1] - (deta_p[i] - deta_m[i]) / (2.0 * h)).abs() < 1e-9);
89            }
90        }
91        let d = tri3_dshape();
92        let (a, b) = (tri3_shape(0.3 + h, 0.2), tri3_shape(0.3 - h, 0.2));
93        let (c, e) = (tri3_shape(0.3, 0.2 + h), tri3_shape(0.3, 0.2 - h));
94        for i in 0..3 {
95            assert!((d[i][0] - (a[i] - b[i]) / (2.0 * h)).abs() < 1e-9);
96            assert!((d[i][1] - (c[i] - e[i]) / (2.0 * h)).abs() < 1e-9);
97        }
98    }
99}
carapace_core/model/continuum/shape.rs

carapace_core/model/continuum/
shape.rs
