carapace_core/model/continuum/
jacobian.rs

1//! The 2D isoparametric Jacobian: `J = [dx/dxi dy/dxi; dx/deta dy/deta]`
2//! (row `a` is the derivative with respect to natural coordinate `a`), its
3//! determinant and inverse, and physical shape-function derivatives.
4
5/// `J`, `det J` and `J^{-1}` at one point.
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct Jacobian {
8    pub matrix: [[f64; 2]; 2],
9    pub det: f64,
10    /// `J^{-1}`; meaningless (non-finite entries) when `det` is zero.
11    pub inverse: [[f64; 2]; 2],
12}
13
14impl Jacobian {
15    /// `J[a][c] = sum_i dN_i/dxi_a * x_i[c]` for `NN` nodes with reference
16    /// derivatives `dn` (`[dN/dxi, dN/deta]` per node).
17    pub fn new<const NN: usize>(coords: &[[f64; 2]; NN], dn: &[[f64; 2]; NN]) -> Self {
18        let mut j = [[0.0; 2]; 2];
19        for (x, d) in coords.iter().zip(dn) {
20            for a in 0..2 {
21                for c in 0..2 {
22                    j[a][c] += d[a] * x[c];
23                }
24            }
25        }
26        let det = j[0][0] * j[1][1] - j[0][1] * j[1][0];
27        let inverse = [
28            [j[1][1] / det, -j[0][1] / det],
29            [-j[1][0] / det, j[0][0] / det],
30        ];
31        Jacobian {
32            matrix: j,
33            det,
34            inverse,
35        }
36    }
37}
38
39/// `dN_i/dx_c = sum_a J^{-1}[c][a] dN_i/dxi_a`, one `[dN/dx, dN/dy]` per node.
40pub fn physical_derivatives<const NN: usize>(
41    jacobian: &Jacobian,
42    dn: &[[f64; 2]; NN],
43) -> [[f64; 2]; NN] {
44    let inv = &jacobian.inverse;
45    std::array::from_fn(|i| {
46        [
47            inv[0][0] * dn[i][0] + inv[0][1] * dn[i][1],
48            inv[1][0] * dn[i][0] + inv[1][1] * dn[i][1],
49        ]
50    })
51}
52
53#[cfg(test)]
54mod tests {
55    use super::*;
56    use crate::model::continuum::shape::{quad4_dshape, tri3_dshape};
57
58    #[test]
59    fn affine_maps_have_constant_determinant_equal_to_the_area_scale() {
60        // Triangle (0,0), (4,0), (1,3): area 6, reference area 1/2 -> detJ = 12.
61        let tri = [[0.0, 0.0], [4.0, 0.0], [1.0, 3.0]];
62        let j = Jacobian::new(&tri, &tri3_dshape());
63        assert!((j.det - 12.0).abs() < 1e-14);
64        // Parallelogram (0,0), (2,0), (3,1.5), (1,1.5): area 3, reference area 4 -> 0.75 everywhere.
65        let par = [[0.0, 0.0], [2.0, 0.0], [3.0, 1.5], [1.0, 1.5]];
66        for [xi, eta] in [[0.0, 0.0], [0.6, -0.3], [-1.0, 1.0]] {
67            let j = Jacobian::new(&par, &quad4_dshape(xi, eta));
68            assert!((j.det - 0.75).abs() < 1e-14);
69        }
70    }
71
72    #[test]
73    fn physical_derivatives_reproduce_a_linear_field_gradient() {
74        // u(x, y) = 3x - 2y + 1 sampled at the nodes of a skewed quad has gradient (3, -2) everywhere.
75        let quad = [[0.0, 0.0], [2.0, 0.3], [2.4, 1.9], [-0.2, 1.5]];
76        let u: Vec<f64> = quad.iter().map(|p| 3.0 * p[0] - 2.0 * p[1] + 1.0).collect();
77        for [xi, eta] in [[0.0, 0.0], [0.57, -0.57], [-0.8, 0.2]] {
78            let dn = quad4_dshape(xi, eta);
79            let dx = physical_derivatives(&Jacobian::new(&quad, &dn), &dn);
80            let grad: [f64; 2] =
81                std::array::from_fn(|c| dx.iter().zip(&u).map(|(d, u)| d[c] * u).sum());
82            assert!((grad[0] - 3.0).abs() < 1e-13 && (grad[1] + 2.0).abs() < 1e-13);
83        }
84    }
85
86    #[test]
87    fn inverse_times_matrix_is_identity() {
88        let quad = [[0.0, 0.0], [2.0, 0.3], [2.4, 1.9], [-0.2, 1.5]];
89        let j = Jacobian::new(&quad, &quad4_dshape(0.3, -0.4));
90        for a in 0..2 {
91            for b in 0..2 {
92                let v: f64 = (0..2).map(|c| j.matrix[a][c] * j.inverse[c][b]).sum();
93                assert!((v - if a == b { 1.0 } else { 0.0 }).abs() < 1e-14);
94            }
95        }
96    }
97}
carapace_core/model/continuum/jacobian.rs

carapace_core/model/continuum/
jacobian.rs
