gemlab/integ/
analytical_tri3.rs

1use crate::shapes::{GeoKind, Scratchpad};
2use crate::StrError;
3use russell_lab::math::SQRT_2;
4use russell_lab::{mat_mat_mul, mat_t_mat_mul, Matrix, Vector};
5use russell_tensor::{LinElasticity, Tensor2};
6
7/// Performs analytical integrations on a Tri3
8pub struct AnalyticalTri3 {
9    /// Holds the area of the triangle
10    pub area: f64,
11
12    /// Holds the gradients (B-matrix)
13    ///
14    /// ```text
15    ///             →
16    /// →  →    dNᵐ(ξ)
17    /// Bᵐ(ξ) = ——————
18    ///            →
19    ///           dx
20    /// ```
21    ///
22    /// Organized as the B matrix (nnode=3, space_ndim=2)
23    pub bb: Matrix,
24
25    /// Holds the (element) Be-matrix (4, 6)
26    pub bbe: Matrix,
27
28    // Holds the lengths of each edge
29    pub ll: Vec<f64>,
30
31    // Coordinates
32    x0: f64,
33    x1: f64,
34    x2: f64,
35}
36
37impl AnalyticalTri3 {
38    /// Allocates a new instance
39    pub fn new(pad: &Scratchpad) -> Self {
40        assert_eq!(pad.kind, GeoKind::Tri3);
41
42        // coefficients
43        let (x0, y0) = (pad.xxt.get(0, 0), pad.xxt.get(1, 0));
44        let (x1, y1) = (pad.xxt.get(0, 1), pad.xxt.get(1, 1));
45        let (x2, y2) = (pad.xxt.get(0, 2), pad.xxt.get(1, 2));
46        let (a0, a1, a2) = (y1 - y2, y2 - y0, y0 - y1);
47        let (b0, b1, b2) = (x2 - x1, x0 - x2, x1 - x0);
48        let (f0, f1, f2) = (x1 * y2 - x2 * y1, x2 * y0 - x0 * y2, x0 * y1 - x1 * y0);
49
50        // area
51        let area = (f0 + f1 + f2) / 2.0;
52
53        // auxiliary
54        let r = 2.0 * area;
55        let s = r * SQRT_2;
56
57        // gradients
58        #[rustfmt::skip]
59        let bb = Matrix::from(&[
60            [a0/r, b0/r],
61            [a1/r, b1/r],
62            [a2/r, b2/r],
63        ]);
64
65        // element Be-matrix (plane-strain or plane-stress)
66        #[rustfmt::skip]
67        let bbe = Matrix::from(&[
68            [a0/r,  0.0, a1/r,  0.0, a2/r,  0.0],
69            [ 0.0, b0/r,  0.0, b1/r,  0.0, b2/r],
70            [ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0],
71            [b0/s, a0/s, b1/s, a1/s, b2/s, a2/s],
72        ]);
73
74        // edges lengths
75        let ll = vec![
76            f64::sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1)), // L01
77            f64::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)), // L12
78            f64::sqrt((x2 - x0) * (x2 - x0) + (y2 - y0) * (y2 - y0)), // L20
79        ];
80
81        // results
82        AnalyticalTri3 {
83            area,
84            bb,
85            bbe,
86            ll,
87            x0,
88            x1,
89            x2,
90        }
91    }
92
93    /// Integrates shape times scalar with constant function s(x)
94    pub fn vec_01_ns(&self, s: f64, axisymmetric: bool) -> Vector {
95        if axisymmetric {
96            let c = s * self.area / 12.0;
97            Vector::from(&[
98                c * (2.0 * self.x0 + self.x1 + self.x2),
99                c * (self.x0 + 2.0 * self.x1 + self.x2),
100                c * (self.x0 + self.x1 + 2.0 * self.x2),
101            ])
102        } else {
103            let c = s * self.area / 3.0;
104            Vector::from(&[c, c, c])
105        }
106    }
107
108    /// Integrates shape times scalar with constant function s(x)
109    ///
110    /// **Important:** `side` must be 0, 1, or 2
111    pub fn vec_01_ns_bry(&self, side: usize, s: f64, axisymmetric: bool) -> Vector {
112        if axisymmetric {
113            let (a0, a1) = match side {
114                0 => (self.x1, self.x0),
115                1 => (self.x2, self.x1),
116                2 => (self.x0, self.x2),
117                _ => panic!("invalid side"),
118            };
119            let c = s * self.ll[side] / 6.0;
120            Vector::from(&[c * (2.0 * a0 + a1), c * (a0 + 2.0 * a1)])
121        } else {
122            let c = s * self.ll[side] / 2.0;
123            Vector::from(&[c, c])
124        }
125    }
126
127    /// Integrates shape times vector with constant vector v(x) = {v₀, v₁}
128    ///
129    /// solution:
130    ///
131    /// ```text
132    /// bᵐ₀ = v₀ A / 3
133    /// bᵐ₁ = v₁ A / 3
134    ///       ┌    ┐
135    ///       │ v0 │
136    ///       │ v1 │
137    ///     A │ v0 │
138    /// b = — │ v1 │
139    ///     3 │ v0 │
140    ///       │ v1 │
141    ///       └    ┘
142    /// ```
143    pub fn vec_02_nv(&self, v0: f64, v1: f64) -> Vector {
144        Vector::from(&[
145            v0 * self.area / 3.0,
146            v1 * self.area / 3.0,
147            v0 * self.area / 3.0,
148            v1 * self.area / 3.0,
149            v0 * self.area / 3.0,
150            v1 * self.area / 3.0,
151        ])
152    }
153
154    /// Integrates gradient dot vector with constant vector function w(x) = {w₀, w₁}
155    ///
156    /// solution:
157    ///
158    /// ```text
159    /// cᵐ = (w₀ Bᵐ₀ + w₁ Bᵐ₁) A
160    /// ```
161    pub fn vec_03_bv(&self, w0: f64, w1: f64) -> Vector {
162        Vector::from(&[
163            (w0 * self.bb.get(0, 0) + w1 * self.bb.get(0, 1)) * self.area,
164            (w0 * self.bb.get(1, 0) + w1 * self.bb.get(1, 1)) * self.area,
165            (w0 * self.bb.get(2, 0) + w1 * self.bb.get(2, 1)) * self.area,
166        ])
167    }
168
169    /// Integrates gradient dot vector with bilinear vector function w(x) = {x, y}
170    ///
171    /// solution:
172    ///
173    /// ```text
174    /// cᵐ = ((x₀+x₁+x₂) Bᵐ₀ + (y₀+y₁+y₂) Bᵐ₁) A / 3
175    /// ```
176    ///
177    /// # Input
178    ///
179    /// * `pad` -- The same shape used in `new` because we need the nodal coordinates here
180    ///            Do not change the coordinates, otherwise the values will be wrong.
181    #[rustfmt::skip]
182    pub fn vec_03_bv_bilinear(&self, pad: &Scratchpad) -> Vector {
183        let (x0, x1, x2) = (pad.xxt.get(0,0), pad.xxt.get(0,1), pad.xxt.get(0,2));
184        let (y0, y1, y2) = (pad.xxt.get(1,0), pad.xxt.get(1,1), pad.xxt.get(1,2));
185        Vector::from(&[
186            ((x0 + x1 + x2) * self.bb.get(0,0) + (y0 + y1 + y2) * self.bb.get(0,1)) * self.area / 3.0,
187            ((x0 + x1 + x2) * self.bb.get(1,0) + (y0 + y1 + y2) * self.bb.get(1,1)) * self.area / 3.0,
188            ((x0 + x1 + x2) * self.bb.get(2,0) + (y0 + y1 + y2) * self.bb.get(2,1)) * self.area / 3.0,
189        ])
190    }
191
192    /// Integrates gradient dot transpose tensor with constant tensor function σ(x) = {σ₀₀, σ₁₁, σ₂₂, σ₀₁√2}
193    #[rustfmt::skip]
194    pub fn vec_04_bt(&self, tt: &Tensor2, axisymmetric: bool) -> Vector {
195        let (x0, x1, x2) = (self.x0, self.x1, self.x2);
196        let (b00, b01) = (self.bb.get(0,0), self.bb.get(0,1));
197        let (b10, b11) = (self.bb.get(1,0), self.bb.get(1,1));
198        let (b20, b21) = (self.bb.get(2,0), self.bb.get(2,1));
199        let (t0, t1, t2, t3) = (tt.vector()[0], tt.vector()[1], tt.vector()[2], tt.vector()[3]);
200        if axisymmetric {
201            let c = self.area / 3.0;
202            Vector::from(&[
203                c * (b00 * t0 + b01 * t3 / SQRT_2) * (x0 + x1 + x2) + c * t2,
204                c * (b01 * t1 + b00 * t3 / SQRT_2) * (x0 + x1 + x2),
205                c * (b10 * t0 + b11 * t3 / SQRT_2) * (x0 + x1 + x2) + c * t2,
206                c * (b11 * t1 + b10 * t3 / SQRT_2) * (x0 + x1 + x2),
207                c * (b20 * t0 + b21 * t3 / SQRT_2) * (x0 + x1 + x2) + c * t2,
208                c * (b21 * t1 + b20 * t3 / SQRT_2) * (x0 + x1 + x2),
209            ])
210        } else {
211            let c = self.area;
212            Vector::from(&[
213                c * (b00 * t0 + b01 * t3 / SQRT_2),
214                c * (b01 * t1 + b00 * t3 / SQRT_2),
215                c * (b10 * t0 + b11 * t3 / SQRT_2),
216                c * (b11 * t1 + b10 * t3 / SQRT_2),
217                c * (b20 * t0 + b21 * t3 / SQRT_2),
218                c * (b21 * t1 + b20 * t3 / SQRT_2),
219            ])
220        }
221    }
222
223    /// Performs the n-s-n integration with constant s(x) field
224    pub fn mat_01_nsn(&self, s: f64, th: f64) -> Matrix {
225        let c = th * s * self.area / 12.0;
226        Matrix::from(&[[2.0 * c, c, c], [c, 2.0 * c, c], [c, c, 2.0 * c]])
227    }
228
229    /// Performs the n-s-n integration with constant s(x) field (boundary integral version)
230    ///
231    /// **Important:** `side` must be 0, 1, or 2
232    pub fn mat_01_nsn_bry(&self, side: usize, s: f64, axisymmetric: bool) -> Matrix {
233        if axisymmetric {
234            let (a0, a1) = match side {
235                0 => (self.x1, self.x0),
236                1 => (self.x2, self.x1),
237                2 => (self.x0, self.x2),
238                _ => panic!("invalid side"),
239            };
240            let c = s * self.ll[side] / 12.0;
241            Matrix::from(&[
242                [c * (3.0 * a0 + a1), c * (a0 + a1)],
243                [c * (a0 + a1), c * (a0 + 3.0 * a1)],
244            ])
245        } else {
246            let c = s * self.ll[side] / 6.0;
247            Matrix::from(&[[2.0 * c, c], [c, 2.0 * c]])
248        }
249    }
250
251    /// Performs the b-v-n integration with constant vector
252    #[rustfmt::skip]
253    pub fn mat_02_bvn(&self, v0: f64, v1: f64) -> Matrix {
254        let aa = self.area;
255        let g = &self.bb;
256        Matrix::from(&[
257            [aa*(g.get(0,0)*v0 + g.get(0,1)*v1)/3.0, aa*(g.get(0,0)*v0 + g.get(0,1)*v1)/3.0, aa*(g.get(0,0)*v0 + g.get(0,1)*v1)/3.0],
258            [aa*(g.get(1,0)*v0 + g.get(1,1)*v1)/3.0, aa*(g.get(1,0)*v0 + g.get(1,1)*v1)/3.0, aa*(g.get(1,0)*v0 + g.get(1,1)*v1)/3.0],
259            [aa*(g.get(2,0)*v0 + g.get(2,1)*v1)/3.0, aa*(g.get(2,0)*v0 + g.get(2,1)*v1)/3.0, aa*(g.get(2,0)*v0 + g.get(2,1)*v1)/3.0],
260        ])
261    }
262
263    /// Performs the b-v-n integration with v = {x, y}
264    #[rustfmt::skip]
265    pub fn mat_02_bvn_bilinear(&self, pad: &Scratchpad) -> Matrix {
266        let aa = self.area;
267        let (b00, b01) = (self.bb.get(0,0), self.bb.get(0,1));
268        let (b10, b11) = (self.bb.get(1,0), self.bb.get(1,1));
269        let (b20, b21) = (self.bb.get(2,0), self.bb.get(2,1));
270        let (x0, y0) = (pad.xxt.get(0,0), pad.xxt.get(1,0));
271        let (x1, y1) = (pad.xxt.get(0,1), pad.xxt.get(1,1));
272        let (x2, y2) = (pad.xxt.get(0,2), pad.xxt.get(1,2));
273        Matrix::from(&[
274            [(aa*(2.0*b00*x0 + b00*x1 + b00*x2 + 2.0*b01*y0 + b01*y1 + b01*y2))/12.0,(aa*(b00*x0 + 2.0*b00*x1 + b00*x2 + b01*y0 + 2.0*b01*y1 + b01*y2))/12.0, (aa*(b00*x0 + b00*x1 + 2.0*b00*x2 + b01*y0 + b01*y1 + 2.0*b01*y2))/12.0],
275            [(aa*(2.0*b10*x0 + b10*x1 + b10*x2 + 2.0*b11*y0 + b11*y1 + b11*y2))/12.0,(aa*(b10*x0 + 2.0*b10*x1 + b10*x2 + b11*y0 + 2.0*b11*y1 + b11*y2))/12.0, (aa*(b10*x0 + b10*x1 + 2.0*b10*x2 + b11*y0 + b11*y1 + 2.0*b11*y2))/12.0],
276            [(aa*(2.0*b20*x0 + b20*x1 + b20*x2 + 2.0*b21*y0 + b21*y1 + b21*y2))/12.0,(aa*(b20*x0 + 2.0*b20*x1 + b20*x2 + b21*y0 + 2.0*b21*y1 + b21*y2))/12.0, (aa*(b20*x0 + b20*x1 + 2.0*b20*x2 + b21*y0 + b21*y1 + 2.0*b21*y2))/12.0],
277        ])
278    }
279
280    /// Performs the b-t-b integration with constant tensor
281    pub fn mat_03_btb(&self, kx: f64, ky: f64, axisymmetric: bool) -> Matrix {
282        let c = if axisymmetric {
283            self.area * (self.x0 + self.x1 + self.x2) / 3.0
284        } else {
285            self.area
286        };
287        let (b00, b01) = (self.bb.get(0, 0), self.bb.get(0, 1));
288        let (b10, b11) = (self.bb.get(1, 0), self.bb.get(1, 1));
289        let (b20, b21) = (self.bb.get(2, 0), self.bb.get(2, 1));
290        let k00 = c * (b00 * b00 * kx + b01 * b01 * ky);
291        let k11 = c * (b10 * b10 * kx + b11 * b11 * ky);
292        let k01 = c * (b00 * b10 * kx + b01 * b11 * ky);
293        let k12 = c * (b10 * b20 * kx + b11 * b21 * ky);
294        let k02 = c * (b00 * b20 * kx + b01 * b21 * ky);
295        let k22 = c * (b20 * b20 * kx + b21 * b21 * ky);
296        Matrix::from(&[[k00, k01, k02], [k01, k11, k12], [k02, k12, k22]])
297    }
298
299    /// Performs the n-t-n integration with constant tensor field rho*delta
300    pub fn mat_08_ntn(&self, rho: f64, th: f64) -> Matrix {
301        let c = th * rho * self.area / 12.0;
302        Matrix::from(&[
303            [c * 2.0, c * 0.0, c * 1.0, c * 0.0, c * 1.0, c * 0.0],
304            [c * 0.0, c * 2.0, c * 0.0, c * 1.0, c * 0.0, c * 1.0],
305            [c * 1.0, c * 0.0, c * 2.0, c * 0.0, c * 1.0, c * 0.0],
306            [c * 0.0, c * 1.0, c * 0.0, c * 2.0, c * 0.0, c * 1.0],
307            [c * 1.0, c * 0.0, c * 1.0, c * 0.0, c * 2.0, c * 0.0],
308            [c * 0.0, c * 1.0, c * 0.0, c * 1.0, c * 0.0, c * 2.0],
309        ])
310    }
311
312    /// Performs the n-v-b integration with constant vector
313    #[rustfmt::skip]
314    pub fn mat_09_nvb(&self, v0: f64, v1: f64) -> Matrix {
315        let (b00, b01) = (self.bb.get(0,0), self.bb.get(0,1));
316        let (b10, b11) = (self.bb.get(1,0), self.bb.get(1,1));
317        let (b20, b21) = (self.bb.get(2,0), self.bb.get(2,1));
318        let c = self.area / 3.0;
319        Matrix::from(&[
320            [c*b00*v0, c*b01*v0, c*b10*v0, c*b11*v0, c*b20*v0, c*b21*v0],
321            [c*b00*v1, c*b01*v1, c*b10*v1, c*b11*v1, c*b20*v1, c*b21*v1],
322            [c*b00*v0, c*b01*v0, c*b10*v0, c*b11*v0, c*b20*v0, c*b21*v0],
323            [c*b00*v1, c*b01*v1, c*b10*v1, c*b11*v1, c*b20*v1, c*b21*v1],
324            [c*b00*v0, c*b01*v0, c*b10*v0, c*b11*v0, c*b20*v0, c*b21*v0],
325            [c*b00*v1, c*b01*v1, c*b10*v1, c*b11*v1, c*b20*v1, c*b21*v1],
326        ])
327    }
328
329    /// Performs the b-d-b integration with constant tensor field (calculates the stiffness matrix)
330    ///
331    /// solution:
332    ///
333    /// ```text
334    /// K = Bᵀ ⋅ D ⋅ B ⋅ th ⋅ area
335    /// ```
336    pub fn mat_10_bdb(&self, young: f64, poisson: f64, plane_stress: bool, th: f64) -> Result<Matrix, StrError> {
337        let ela = LinElasticity::new(young, poisson, true, plane_stress);
338        let dd = ela.get_modulus();
339        let dim_dd = 4;
340        let dim_kk = 6;
341        let mut bb_t_dd = Matrix::new(dim_kk, dim_dd);
342        let mut kk = Matrix::new(dim_kk, dim_kk);
343        mat_t_mat_mul(&mut bb_t_dd, 1.0, &self.bbe, dd.matrix(), 0.0).unwrap();
344        mat_mat_mul(&mut kk, th * self.area, &bb_t_dd, &self.bbe, 0.0).unwrap();
345        Ok(kk)
346    }
347}
348
349////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
350
351#[cfg(test)]
352mod tests {
353    use super::AnalyticalTri3;
354    use crate::shapes::{GeoKind, Scratchpad};
355    use russell_lab::approx_eq;
356
357    #[test]
358    fn new_works() {
359        // Tri3 # 1 from Figure 1.18, page 29 of @bhatti
360        // @bhatti Bhatti, M.A. (2005) Fundamental Finite Element Analysis and Applications, Wiley, 700p.
361        let space_ndim = 2;
362        let mut pad = Scratchpad::new(space_ndim, GeoKind::Tri3).unwrap();
363        pad.set_xx(0, 0, 0.0);
364        pad.set_xx(0, 1, 0.0);
365        pad.set_xx(1, 0, 0.2);
366        pad.set_xx(1, 1, 0.0);
367        pad.set_xx(2, 0, 0.1);
368        pad.set_xx(2, 1, 0.1);
369        let ana = AnalyticalTri3::new(&pad);
370        approx_eq(ana.area, 0.01, 1e-15);
371        // a = [-0.1, 0.1, 0.0]
372        // b = [-0.1, -0.1, 0.2]
373        // A = 0.01, Gmi = ai/(2 A) = -0.1 / 0.02 = -5
374        assert_eq!(
375            format!("{:.2}", ana.bb),
376            "┌             ┐\n\
377             │ -5.00 -5.00 │\n\
378             │  5.00 -5.00 │\n\
379             │  0.00 10.00 │\n\
380             └             ┘"
381        );
382    }
383}
gemlab/integ/analytical_tri3.rs

gemlab/integ/
analytical_tri3.rs
