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 = 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 scaled_jacobians<const I: usize>(
80 _nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
81 ) -> ScalarList<P> {
82 todo!()
83 }
84 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
85 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
86 [
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.125 * xi_1 * xi_2 * xi_3 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3),
92 -0.125 * xi_1 * xi_2 * xi_3 * (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3),
93 0.125 * xi_1 * xi_2 * xi_3 * (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
94 -0.125 * xi_1 * xi_2 * xi_3 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
95 0.25 * xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 - xi_3),
96 -0.25 * xi_1 * xi_3 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
97 -0.25 * xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 - xi_3),
98 0.25 * xi_1 * xi_3 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
99 0.25 * xi_1 * xi_2 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
100 -0.25 * xi_1 * xi_2 * (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
101 0.25 * xi_1 * xi_2 * (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
102 -0.25 * xi_1 * xi_2 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
103 -0.25 * xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 + xi_3),
104 0.25 * xi_1 * xi_3 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
105 0.25 * xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
106 -0.25 * xi_1 * xi_3 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
107 (1.0 - xi_1 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3 * xi_3),
108 -0.5 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
109 0.5 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
110 -0.5 * xi_1 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3 * xi_3),
111 0.5 * xi_1 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3 * xi_3),
112 -0.5 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
113 0.5 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
114 ]
115 .into()
116 }
117 fn shape_functions_gradients(
118 parametric_coordinate: ParametricCoordinate<M>,
119 ) -> ShapeFunctionsGradients<M, N> {
120 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
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.125 * xi_2 * xi_3 * (1.0 - 2.0 * xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
159 -0.125 * xi_1 * xi_3 * (1.0 - xi_1) * (1.0 + 2.0 * xi_2) * (1.0 + xi_3),
160 -0.125 * xi_1 * xi_2 * (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + 2.0 * xi_3),
161 ],
162 [
163 -0.5 * xi_1 * xi_2 * xi_3 * (1.0 - xi_2) * (1.0 - xi_3),
164 0.25 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - 2.0 * xi_2) * (1.0 - xi_3),
165 0.25 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 - 2.0 * xi_3),
166 ],
167 [
168 -0.25 * xi_3 * (1.0 + 2.0 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
169 0.5 * xi_1 * xi_2 * xi_3 * (1.0 + xi_1) * (1.0 - xi_3),
170 -0.25 * xi_1 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - 2.0 * xi_3),
171 ],
172 [
173 0.5 * xi_1 * xi_2 * xi_3 * (1.0 + xi_2) * (1.0 - xi_3),
174 -0.25 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + 2.0 * xi_2) * (1.0 - xi_3),
175 -0.25 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 - 2.0 * xi_3),
176 ],
177 [
178 0.25 * xi_3 * (1.0 - 2.0 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
179 -0.5 * xi_1 * xi_2 * xi_3 * (1.0 - xi_1) * (1.0 - xi_3),
180 0.25 * xi_1 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - 2.0 * xi_3),
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.25 * xi_2 * (1.0 - 2.0 * xi_1) * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
199 -0.25 * xi_1 * (1.0 - xi_1) * (1.0 + 2.0 * xi_2) * (1.0 - xi_3 * xi_3),
200 0.5 * xi_1 * xi_2 * xi_3 * (1.0 - xi_1) * (1.0 + xi_2),
201 ],
202 [
203 0.5 * xi_1 * xi_2 * xi_3 * (1.0 - xi_2) * (1.0 + xi_3),
204 -0.25 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - 2.0 * xi_2) * (1.0 + xi_3),
205 -0.25 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 + 2.0 * xi_3),
206 ],
207 [
208 0.25 * xi_3 * (1.0 + 2.0 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
209 -0.5 * xi_1 * xi_2 * xi_3 * (1.0 + xi_1) * (1.0 + xi_3),
210 0.25 * xi_1 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + 2.0 * xi_3),
211 ],
212 [
213 -0.5 * xi_1 * xi_2 * xi_3 * (1.0 + xi_2) * (1.0 + xi_3),
214 0.25 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + 2.0 * xi_2) * (1.0 + xi_3),
215 0.25 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 + 2.0 * xi_3),
216 ],
217 [
218 -0.25 * xi_3 * (1.0 - 2.0 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
219 0.5 * xi_1 * xi_2 * xi_3 * (1.0 - xi_1) * (1.0 + xi_3),
220 -0.25 * xi_1 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + 2.0 * xi_3),
221 ],
222 [
223 -2.0 * xi_1 * (1.0 - xi_2 * xi_2) * (1.0 - xi_3 * xi_3),
224 -2.0 * xi_2 * (1.0 - xi_1 * xi_1) * (1.0 - xi_3 * xi_3),
225 -2.0 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2 * xi_2),
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 -xi_1 * xi_3 * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
234 -xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + xi_3),
235 0.5 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + 2.0 * xi_3),
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 0.5 * (1.0 + 2.0 * xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3 * xi_3),
244 -xi_1 * xi_2 * (1.0 + xi_1) * (1.0 - xi_3 * xi_3),
245 -xi_1 * xi_3 * (1.0 + xi_1) * (1.0 - xi_2 * 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 -xi_1 * xi_2 * (1.0 + xi_2) * (1.0 - xi_3 * xi_3),
254 0.5 * (1.0 - xi_1 * xi_1) * (1.0 + 2.0 * xi_2) * (1.0 - xi_3 * xi_3),
255 -xi_2 * xi_3 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2),
256 ],
257 ]
258 .into()
259 }
260}
261
262impl QuadraticFiniteElement<G, N> for Hexahedron {}