conspire/domain/fem/block/element/linear/pyramid/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    fem::block::element::{
6        FRAC_1_SQRT_3, FiniteElement, ParametricCoordinate, ParametricCoordinates,
7        ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8        linear::{LinearElement, LinearFiniteElement, M},
9    },
10    math::{Scalar, ScalarList},
11};
12
13const G: usize = 8;
14const N: usize = 5;
15const P: usize = N;
16
17pub type Pyramid = LinearElement<G, N>;
18
19impl FiniteElement<G, M, N, P> for Pyramid {
20    fn integration_points() -> ParametricCoordinates<G, M> {
21        integration_points_and_weights().0
22    }
23    fn integration_weights(&self) -> &ScalarList<G> {
24        &self.integration_weights
25    }
26    fn parametric_reference() -> ParametricReference<M, N> {
27        [
28            [-1.0, -1.0, 0.0],
29            [1.0, -1.0, 0.0],
30            [1.0, 1.0, 0.0],
31            [-1.0, 1.0, 0.0],
32            [0.0, 0.0, 1.0],
33        ]
34        .into()
35    }
36    fn parametric_weights() -> ScalarList<G> {
37        integration_points_and_weights().1
38    }
39    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
40        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
41        let bottom = bottom(xi_3);
42        [
43            ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom) / 4.0,
44            ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom) / 4.0,
45            ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom) / 4.0,
46            ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom) / 4.0,
47            xi_3,
48        ]
49        .into()
50    }
51    fn shape_functions_gradients(
52        parametric_coordinate: ParametricCoordinate<M>,
53    ) -> ShapeFunctionsGradients<M, N> {
54        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
55        let bottom = bottom(xi_3);
56        let bottom_squared = bottom * bottom;
57        [
58            [
59                (-(1.0 - xi_2) + xi_2 * xi_3 / bottom) / 4.0,
60                (-(1.0 - xi_1) + xi_1 * xi_3 / bottom) / 4.0,
61                (-1.0 + xi_1 * xi_2 / bottom_squared) / 4.0,
62            ],
63            [
64                ((1.0 - xi_2) - xi_2 * xi_3 / bottom) / 4.0,
65                (-(1.0 + xi_1) - xi_1 * xi_3 / bottom) / 4.0,
66                (-1.0 - xi_1 * xi_2 / bottom_squared) / 4.0,
67            ],
68            [
69                ((1.0 + xi_2) + xi_2 * xi_3 / bottom) / 4.0,
70                ((1.0 + xi_1) + xi_1 * xi_3 / bottom) / 4.0,
71                (-1.0 + xi_1 * xi_2 / bottom_squared) / 4.0,
72            ],
73            [
74                (-(1.0 + xi_2) - xi_2 * xi_3 / bottom) / 4.0,
75                ((1.0 - xi_1) - xi_1 * xi_3 / bottom) / 4.0,
76                (-1.0 - xi_1 * xi_2 / bottom_squared) / 4.0,
77            ],
78            [0.0, 0.0, 1.0],
79        ]
80        .into()
81    }
82}
83
84fn bottom(xi_3: Scalar) -> Scalar {
85    const SMALL: Scalar = 4e1 * f64::EPSILON;
86    if (1.0 - xi_3).abs() > SMALL {
87        1.0 - xi_3
88    } else {
89        SMALL
90    }
91}
92
93fn integration_points_and_weights() -> (ParametricCoordinates<G, M>, ScalarList<G>) {
94    const X: [Scalar; 2] = [0.455_848_155_988_775, 0.877_485_177_344_559];
95    const B: [Scalar; 2] = [0.100_785_882_079_825, 0.232_547_451_253_508];
96    const U1_2D: [Scalar; 4] = [-FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3, -FRAC_1_SQRT_3];
97    const U2_2D: [Scalar; 4] = [-FRAC_1_SQRT_3, -FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3];
98    const W_2D: [Scalar; 4] = [1.0; _];
99    let mut points = [[0.0; M]; G];
100    let mut weights = [0.0; G];
101    let mut i = 0;
102    X.into_iter().zip(B).for_each(|(x, b)| {
103        U1_2D
104            .into_iter()
105            .zip(U2_2D)
106            .zip(W_2D)
107            .for_each(|((u1, u2), w)| {
108                points[i][0] = x * u1;
109                points[i][1] = x * u2;
110                points[i][2] = 1.0 - x;
111                weights[i] = w * b;
112                i += 1;
113            })
114    });
115    (points.into(), weights.into())
116}
117
118impl LinearFiniteElement<G, N> for Pyramid {}
conspire/domain/fem/block/element/linear/pyramid/mod.rs

conspire/domain/fem/block/element/linear/pyramid/
mod.rs
