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