mesh_geometry/utils/
jacobian.rs

1use crate::{Float, Point2};
2
3/// 2×2 Jacobian matrix for a bilinear quad with nodes [A,B,C,D] at
4/// reference coords (-1,-1), (1,-1), (1,1), (-1,1).
5#[derive(Debug, Clone, Copy, PartialEq)]
6pub struct Jacobian2x2<T: Float> {
7    /// Entry (1,1) of the Jacobian matrix
8    pub m11: T,
9    /// Entry (1,2) of the Jacobian matrix
10    pub m12: T,
11    /// Entry (2,1) of the Jacobian matrix
12    pub m21: T,
13    /// Entry (2,2) of the Jacobian matrix
14    pub m22: T,
15}
16
17impl<T: Float> Jacobian2x2<T> {
18    /// Compute J(xi,eta) = ∑ Ni,xi * Xi for i=1..4
19    pub fn for_quad(
20        xi: T,
21        eta: T,
22        a: Point2<T>,
23        b: Point2<T>,
24        c: Point2<T>,
25        d: Point2<T>,
26    ) -> Self {
27        // shape‐fn derivatives
28        let one = T::one();
29        let d_n1_dxi = -one + eta; let d_n1_deta = -one + xi; // node A
30        let d_n2_dxi =  one - eta; let d_n2_deta = -one - xi; // node B
31        let d_n3_dxi =  one + eta; let d_n3_deta =  one + xi; // node C
32        let d_n4_dxi = -one - eta; let d_n4_deta =  one - xi; // node D
33        // physical grads
34        let m11 = d_n1_dxi * a.x + d_n2_dxi * b.x + d_n3_dxi * c.x + d_n4_dxi * d.x;
35        let m12 = d_n1_deta * a.x + d_n2_deta * b.x + d_n3_deta * c.x + d_n4_deta * d.x;
36        let m21 = d_n1_dxi * a.y + d_n2_dxi * b.y + d_n3_dxi * c.y + d_n4_dxi * d.y;
37        let m22 = d_n1_deta * a.y + d_n2_deta * b.y + d_n3_deta * c.y + d_n4_deta * d.y;
38        Jacobian2x2 { m11, m12, m21, m22 }
39    }
40
41    /// Determinant det(J)
42    pub fn det(self) -> T {
43        self.m11 * self.m22 - self.m12 * self.m21
44    }
45
46    /// Inverse J⁻¹
47    pub fn inverse(self) -> Option<Jacobian2x2<T>> {
48        let d = self.det();
49        if d.abs() < T::epsilon() { return None; }
50        let inv = T::one() / d;
51        Some(Jacobian2x2 {
52            m11:  self.m22 * inv,
53            m12: -self.m12 * inv,
54            m21: -self.m21 * inv,
55            m22:  self.m11 * inv,
56        })
57    }
58}
59
60/// Given physical point `p` and quad corners `a,b,c,d`, find (xi,eta) via Newton:
61pub fn invert_quad_mapping<T: Float>(
62    mut xi: T,
63    mut eta: T,
64    p: Point2<T>,
65    a: Point2<T>,
66    b: Point2<T>,
67    c: Point2<T>,
68    d: Point2<T>,
69    tol: T,
70    max_iters: usize,
71) -> Option<(T, T)> {
72    for _ in 0..max_iters {
73        // shape‐fn values
74        let n1 = (T::one()-xi)*(T::one()-eta)/T::from(4.0).unwrap();
75        let n2 = (T::one()+xi)*(T::one()-eta)/T::from(4.0).unwrap();
76        let n3 = (T::one()+xi)*(T::one()+eta)/T::from(4.0).unwrap();
77        let n4 = (T::one()-xi)*(T::one()+eta)/T::from(4.0).unwrap();
78        // physical point at (xi,eta)
79        let x = a.x*n1 + b.x*n2 + c.x*n3 + d.x*n4;
80        let y = a.y*n1 + b.y*n2 + c.y*n3 + d.y*n4;
81        // residual
82        let rx = x - p.x;
83        let ry = y - p.y;
84        if rx.abs() < tol && ry.abs() < tol {
85            return Some((xi, eta));
86        }
87        // Jacobian & inverse
88        let j = Jacobian2x2::for_quad(xi, eta, a, b, c, d);
89        let inv_j = match j.inverse() {
90            Some(inv) => inv,
91            None => return None,
92        };
93        // Newton update: [Δxi; Δeta] = invJ * [rx; ry]
94        let dxi = inv_j.m11 * rx + inv_j.m12 * ry;
95        let deta = inv_j.m21 * rx + inv_j.m22 * ry;
96        xi = xi - dxi;
97        eta = eta - deta;
98    }
99    None
100}
101
102#[cfg(test)]
103mod tests {
104    use super::*;
105    use crate::Point2;
106    use num_traits::Zero;
107
108    #[test]
109    fn jacobian_quad_det_nonzero() {
110        let a = Point2::new(0.0,0.0);
111        let b = Point2::new(2.0,0.0);
112        let c = Point2::new(2.0,1.0);
113        let d = Point2::new(0.0,1.0);
114        let j = Jacobian2x2::for_quad(0.0_f64, 0.0, a,b,c,d );
115        assert!(!j.det().is_zero());
116    }
117
118    #[test]
119    fn invert_quad_identity() {
120        // unit square maps xi,eta -> x,y one-to-one
121        let a=Point2::new(0.0,0.0);
122        let b=Point2::new(1.0,0.0);
123        let c=Point2::new(1.0,1.0);
124        let d=Point2::new(0.0,1.0);
125        let p = Point2::new(0.3, 0.7);
126        let got = invert_quad_mapping(0.0_f64, 0.0, p, a,b,c,d, 1e-6, 50 );
127        assert!(got.is_some());
128        let (xi,eta)=got.unwrap();
129        assert!((xi - (-0.4)).abs() < 1e-3);
130        assert!((eta - 0.4).abs() < 1e-3);
131    }
132}
mesh_geometry/utils/jacobian.rs

mesh_geometry/utils/
jacobian.rs
