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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    fem::block::element::{
6        ElementNodalEitherCoordinates, FRAC_SQRT_3_5, FiniteElement, ParametricCoordinate,
7        ParametricCoordinates, ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8        quadratic::{M, QuadraticElement, QuadraticFiniteElement},
9    },
10    math::{Scalar, ScalarList, TensorRank1},
11};
12
13const G: usize = 27;
14const N: usize = 13;
15const P: usize = N;
16
17pub type Pyramid = QuadraticElement<G, N>;
18
19impl FiniteElement<G, M, N, P> for Pyramid {
20    fn integration_points() -> ParametricCoordinates<G, M> {
21        const X: [f64; 3] = [0.294997790111502, 0.652996233961648, 0.927005975926850];
22        let u1_2d = [
23            -FRAC_SQRT_3_5,
24            0.0,
25            FRAC_SQRT_3_5,
26            -FRAC_SQRT_3_5,
27            0.0,
28            FRAC_SQRT_3_5,
29            -FRAC_SQRT_3_5,
30            0.0,
31            FRAC_SQRT_3_5,
32        ];
33        let u2_2d = [
34            -FRAC_SQRT_3_5,
35            -FRAC_SQRT_3_5,
36            -FRAC_SQRT_3_5,
37            0.0,
38            0.0,
39            0.0,
40            FRAC_SQRT_3_5,
41            FRAC_SQRT_3_5,
42            FRAC_SQRT_3_5,
43        ];
44        X.into_iter()
45            .flat_map(|x| {
46                u1_2d
47                    .into_iter()
48                    .zip(u2_2d)
49                    .map(move |(u1, u2)| TensorRank1::from([x * u1, x * u2, 1.0 - x]))
50            })
51            .collect()
52    }
53    fn integration_weights(&self) -> &ScalarList<G> {
54        &self.integration_weights
55    }
56    fn parametric_reference() -> ParametricReference<M, N> {
57        [
58            [-1.0, -1.0, 0.0],
59            [1.0, -1.0, 0.0],
60            [1.0, 1.0, 0.0],
61            [-1.0, 1.0, 0.0],
62            [0.0, 0.0, 1.0],
63            [0.0, -1.0, 0.0],
64            [1.0, 0.0, 0.0],
65            [0.0, 1.0, 0.0],
66            [-1.0, 0.0, 0.0],
67            [-0.5, -0.5, 0.5],
68            [0.5, -0.5, 0.5],
69            [0.5, 0.5, 0.5],
70            [-0.5, 0.5, 0.5],
71        ]
72        .into()
73    }
74    fn parametric_weights() -> ScalarList<G> {
75        const B: [f64; 3] = [0.029950703008581, 0.146246269259866, 0.157136361064887];
76        const W1: f64 = 5.0 / 9.0;
77        const W2: f64 = 8.0 / 9.0;
78        let w_2d = [
79            W1 * W1,
80            W2 * W1,
81            W1 * W1,
82            W1 * W2,
83            W2 * W2,
84            W1 * W2,
85            W1 * W1,
86            W2 * W1,
87            W1 * W1,
88        ];
89        B.into_iter()
90            .flat_map(|b| w_2d.into_iter().map(move |w| w * b))
91            .collect()
92    }
93    fn scaled_jacobians<const I: usize>(
94        _nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
95    ) -> ScalarList<P> {
96        todo!()
97    }
98    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
99        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
100        let bottom = bottom(xi_3);
101        [
102            0.25 * (-xi_1 - xi_2 - 1.0)
103                * ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
104            0.25 * (xi_1 - xi_2 - 1.0)
105                * ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
106            0.25 * (xi_1 + xi_2 - 1.0)
107                * ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom),
108            0.25 * (-xi_1 + xi_2 - 1.0)
109                * ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom),
110            xi_3 * (2.0 * xi_3 - 1.0),
111            0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
112            0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
113            0.5 * (1.0 + xi_1 - xi_3) * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
114            0.5 * (1.0 + xi_2 - xi_3) * (1.0 - xi_2 - xi_3) * (1.0 - xi_1 - xi_3) / bottom,
115            xi_3 * (1.0 - xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
116            xi_3 * (1.0 + xi_1 - xi_3) * (1.0 - xi_2 - xi_3) / bottom,
117            xi_3 * (1.0 + xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
118            xi_3 * (1.0 - xi_1 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
119        ]
120        .into()
121    }
122    fn shape_functions_gradients(
123        parametric_coordinate: ParametricCoordinate<M>,
124    ) -> ShapeFunctionsGradients<M, N> {
125        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
126        let bottom = bottom(xi_3);
127        let bottom_squared = bottom * bottom;
128        [
129            [
130                0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_2 + xi_2 * xi_3 / bottom)
131                    - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
132                0.25 * ((-xi_1 - xi_2 - 1.0) * (-1.0 + xi_1 + xi_1 * xi_3 / bottom)
133                    - ((1.0 - xi_1) * (1.0 - xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
134                0.25 * (-xi_1 - xi_2 - 1.0)
135                    * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
136            ],
137            [
138                0.25 * ((xi_1 - xi_2 - 1.0) * (1.0 - xi_2 - xi_2 * xi_3 / bottom)
139                    + ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
140                0.25 * ((xi_1 - xi_2 - 1.0) * (-1.0 - xi_1 - xi_1 * xi_3 / bottom)
141                    - ((1.0 + xi_1) * (1.0 - xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
142                0.25 * (xi_1 - xi_2 - 1.0)
143                    * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
144            ],
145            [
146                0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_2 + xi_2 * xi_3 / bottom)
147                    + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
148                0.25 * ((xi_1 + xi_2 - 1.0) * (1.0 + xi_1 + xi_1 * xi_3 / bottom)
149                    + ((1.0 + xi_1) * (1.0 + xi_2) - xi_3 + xi_1 * xi_2 * xi_3 / bottom)),
150                0.25 * (xi_1 + xi_2 - 1.0)
151                    * (-1.0 + xi_1 * xi_2 / bottom + xi_1 * xi_2 * xi_3 / bottom_squared),
152            ],
153            [
154                0.25 * ((-xi_1 + xi_2 - 1.0) * (-1.0 - xi_2 - xi_2 * xi_3 / bottom)
155                    - ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
156                0.25 * ((-xi_1 + xi_2 - 1.0) * (1.0 - xi_1 - xi_1 * xi_3 / bottom)
157                    + ((1.0 - xi_1) * (1.0 + xi_2) - xi_3 - xi_1 * xi_2 * xi_3 / bottom)),
158                0.25 * (-xi_1 + xi_2 - 1.0)
159                    * (-1.0 - xi_1 * xi_2 / bottom - xi_1 * xi_2 * xi_3 / bottom_squared),
160            ],
161            [0.0, 0.0, 4.0 * xi_3 - 1.0],
162            [
163                -xi_1 * (1.0 - xi_2 - xi_3) / bottom,
164                -0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
165                0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
166            ],
167            [
168                0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
169                -xi_2 * (1.0 + xi_1 - xi_3) / bottom,
170                -0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
171            ],
172            [
173                -xi_1 * (1.0 + xi_2 - xi_3) / bottom,
174                0.5 * (1.0 - xi_1 - xi_3) * (1.0 + xi_1 - xi_3) / bottom,
175                -0.5 * (xi_1 * xi_1 * xi_2 / bottom_squared + xi_2) - 1.0 + xi_3,
176            ],
177            [
178                -0.5 * (1.0 - xi_2 - xi_3) * (1.0 + xi_2 - xi_3) / bottom,
179                -xi_2 * (1.0 - xi_1 - xi_3) / bottom,
180                0.5 * (xi_1 * xi_2 * xi_2 / bottom_squared + xi_1) - 1.0 + xi_3,
181            ],
182            [
183                -(1.0 - xi_2 - xi_3) * xi_3 / bottom,
184                -(1.0 - xi_1 - xi_3) * xi_3 / bottom,
185                xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 - xi_2 - 2.0 * xi_3,
186            ],
187            [
188                (1.0 - xi_2 - xi_3) * xi_3 / bottom,
189                -(1.0 + xi_1 - xi_3) * xi_3 / bottom,
190                -xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 - xi_2 - 2.0 * xi_3,
191            ],
192            [
193                (1.0 + xi_2 - xi_3) * xi_3 / bottom,
194                (1.0 + xi_1 - xi_3) * xi_3 / bottom,
195                xi_1 * xi_2 / bottom_squared + 1.0 + xi_1 + xi_2 - 2.0 * xi_3,
196            ],
197            [
198                -(1.0 + xi_2 - xi_3) * xi_3 / bottom,
199                (1.0 - xi_1 - xi_3) * xi_3 / bottom,
200                -xi_1 * xi_2 / bottom_squared + 1.0 - xi_1 + xi_2 - 2.0 * xi_3,
201            ],
202        ]
203        .into()
204    }
205}
206
207fn bottom(xi_3: Scalar) -> Scalar {
208    const SMALL: Scalar = 4e1 * f64::EPSILON;
209    if (1.0 - xi_3).abs() > SMALL {
210        1.0 - xi_3
211    } else {
212        SMALL
213    }
214}
215
216impl QuadraticFiniteElement<G, N> for Pyramid {}
conspire/domain/fem/block/element/quadratic/pyramid/mod.rs

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