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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    fem::block::element::{
6        ElementNodalEitherCoordinates, FRAC_1_SQRT_3, FiniteElement, ParametricCoordinate,
7        ParametricCoordinates, 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 scaled_jacobians<const I: usize>(
40        _nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
41    ) -> ScalarList<P> {
42        todo!()
43    }
44    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
45        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
46        let bottom = bottom(xi_3);
47        [
48            ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom) / 4.0,
49            ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom) / 4.0,
50            ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom) / 4.0,
51            ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom) / 4.0,
52            xi_3,
53        ]
54        .into()
55    }
56    fn shape_functions_gradients(
57        parametric_coordinate: ParametricCoordinate<M>,
58    ) -> ShapeFunctionsGradients<M, N> {
59        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
60        let bottom = bottom(xi_3);
61        let bottom_squared = bottom * bottom;
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            [
79                (-(1.0 + xi_2) - xi_2 * xi_3 / bottom) / 4.0,
80                ((1.0 - xi_1) - xi_1 * xi_3 / bottom) / 4.0,
81                (-1.0 - xi_1 * xi_2 / bottom_squared) / 4.0,
82            ],
83            [0.0, 0.0, 1.0],
84        ]
85        .into()
86    }
87}
88
89fn bottom(xi_3: Scalar) -> Scalar {
90    const SMALL: Scalar = 4e1 * f64::EPSILON;
91    if (1.0 - xi_3).abs() > SMALL {
92        1.0 - xi_3
93    } else {
94        SMALL
95    }
96}
97
98fn integration_points_and_weights() -> (ParametricCoordinates<G, M>, ScalarList<G>) {
99    const X: [Scalar; 2] = [0.455_848_155_988_775, 0.877_485_177_344_559];
100    const B: [Scalar; 2] = [0.100_785_882_079_825, 0.232_547_451_253_508];
101    const U1_2D: [Scalar; 4] = [-FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3, -FRAC_1_SQRT_3];
102    const U2_2D: [Scalar; 4] = [-FRAC_1_SQRT_3, -FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3];
103    const W_2D: [Scalar; 4] = [1.0; _];
104    let mut points = [[0.0; M]; G];
105    let mut weights = [0.0; G];
106    let mut i = 0;
107    X.into_iter().zip(B).for_each(|(x, b)| {
108        U1_2D
109            .into_iter()
110            .zip(U2_2D)
111            .zip(W_2D)
112            .for_each(|((u1, u2), w)| {
113                points[i][0] = x * u1;
114                points[i][1] = x * u2;
115                points[i][2] = 1.0 - x;
116                weights[i] = w * b;
117                i += 1;
118            })
119    });
120    (points.into(), weights.into())
121}
122
123impl LinearFiniteElement<G, N> for Pyramid {}
conspire/domain/fem/block/element/linear/pyramid/mod.rs

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