1#[cfg(test)]
2mod test;
3
4use crate::{
5 fem::block::element::{
6 Element, ElementNodalReferenceCoordinates, FiniteElement, GradientVectors,
7 ShapeFunctionsAtIntegrationPoints, StandardGradientOperators,
8 StandardGradientOperatorsTransposed,
9 composite::{
10 NormalizedProjectionMatrix, ParametricGradientOperators, ProjectionMatrix,
11 ShapeFunctionIntegrals, ShapeFunctionIntegralsProducts,
12 },
13 },
14 math::{Scalar, Scalars, Tensor, TensorRank1},
15};
16
17const G: usize = 4;
18const M: usize = 3;
19const N: usize = 10;
20const P: usize = 12;
21const Q: usize = 4;
22
23pub type Tetrahedron = Element<G, N>;
24
25impl FiniteElement<G, N> for Tetrahedron {
26 fn initialize(
27 reference_nodal_coordinates: ElementNodalReferenceCoordinates<N>,
28 ) -> (GradientVectors<G, N>, Scalars<G>) {
29 let gradient_vectors = Self::projected_gradient_vectors(&reference_nodal_coordinates);
30 let integration_weights =
31 Self::reference_jacobians(&reference_nodal_coordinates) * Self::integration_weight();
32 (gradient_vectors, integration_weights)
33 }
34 fn reset(&mut self) {
35 let (gradient_vectors, integration_weights) = Self::initialize(Self::reference());
36 self.gradient_vectors = gradient_vectors;
37 self.integration_weights = integration_weights;
38 }
39}
40
41impl Tetrahedron {
42 const fn integration_weight() -> Scalar {
43 1.0 / 24.0
44 }
45 const fn inverse_normalized_projection_matrix() -> NormalizedProjectionMatrix<Q> {
46 const DIAG: Scalar = 4.0 / 640.0;
47 const OFF: Scalar = -1.0 / 640.0;
48 NormalizedProjectionMatrix::<Q>::const_from([
49 [DIAG, OFF, OFF, OFF],
50 [OFF, DIAG, OFF, OFF],
51 [OFF, OFF, DIAG, OFF],
52 [OFF, OFF, OFF, DIAG],
53 ])
54 }
55 fn inverse_projection_matrix(
56 reference_jacobians_subelements: &Scalars<P>,
57 ) -> NormalizedProjectionMatrix<Q> {
58 Self::shape_function_integrals_products()
59 .iter()
60 .zip(reference_jacobians_subelements.iter())
61 .map(
62 |(shape_function_integrals_products, reference_jacobian_subelement)| {
63 shape_function_integrals_products * reference_jacobian_subelement
64 },
65 )
66 .sum::<ProjectionMatrix<Q>>()
67 .inverse()
68 }
69 fn projected_gradient_vectors(
70 reference_nodal_coordinates: &ElementNodalReferenceCoordinates<N>,
71 ) -> GradientVectors<G, N> {
72 let parametric_gradient_operators = Self::standard_gradient_operators()
73 .iter()
74 .map(|standard_gradient_operator| {
75 reference_nodal_coordinates * standard_gradient_operator
76 })
77 .collect::<ParametricGradientOperators<P>>();
78 let reference_jacobians_subelements =
79 Self::reference_jacobians_subelements(reference_nodal_coordinates);
80 let inverse_projection_matrix =
81 Self::inverse_projection_matrix(&reference_jacobians_subelements);
82 Self::shape_functions_at_integration_points()
83 .iter()
84 .map(|shape_functions_at_integration_point| {
85 Self::standard_gradient_operators_transposed()
86 .iter()
87 .map(|standard_gradient_operators_a| {
88 Self::shape_function_integrals()
89 .iter()
90 .zip(
91 standard_gradient_operators_a.iter().zip(
92 parametric_gradient_operators
93 .iter()
94 .zip(reference_jacobians_subelements.iter()),
95 ),
96 )
97 .map(
98 |(
99 shape_function_integral,
100 (
101 standard_gradient_operator,
102 (
103 parametric_gradient_operator,
104 reference_jacobian_subelement,
105 ),
106 ),
107 )| {
108 (parametric_gradient_operator.inverse_transpose()
109 * standard_gradient_operator)
110 * reference_jacobian_subelement
111 * (shape_functions_at_integration_point
112 * (&inverse_projection_matrix
113 * shape_function_integral))
114 },
115 )
116 .sum()
117 })
118 .collect()
119 })
120 .collect()
121 }
122 const fn reference() -> ElementNodalReferenceCoordinates<N> {
123 ElementNodalReferenceCoordinates::<N>::const_from([
124 [0.0, 0.0, 0.0],
125 [1.0, 0.0, 0.0],
126 [0.0, 1.0, 0.0],
127 [0.0, 0.0, 1.0],
128 [0.25, 0.25, 0.0],
129 [0.25, 0.0, 0.25],
130 [0.0, 0.25, 0.25],
131 [0.5, 0.5, 0.0],
132 [0.5, 0.0, 0.5],
133 [0.0, 0.5, 0.5],
134 ])
135 }
136 fn reference_jacobians(
137 reference_nodal_coordinates: &ElementNodalReferenceCoordinates<N>,
138 ) -> Scalars<G> {
139 let vector = Self::inverse_normalized_projection_matrix()
140 * Self::shape_function_integrals()
141 .iter()
142 .zip(Self::reference_jacobians_subelements(reference_nodal_coordinates).iter())
143 .map(|(shape_function_integral, reference_jacobian_subelement)| {
144 shape_function_integral * reference_jacobian_subelement
145 })
146 .sum::<TensorRank1<Q, 9>>();
147 Self::shape_functions_at_integration_points()
148 .iter()
149 .map(|shape_functions_at_integration_point| {
150 shape_functions_at_integration_point * &vector
151 })
152 .collect()
153 }
154 fn reference_jacobians_subelements(
155 reference_nodal_coordinates: &ElementNodalReferenceCoordinates<N>,
156 ) -> Scalars<P> {
157 Self::standard_gradient_operators()
158 .iter()
159 .map(|standard_gradient_operator| {
160 reference_nodal_coordinates * standard_gradient_operator
161 })
162 .collect::<ParametricGradientOperators<P>>()
163 .iter()
164 .map(|parametric_gradient_operator| parametric_gradient_operator.determinant())
165 .collect()
166 }
167 const fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, Q> {
168 const DIAG: Scalar = 0.585_410_196_624_968_5;
169 const OFF: Scalar = 0.138_196_601_125_010_5;
170 ShapeFunctionsAtIntegrationPoints::<G, Q>::const_from([
171 [DIAG, OFF, OFF, OFF],
172 [OFF, DIAG, OFF, OFF],
173 [OFF, OFF, DIAG, OFF],
174 [OFF, OFF, OFF, DIAG],
175 ])
176 }
177 const fn shape_function_integrals() -> ShapeFunctionIntegrals<P, Q> {
178 ShapeFunctionIntegrals::<P, Q>::const_from([
179 [200.0, 40.0, 40.0, 40.0],
180 [40.0, 200.0, 40.0, 40.0],
181 [40.0, 40.0, 200.0, 40.0],
182 [40.0, 40.0, 40.0, 200.0],
183 [30.0, 70.0, 30.0, 30.0],
184 [10.0, 50.0, 50.0, 50.0],
185 [30.0, 30.0, 30.0, 70.0],
186 [50.0, 50.0, 10.0, 50.0],
187 [50.0, 50.0, 50.0, 10.0],
188 [30.0, 30.0, 70.0, 30.0],
189 [50.0, 10.0, 50.0, 50.0],
190 [70.0, 30.0, 30.0, 30.0],
191 ])
192 }
193 const fn shape_function_integrals_products() -> ShapeFunctionIntegralsProducts<P, Q> {
194 ShapeFunctionIntegralsProducts::<P, Q>::const_from([
195 [
196 [128.0, 24.0, 24.0, 24.0],
197 [24.0, 8.0, 4.0, 4.0],
198 [24.0, 4.0, 8.0, 4.0],
199 [24.0, 4.0, 4.0, 8.0],
200 ],
201 [
202 [8.0, 24.0, 4.0, 4.0],
203 [24.0, 128.0, 24.0, 24.0],
204 [4.0, 24.0, 8.0, 4.0],
205 [4.0, 24.0, 4.0, 8.0],
206 ],
207 [
208 [8.0, 4.0, 24.0, 4.0],
209 [4.0, 8.0, 24.0, 4.0],
210 [24.0, 24.0, 128.0, 24.0],
211 [4.0, 4.0, 24.0, 8.0],
212 ],
213 [
214 [8.0, 4.0, 4.0, 24.0],
215 [4.0, 8.0, 4.0, 24.0],
216 [4.0, 4.0, 8.0, 24.0],
217 [24.0, 24.0, 24.0, 128.0],
218 ],
219 [
220 [7.0, 13.0, 5.0, 5.0],
221 [13.0, 31.0, 13.0, 13.0],
222 [5.0, 13.0, 7.0, 5.0],
223 [5.0, 13.0, 5.0, 7.0],
224 ],
225 [
226 [1.0, 3.0, 3.0, 3.0],
227 [3.0, 17.0, 15.0, 15.0],
228 [3.0, 15.0, 17.0, 15.0],
229 [3.0, 15.0, 15.0, 17.0],
230 ],
231 [
232 [7.0, 5.0, 5.0, 13.0],
233 [5.0, 7.0, 5.0, 13.0],
234 [5.0, 5.0, 7.0, 13.0],
235 [13.0, 13.0, 13.0, 31.0],
236 ],
237 [
238 [17.0, 15.0, 3.0, 15.0],
239 [15.0, 17.0, 3.0, 15.0],
240 [3.0, 3.0, 1.0, 3.0],
241 [15.0, 15.0, 3.0, 17.0],
242 ],
243 [
244 [17.0, 15.0, 15.0, 3.0],
245 [15.0, 17.0, 15.0, 3.0],
246 [15.0, 15.0, 17.0, 3.0],
247 [3.0, 3.0, 3.0, 1.0],
248 ],
249 [
250 [7.0, 5.0, 13.0, 5.0],
251 [5.0, 7.0, 13.0, 5.0],
252 [13.0, 13.0, 31.0, 13.0],
253 [5.0, 5.0, 13.0, 7.0],
254 ],
255 [
256 [17.0, 3.0, 15.0, 15.0],
257 [3.0, 1.0, 3.0, 3.0],
258 [15.0, 3.0, 17.0, 15.0],
259 [15.0, 3.0, 15.0, 17.0],
260 ],
261 [
262 [31.0, 13.0, 13.0, 13.0],
263 [13.0, 7.0, 5.0, 5.0],
264 [13.0, 5.0, 7.0, 5.0],
265 [13.0, 5.0, 5.0, 7.0],
266 ],
267 ])
268 }
269 const fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
270 StandardGradientOperators::<M, N, P>::const_from([
271 [
272 [-2.0, -2.0, -2.0],
273 [0.0, 0.0, 0.0],
274 [0.0, 0.0, 0.0],
275 [0.0, 0.0, 0.0],
276 [2.0, 0.0, 0.0],
277 [0.0, 0.0, 0.0],
278 [0.0, 2.0, 0.0],
279 [0.0, 0.0, 2.0],
280 [0.0, 0.0, 0.0],
281 [0.0, 0.0, 0.0],
282 ],
283 [
284 [0.0, 0.0, 0.0],
285 [2.0, 0.0, 0.0],
286 [0.0, 0.0, 0.0],
287 [0.0, 0.0, 0.0],
288 [-2.0, -2.0, -2.0],
289 [0.0, 2.0, 0.0],
290 [0.0, 0.0, 0.0],
291 [0.0, 0.0, 0.0],
292 [0.0, 0.0, 2.0],
293 [0.0, 0.0, 0.0],
294 ],
295 [
296 [0.0, 0.0, 0.0],
297 [0.0, 0.0, 0.0],
298 [0.0, 2.0, 0.0],
299 [0.0, 0.0, 0.0],
300 [0.0, 0.0, 0.0],
301 [2.0, 0.0, 0.0],
302 [-2.0, -2.0, -2.0],
303 [0.0, 0.0, 0.0],
304 [0.0, 0.0, 0.0],
305 [0.0, 0.0, 2.0],
306 ],
307 [
308 [0.0, 0.0, 0.0],
309 [0.0, 0.0, 0.0],
310 [0.0, 0.0, 0.0],
311 [0.0, 0.0, 2.0],
312 [0.0, 0.0, 0.0],
313 [0.0, 0.0, 0.0],
314 [0.0, 0.0, 0.0],
315 [-2.0, -2.0, -2.0],
316 [2.0, 0.0, 0.0],
317 [0.0, 2.0, 0.0],
318 ],
319 [
320 [0.0, 0.0, 0.0],
321 [0.0, 0.0, 0.0],
322 [0.0, 0.0, 0.0],
323 [0.0, 0.0, 0.0],
324 [-2.0 / 3.0, -2.0, -2.0],
325 [4.0 / 3.0, 2.0, 0.0],
326 [-2.0 / 3.0, 0.0, 0.0],
327 [-2.0 / 3.0, 0.0, 0.0],
328 [4.0 / 3.0, 0.0, 2.0],
329 [-2.0 / 3.0, 0.0, 0.0],
330 ],
331 [
332 [0.0, 0.0, 0.0],
333 [0.0, 0.0, 0.0],
334 [0.0, 0.0, 0.0],
335 [0.0, 0.0, 0.0],
336 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
337 [4.0 / 3.0, 4.0 / 3.0, -2.0 / 3.0],
338 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
339 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
340 [4.0 / 3.0, -2.0 / 3.0, 4.0 / 3.0],
341 [-2.0 / 3.0, 4.0 / 3.0, 4.0 / 3.0],
342 ],
343 [
344 [0.0, 0.0, 0.0],
345 [0.0, 0.0, 0.0],
346 [0.0, 0.0, 0.0],
347 [0.0, 0.0, 0.0],
348 [0.0, 0.0, -2.0 / 3.0],
349 [0.0, 0.0, -2.0 / 3.0],
350 [0.0, 0.0, -2.0 / 3.0],
351 [-2.0, -2.0, -2.0 / 3.0],
352 [2.0, 0.0, 4.0 / 3.0],
353 [0.0, 2.0, 4.0 / 3.0],
354 ],
355 [
356 [0.0, 0.0, 0.0],
357 [0.0, 0.0, 0.0],
358 [0.0, 0.0, 0.0],
359 [0.0, 0.0, 0.0],
360 [0.0, -4.0 / 3.0, -2.0],
361 [0.0, 2.0 / 3.0, 0.0],
362 [0.0, 2.0 / 3.0, 0.0],
363 [-2.0, -4.0 / 3.0, 0.0],
364 [2.0, 2.0 / 3.0, 2.0],
365 [0.0, 2.0 / 3.0, 0.0],
366 ],
367 [
368 [0.0, 0.0, 0.0],
369 [0.0, 0.0, 0.0],
370 [0.0, 0.0, 0.0],
371 [0.0, 0.0, 0.0],
372 [0.0, -2.0, -4.0 / 3.0],
373 [2.0, 2.0, 2.0 / 3.0],
374 [-2.0, 0.0, -4.0 / 3.0],
375 [0.0, 0.0, 2.0 / 3.0],
376 [0.0, 0.0, 2.0 / 3.0],
377 [0.0, 0.0, 2.0 / 3.0],
378 ],
379 [
380 [0.0, 0.0, 0.0],
381 [0.0, 0.0, 0.0],
382 [0.0, 0.0, 0.0],
383 [0.0, 0.0, 0.0],
384 [0.0, -2.0 / 3.0, 0.0],
385 [2.0, 4.0 / 3.0, 0.0],
386 [-2.0, -2.0 / 3.0, -2.0],
387 [0.0, -2.0 / 3.0, 0.0],
388 [0.0, -2.0 / 3.0, 0.0],
389 [0.0, 4.0 / 3.0, 2.0],
390 ],
391 [
392 [0.0, 0.0, 0.0],
393 [0.0, 0.0, 0.0],
394 [0.0, 0.0, 0.0],
395 [0.0, 0.0, 0.0],
396 [2.0 / 3.0, 0.0, 0.0],
397 [2.0 / 3.0, 0.0, 0.0],
398 [-4.0 / 3.0, 0.0, -2.0],
399 [-4.0 / 3.0, -2.0, 0.0],
400 [2.0 / 3.0, 0.0, 0.0],
401 [2.0 / 3.0, 2.0, 2.0],
402 ],
403 [
404 [0.0, 0.0, 0.0],
405 [0.0, 0.0, 0.0],
406 [0.0, 0.0, 0.0],
407 [0.0, 0.0, 0.0],
408 [2.0 / 3.0, -4.0 / 3.0, -4.0 / 3.0],
409 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
410 [-4.0 / 3.0, 2.0 / 3.0, -4.0 / 3.0],
411 [-4.0 / 3.0, -4.0 / 3.0, 2.0 / 3.0],
412 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
413 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
414 ],
415 ])
416 }
417 const fn standard_gradient_operators_transposed() -> StandardGradientOperatorsTransposed<M, N, P>
418 {
419 StandardGradientOperatorsTransposed::<M, N, P>::const_from([
420 [
421 [-2.0, -2.0, -2.0],
422 [0.0, 0.0, 0.0],
423 [0.0, 0.0, 0.0],
424 [0.0, 0.0, 0.0],
425 [0.0, 0.0, 0.0],
426 [0.0, 0.0, 0.0],
427 [0.0, 0.0, 0.0],
428 [0.0, 0.0, 0.0],
429 [0.0, 0.0, 0.0],
430 [0.0, 0.0, 0.0],
431 [0.0, 0.0, 0.0],
432 [0.0, 0.0, 0.0],
433 ],
434 [
435 [0.0, 0.0, 0.0],
436 [2.0, 0.0, 0.0],
437 [0.0, 0.0, 0.0],
438 [0.0, 0.0, 0.0],
439 [0.0, 0.0, 0.0],
440 [0.0, 0.0, 0.0],
441 [0.0, 0.0, 0.0],
442 [0.0, 0.0, 0.0],
443 [0.0, 0.0, 0.0],
444 [0.0, 0.0, 0.0],
445 [0.0, 0.0, 0.0],
446 [0.0, 0.0, 0.0],
447 ],
448 [
449 [0.0, 0.0, 0.0],
450 [0.0, 0.0, 0.0],
451 [0.0, 2.0, 0.0],
452 [0.0, 0.0, 0.0],
453 [0.0, 0.0, 0.0],
454 [0.0, 0.0, 0.0],
455 [0.0, 0.0, 0.0],
456 [0.0, 0.0, 0.0],
457 [0.0, 0.0, 0.0],
458 [0.0, 0.0, 0.0],
459 [0.0, 0.0, 0.0],
460 [0.0, 0.0, 0.0],
461 ],
462 [
463 [0.0, 0.0, 0.0],
464 [0.0, 0.0, 0.0],
465 [0.0, 0.0, 0.0],
466 [0.0, 0.0, 2.0],
467 [0.0, 0.0, 0.0],
468 [0.0, 0.0, 0.0],
469 [0.0, 0.0, 0.0],
470 [0.0, 0.0, 0.0],
471 [0.0, 0.0, 0.0],
472 [0.0, 0.0, 0.0],
473 [0.0, 0.0, 0.0],
474 [0.0, 0.0, 0.0],
475 ],
476 [
477 [2.0, 0.0, 0.0],
478 [-2.0, -2.0, -2.0],
479 [0.0, 0.0, 0.0],
480 [0.0, 0.0, 0.0],
481 [-2.0 / 3.0, -2.0, -2.0],
482 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
483 [0.0, 0.0, -2.0 / 3.0],
484 [0.0, -4.0 / 3.0, -2.0],
485 [0.0, -2.0, -4.0 / 3.0],
486 [0.0, -2.0 / 3.0, 0.0],
487 [2.0 / 3.0, 0.0, 0.0],
488 [2.0 / 3.0, -4.0 / 3.0, -4.0 / 3.0],
489 ],
490 [
491 [0.0, 0.0, 0.0],
492 [0.0, 2.0, 0.0],
493 [2.0, 0.0, 0.0],
494 [0.0, 0.0, 0.0],
495 [4.0 / 3.0, 2.0, 0.0],
496 [4.0 / 3.0, 4.0 / 3.0, -2.0 / 3.0],
497 [0.0, 0.0, -2.0 / 3.0],
498 [0.0, 2.0 / 3.0, 0.0],
499 [2.0, 2.0, 2.0 / 3.0],
500 [2.0, 4.0 / 3.0, 0.0],
501 [2.0 / 3.0, 0.0, 0.0],
502 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
503 ],
504 [
505 [0.0, 2.0, 0.0],
506 [0.0, 0.0, 0.0],
507 [-2.0, -2.0, -2.0],
508 [0.0, 0.0, 0.0],
509 [-2.0 / 3.0, 0.0, 0.0],
510 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
511 [0.0, 0.0, -2.0 / 3.0],
512 [0.0, 2.0 / 3.0, 0.0],
513 [-2.0, 0.0, -4.0 / 3.0],
514 [-2.0, -2.0 / 3.0, -2.0],
515 [-4.0 / 3.0, 0.0, -2.0],
516 [-4.0 / 3.0, 2.0 / 3.0, -4.0 / 3.0],
517 ],
518 [
519 [0.0, 0.0, 2.0],
520 [0.0, 0.0, 0.0],
521 [0.0, 0.0, 0.0],
522 [-2.0, -2.0, -2.0],
523 [-2.0 / 3.0, 0.0, 0.0],
524 [-2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0],
525 [-2.0, -2.0, -2.0 / 3.0],
526 [-2.0, -4.0 / 3.0, 0.0],
527 [0.0, 0.0, 2.0 / 3.0],
528 [0.0, -2.0 / 3.0, 0.0],
529 [-4.0 / 3.0, -2.0, 0.0],
530 [-4.0 / 3.0, -4.0 / 3.0, 2.0 / 3.0],
531 ],
532 [
533 [0.0, 0.0, 0.0],
534 [0.0, 0.0, 2.0],
535 [0.0, 0.0, 0.0],
536 [2.0, 0.0, 0.0],
537 [4.0 / 3.0, 0.0, 2.0],
538 [4.0 / 3.0, -2.0 / 3.0, 4.0 / 3.0],
539 [2.0, 0.0, 4.0 / 3.0],
540 [2.0, 2.0 / 3.0, 2.0],
541 [0.0, 0.0, 2.0 / 3.0],
542 [0.0, -2.0 / 3.0, 0.0],
543 [2.0 / 3.0, 0.0, 0.0],
544 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
545 ],
546 [
547 [0.0, 0.0, 0.0],
548 [0.0, 0.0, 0.0],
549 [0.0, 0.0, 2.0],
550 [0.0, 2.0, 0.0],
551 [-2.0 / 3.0, 0.0, 0.0],
552 [-2.0 / 3.0, 4.0 / 3.0, 4.0 / 3.0],
553 [0.0, 2.0, 4.0 / 3.0],
554 [0.0, 2.0 / 3.0, 0.0],
555 [0.0, 0.0, 2.0 / 3.0],
556 [0.0, 4.0 / 3.0, 2.0],
557 [2.0 / 3.0, 2.0, 2.0],
558 [2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0],
559 ],
560 ])
561 }
562}