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