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

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