1#[cfg(test)]
2mod test;
3
4use crate::{
5 fem::block::element::{
6 FiniteElement, ParametricCoordinate, ParametricCoordinates, ParametricReference,
7 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 shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
56 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
57 [
58 0.125 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0),
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.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2) * (1.0 - xi_3),
67 0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
68 0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 - xi_3),
69 0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
70 0.25 * (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3 * 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 * xi_1) * (1.0 - xi_2) * (1.0 + xi_3),
75 0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
76 0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2) * (1.0 + xi_3),
77 0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
78 ]
79 .into()
80 }
81 fn shape_functions_gradients(
82 parametric_coordinate: ParametricCoordinate<M>,
83 ) -> ShapeFunctionsGradients<M, N> {
84 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
85 [
86 [
87 0.125
88 * (-(1.0 - xi_2) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0)
89 - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
90 0.125
91 * (-(1.0 - xi_1) * (1.0 - xi_3) * (-xi_1 - xi_2 - xi_3 - 2.0)
92 - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
93 0.125
94 * (-(1.0 - xi_1) * (1.0 - xi_2) * (-xi_1 - xi_2 - xi_3 - 2.0)
95 - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
96 ],
97 [
98 0.125
99 * ((1.0 - xi_2) * (1.0 - xi_3) * (xi_1 - xi_2 - xi_3 - 2.0)
100 + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
101 0.125
102 * (-(1.0 + xi_1) * (1.0 - xi_3) * (xi_1 - xi_2 - xi_3 - 2.0)
103 - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
104 0.125
105 * (-(1.0 + xi_1) * (1.0 - xi_2) * (xi_1 - xi_2 - xi_3 - 2.0)
106 - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3)),
107 ],
108 [
109 0.125
110 * ((1.0 + xi_2) * (1.0 - xi_3) * (xi_1 + xi_2 - xi_3 - 2.0)
111 + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
112 0.125
113 * ((1.0 + xi_1) * (1.0 - xi_3) * (xi_1 + xi_2 - xi_3 - 2.0)
114 + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
115 0.125
116 * (-(1.0 + xi_1) * (1.0 + xi_2) * (xi_1 + xi_2 - xi_3 - 2.0)
117 - (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
118 ],
119 [
120 0.125
121 * (-(1.0 + xi_2) * (1.0 - xi_3) * (-xi_1 + xi_2 - xi_3 - 2.0)
122 - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
123 0.125
124 * ((1.0 - xi_1) * (1.0 - xi_3) * (-xi_1 + xi_2 - xi_3 - 2.0)
125 + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
126 0.125
127 * (-(1.0 - xi_1) * (1.0 + xi_2) * (-xi_1 + xi_2 - xi_3 - 2.0)
128 - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3)),
129 ],
130 [
131 0.125
132 * (-(1.0 - xi_2) * (1.0 + xi_3) * (-xi_1 - xi_2 + xi_3 - 2.0)
133 - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
134 0.125
135 * (-(1.0 - xi_1) * (1.0 + xi_3) * (-xi_1 - xi_2 + xi_3 - 2.0)
136 - (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
137 0.125
138 * ((1.0 - xi_1) * (1.0 - xi_2) * (-xi_1 - xi_2 + xi_3 - 2.0)
139 + (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
140 ],
141 [
142 0.125
143 * ((1.0 - xi_2) * (1.0 + xi_3) * (xi_1 - xi_2 + xi_3 - 2.0)
144 + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
145 0.125
146 * (-(1.0 + xi_1) * (1.0 + xi_3) * (xi_1 - xi_2 + xi_3 - 2.0)
147 - (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
148 0.125
149 * ((1.0 + xi_1) * (1.0 - xi_2) * (xi_1 - xi_2 + xi_3 - 2.0)
150 + (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3)),
151 ],
152 [
153 0.125
154 * ((1.0 + xi_2) * (1.0 + xi_3) * (xi_1 + xi_2 + xi_3 - 2.0)
155 + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
156 0.125
157 * ((1.0 + xi_1) * (1.0 + xi_3) * (xi_1 + xi_2 + xi_3 - 2.0)
158 + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
159 0.125
160 * ((1.0 + xi_1) * (1.0 + xi_2) * (xi_1 + xi_2 + xi_3 - 2.0)
161 + (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
162 ],
163 [
164 0.125
165 * (-(1.0 + xi_2) * (1.0 + xi_3) * (-xi_1 + xi_2 + xi_3 - 2.0)
166 - (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
167 0.125
168 * ((1.0 - xi_1) * (1.0 + xi_3) * (-xi_1 + xi_2 + xi_3 - 2.0)
169 + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
170 0.125
171 * ((1.0 - xi_1) * (1.0 + xi_2) * (-xi_1 + xi_2 + xi_3 - 2.0)
172 + (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3)),
173 ],
174 [
175 -0.5 * xi_1 * (1.0 - xi_2) * (1.0 - xi_3),
176 -0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_3),
177 -0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2),
178 ],
179 [
180 0.25 * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
181 -0.5 * xi_2 * (1.0 + xi_1) * (1.0 - xi_3),
182 -0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2),
183 ],
184 [
185 -0.5 * xi_1 * (1.0 + xi_2) * (1.0 - xi_3),
186 0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_3),
187 -0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2),
188 ],
189 [
190 -0.25 * (1.0 - xi_2 * xi_2) * (1.0 - xi_3),
191 -0.5 * xi_2 * (1.0 - xi_1) * (1.0 - xi_3),
192 -0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2),
193 ],
194 [
195 -0.25 * (1.0 - xi_2) * (1.0 - xi_3 * xi_3),
196 -0.25 * (1.0 - xi_1) * (1.0 - xi_3 * xi_3),
197 -0.5 * xi_3 * (1.0 - xi_1) * (1.0 - 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.5 * xi_1 * (1.0 - xi_2) * (1.0 + xi_3),
216 -0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_3),
217 0.25 * (1.0 - xi_1 * xi_1) * (1.0 - xi_2),
218 ],
219 [
220 0.25 * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
221 -0.5 * xi_2 * (1.0 + xi_1) * (1.0 + xi_3),
222 0.25 * (1.0 + xi_1) * (1.0 - xi_2 * xi_2),
223 ],
224 [
225 -0.5 * xi_1 * (1.0 + xi_2) * (1.0 + xi_3),
226 0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_3),
227 0.25 * (1.0 - xi_1 * xi_1) * (1.0 + xi_2),
228 ],
229 [
230 -0.25 * (1.0 - xi_2 * xi_2) * (1.0 + xi_3),
231 -0.5 * xi_2 * (1.0 - xi_1) * (1.0 + xi_3),
232 0.25 * (1.0 - xi_1) * (1.0 - xi_2 * xi_2),
233 ],
234 ]
235 .into()
236 }
237}
238
239impl QuadraticFiniteElement<G, N> for Hexahedron {}