conspire/domain/fem/block/element/linear/tetrahedron/
mod.rs1#[cfg(test)]
2mod test;
3
4use crate::{
5 fem::block::element::{
6 ElementNodalEitherCoordinates, FiniteElement, ParametricCoordinate, ParametricCoordinates,
7 ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8 linear::{LinearElement, LinearFiniteElement, M},
9 },
10 math::{ScalarList, Tensor},
11};
12use std::f64::consts::SQRT_2;
13
14const G: usize = 1;
15const N: usize = 4;
16const P: usize = N;
17
18const EDGES: usize = 6;
19
20pub type Tetrahedron = LinearElement<G, N>;
21
22impl FiniteElement<G, M, N, P> for Tetrahedron {
23 fn integration_points() -> ParametricCoordinates<G, M> {
24 [[0.25; M]].into()
25 }
26 fn integration_weights(&self) -> &ScalarList<G> {
27 &self.integration_weights
28 }
29 fn parametric_reference() -> ParametricReference<M, N> {
30 [
31 [0.0, 0.0, 0.0],
32 [1.0, 0.0, 0.0],
33 [0.0, 1.0, 0.0],
34 [0.0, 0.0, 1.0],
35 ]
36 .into()
37 }
38 fn parametric_weights() -> ScalarList<G> {
39 [1.0 / 6.0; G].into()
40 }
41 fn scaled_jacobians<const I: usize>(
42 nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
43 ) -> ScalarList<P> {
44 let numerator = ((&nodal_coordinates[1] - &nodal_coordinates[0])
45 .cross(&(&nodal_coordinates[2] - &nodal_coordinates[0]))
46 * (&nodal_coordinates[3] - &nodal_coordinates[0]))
47 * SQRT_2;
48 let lengths = lengths(nodal_coordinates);
49 [
50 numerator / (lengths[0] * lengths[2] * lengths[3]),
51 numerator / (lengths[0] * lengths[1] * lengths[4]),
52 numerator / (lengths[1] * lengths[2] * lengths[5]),
53 numerator / (lengths[3] * lengths[4] * lengths[5]),
54 ]
55 .into()
56 }
57 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
58 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
59 [1.0 - xi_1 - xi_2 - xi_3, xi_1, xi_2, xi_3].into()
60 }
61 fn shape_functions_gradients(
62 _parametric_coordinate: ParametricCoordinate<M>,
63 ) -> ShapeFunctionsGradients<M, N> {
64 [
65 [-1.0, -1.0, -1.0],
66 [1.0, 0.0, 0.0],
67 [0.0, 1.0, 0.0],
68 [0.0, 0.0, 1.0],
69 ]
70 .into()
71 }
72}
73
74impl LinearFiniteElement<G, N> for Tetrahedron {}
75
76fn edges<const I: usize>(
77 nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
78) -> ElementNodalEitherCoordinates<I, 6> {
79 [
80 &nodal_coordinates[1] - &nodal_coordinates[0],
81 &nodal_coordinates[2] - &nodal_coordinates[1],
82 &nodal_coordinates[0] - &nodal_coordinates[2],
83 &nodal_coordinates[3] - &nodal_coordinates[0],
84 &nodal_coordinates[3] - &nodal_coordinates[1],
85 &nodal_coordinates[3] - &nodal_coordinates[2],
86 ]
87 .into()
88}
89
90fn lengths<const I: usize>(
91 nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
92) -> ScalarList<EDGES> {
93 edges(nodal_coordinates)
94 .into_iter()
95 .map(|edge| edge.norm())
96 .collect()
97}