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