Skip to main content

conspire/domain/fem/block/element/serendipity/hexahedron/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    fem::block::element::{
6        FiniteElement, IntegrationWeights, ParametricCoordinate, ParametricCoordinates,
7        ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8        quadratic::{Hexahedron as QuadraticHexahedron, QuadraticElement, QuadraticFiniteElement},
9        serendipity::M,
10    },
11    math::ScalarList,
12    units::Volume,
13};
14
15const G: usize = 27;
16const N: usize = 20;
17const P: usize = N;
18
19pub type Hexahedron = QuadraticElement<G, N>;
20
21impl FiniteElement<G, M, N, P> for Hexahedron {
22    fn integration_points() -> ParametricCoordinates<G, M> {
23        QuadraticHexahedron::integration_points()
24    }
25    fn integration_weights(&self) -> &IntegrationWeights<G, Volume> {
26        &self.integration_weights
27    }
28    fn parametric_reference() -> ParametricReference<M, N> {
29        [
30            [-1.0, -1.0, -1.0],
31            [1.0, -1.0, -1.0],
32            [1.0, 1.0, -1.0],
33            [-1.0, 1.0, -1.0],
34            [-1.0, -1.0, 1.0],
35            [1.0, -1.0, 1.0],
36            [1.0, 1.0, 1.0],
37            [-1.0, 1.0, 1.0],
38            [0.0, -1.0, -1.0],
39            [1.0, 0.0, -1.0],
40            [0.0, 1.0, -1.0],
41            [-1.0, 0.0, -1.0],
42            [-1.0, -1.0, 0.0],
43            [1.0, -1.0, 0.0],
44            [1.0, 1.0, 0.0],
45            [-1.0, 1.0, 0.0],
46            [0.0, -1.0, 1.0],
47            [1.0, 0.0, 1.0],
48            [0.0, 1.0, 1.0],
49            [-1.0, 0.0, 1.0],
50        ]
51        .into()
52    }
53    fn parametric_weights() -> ScalarList<G> {
54        QuadraticHexahedron::parametric_weights()
55    }
56    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
57        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
58        [
59            0.125 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0),
60            0.125 * (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3) * (xi_1 - xi_2 - xi_3 - 2.0),
61            0.125 * (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3) * (xi_1 + xi_2 - xi_3 - 2.0),
62            0.125 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3) * (-xi_1 + xi_2 - xi_3 - 2.0),
63            0.125 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3) * (-xi_1 - xi_2 + xi_3 - 2.0),
64            0.125 * (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3) * (xi_1 - xi_2 + xi_3 - 2.0),
65            0.125 * (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3) * (xi_1 + xi_2 + xi_3 - 2.0),
66            0.125 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3) * (-xi_1 + xi_2 + xi_3 - 2.0),
67            0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 - xi_3),
68            0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
69            0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 - xi_3),
70            0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
71            0.25 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
72            0.25 * (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
73            0.25 * (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
74            0.25 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
75            0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 + xi_3),
76            0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
77            0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
78            0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
79        ]
80        .into()
81    }
82    fn shape_functions_gradients(
83        parametric_coordinate: ParametricCoordinate<M>,
84    ) -> ShapeFunctionsGradients<M, N> {
85        let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
86        [
87            [
88                0.125
89                    * (-(1.0 - xi_2) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0)
90                        - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
91                0.125
92                    * (-(1.0 - xi_1) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0)
93                        - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
94                0.125
95                    * (-(1.0 - xi_1) * (1.0 - xi_2) * (-xi_1 - xi_2 - xi_3 - 2.0)
96                        - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
97            ],
98            [
99                0.125
100                    * ((1.0 - xi_2) * (1.0 - xi_3) * (xi_1 - xi_2 - xi_3 - 2.0)
101                        + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
102                0.125
103                    * (-(1.0 + xi_1) * (1.0 - xi_3) * (xi_1 - xi_2 - xi_3 - 2.0)
104                        - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
105                0.125
106                    * (-(1.0 + xi_1) * (1.0 - xi_2) * (xi_1 - xi_2 - xi_3 - 2.0)
107                        - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
108            ],
109            [
110                0.125
111                    * ((1.0 + xi_2) * (1.0 - xi_3) * (xi_1 + xi_2 - xi_3 - 2.0)
112                        + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
113                0.125
114                    * ((1.0 + xi_1) * (1.0 - xi_3) * (xi_1 + xi_2 - xi_3 - 2.0)
115                        + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
116                0.125
117                    * (-(1.0 + xi_1) * (1.0 + xi_2) * (xi_1 + xi_2 - xi_3 - 2.0)
118                        - (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
119            ],
120            [
121                0.125
122                    * (-(1.0 + xi_2) * (1.0 - xi_3) * (-xi_1 + xi_2 - xi_3 - 2.0)
123                        - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
124                0.125
125                    * ((1.0 - xi_1) * (1.0 - xi_3) * (-xi_1 + xi_2 - xi_3 - 2.0)
126                        + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
127                0.125
128                    * (-(1.0 - xi_1) * (1.0 + xi_2) * (-xi_1 + xi_2 - xi_3 - 2.0)
129                        - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
130            ],
131            [
132                0.125
133                    * (-(1.0 - xi_2) * (1.0 + xi_3) * (-xi_1 - xi_2 + xi_3 - 2.0)
134                        - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
135                0.125
136                    * (-(1.0 - xi_1) * (1.0 + xi_3) * (-xi_1 - xi_2 + xi_3 - 2.0)
137                        - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
138                0.125
139                    * ((1.0 - xi_1) * (1.0 - xi_2) * (-xi_1 - xi_2 + xi_3 - 2.0)
140                        + (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
141            ],
142            [
143                0.125
144                    * ((1.0 - xi_2) * (1.0 + xi_3) * (xi_1 - xi_2 + xi_3 - 2.0)
145                        + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
146                0.125
147                    * (-(1.0 + xi_1) * (1.0 + xi_3) * (xi_1 - xi_2 + xi_3 - 2.0)
148                        - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
149                0.125
150                    * ((1.0 + xi_1) * (1.0 - xi_2) * (xi_1 - xi_2 + xi_3 - 2.0)
151                        + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
152            ],
153            [
154                0.125
155                    * ((1.0 + xi_2) * (1.0 + xi_3) * (xi_1 + xi_2 + xi_3 - 2.0)
156                        + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
157                0.125
158                    * ((1.0 + xi_1) * (1.0 + xi_3) * (xi_1 + xi_2 + xi_3 - 2.0)
159                        + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
160                0.125
161                    * ((1.0 + xi_1) * (1.0 + xi_2) * (xi_1 + xi_2 + xi_3 - 2.0)
162                        + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
163            ],
164            [
165                0.125
166                    * (-(1.0 + xi_2) * (1.0 + xi_3) * (-xi_1 + xi_2 + xi_3 - 2.0)
167                        - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
168                0.125
169                    * ((1.0 - xi_1) * (1.0 + xi_3) * (-xi_1 + xi_2 + xi_3 - 2.0)
170                        + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
171                0.125
172                    * ((1.0 - xi_1) * (1.0 + xi_2) * (-xi_1 + xi_2 + xi_3 - 2.0)
173                        + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
174            ],
175            [
176                -0.5 * xi_1 * (1.0 - xi_2) * (1.0 - xi_3),
177                -0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_3),
178                -0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2),
179            ],
180            [
181                0.25 * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
182                -0.5 * xi_2 * (1.0 + xi_1) * (1.0 - xi_3),
183                -0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2),
184            ],
185            [
186                -0.5 * xi_1 * (1.0 + xi_2) * (1.0 - xi_3),
187                0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_3),
188                -0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2),
189            ],
190            [
191                -0.25 * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
192                -0.5 * xi_2 * (1.0 - xi_1) * (1.0 - xi_3),
193                -0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2),
194            ],
195            [
196                -0.25 * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
197                -0.25 * (1.0 - xi_1) * (1.0 - xi_3 * xi_3),
198                -0.5 * xi_3 * (1.0 - xi_1) * (1.0 - xi_2),
199            ],
200            [
201                0.25 * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
202                -0.25 * (1.0 + xi_1) * (1.0 - xi_3 * xi_3),
203                -0.5 * xi_3 * (1.0 + xi_1) * (1.0 - xi_2),
204            ],
205            [
206                0.25 * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
207                0.25 * (1.0 + xi_1) * (1.0 - xi_3 * xi_3),
208                -0.5 * xi_3 * (1.0 + xi_1) * (1.0 + xi_2),
209            ],
210            [
211                -0.25 * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
212                0.25 * (1.0 - xi_1) * (1.0 - xi_3 * xi_3),
213                -0.5 * xi_3 * (1.0 - xi_1) * (1.0 + xi_2),
214            ],
215            [
216                -0.5 * xi_1 * (1.0 - xi_2) * (1.0 + xi_3),
217                -0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_3),
218                0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2),
219            ],
220            [
221                0.25 * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
222                -0.5 * xi_2 * (1.0 + xi_1) * (1.0 + xi_3),
223                0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2),
224            ],
225            [
226                -0.5 * xi_1 * (1.0 + xi_2) * (1.0 + xi_3),
227                0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_3),
228                0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2),
229            ],
230            [
231                -0.25 * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
232                -0.5 * xi_2 * (1.0 - xi_1) * (1.0 + xi_3),
233                0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2),
234            ],
235        ]
236        .into()
237    }
238}
239
240impl QuadraticFiniteElement<G, N> for Hexahedron {}