conspire/fem/block/element/linear/hexahedron/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::*;
5use crate::{math::tensor_rank_1, mechanics::Scalar};
6
7const G: usize = 8;
8const M: usize = 3;
9const N: usize = 8;
10const P: usize = N;
11
12#[cfg(test)]
13const Q: usize = N;
14
15const SQRT_3: Scalar = 1.732_050_807_568_877_2;
16
17pub type Hexahedron<'a, C> = Element<'a, C, G, N>;
18
19impl<'a, C> FiniteElement<'a, C, G, N> for Hexahedron<'a, C> {
20    fn new(
21        constitutive_model: &'a C,
22        reference_nodal_coordinates: ReferenceNodalCoordinates<N>,
23    ) -> Self {
24        let (gradient_vectors, integration_weights) = Self::initialize(reference_nodal_coordinates);
25        Self {
26            constitutive_model,
27            gradient_vectors,
28            integration_weights,
29        }
30    }
31    fn reference() -> ReferenceNodalCoordinates<N> {
32        ReferenceNodalCoordinates::const_from([
33            tensor_rank_1([-1.0, -1.0, -1.0]),
34            tensor_rank_1([1.0, -1.0, -1.0]),
35            tensor_rank_1([1.0, 1.0, -1.0]),
36            tensor_rank_1([-1.0, 1.0, -1.0]),
37            tensor_rank_1([-1.0, -1.0, 1.0]),
38            tensor_rank_1([1.0, -1.0, 1.0]),
39            tensor_rank_1([1.0, 1.0, 1.0]),
40            tensor_rank_1([-1.0, 1.0, 1.0]),
41        ])
42    }
43    fn reset(&mut self) {
44        let (gradient_vectors, integration_weights) = Self::initialize(Self::reference());
45        self.gradient_vectors = gradient_vectors;
46        self.integration_weights = integration_weights;
47    }
48}
49
50impl<'a, C> Hexahedron<'a, C> {
51    fn initialize(
52        reference_nodal_coordinates: ReferenceNodalCoordinates<N>,
53    ) -> (GradientVectors<G, N>, Scalars<G>) {
54        let standard_gradient_operators = Self::standard_gradient_operators();
55        let gradient_vectors = standard_gradient_operators
56            .iter()
57            .map(|standard_gradient_operator| {
58                (&reference_nodal_coordinates * standard_gradient_operator).inverse_transpose()
59                    * standard_gradient_operator
60            })
61            .collect();
62        let integration_weights = standard_gradient_operators
63            .iter()
64            .map(|standard_gradient_operator| {
65                (&reference_nodal_coordinates * standard_gradient_operator).determinant()
66                    * Self::integration_weight()
67            })
68            .collect();
69        (gradient_vectors, integration_weights)
70    }
71    const fn integration_point(point: usize) -> [Scalar; M] {
72        match point {
73            0 => [-SQRT_3 / 3.0, -SQRT_3 / 3.0, -SQRT_3 / 3.0],
74            1 => [-SQRT_3 / 3.0, -SQRT_3 / 3.0, SQRT_3 / 3.0],
75            2 => [-SQRT_3 / 3.0, SQRT_3 / 3.0, -SQRT_3 / 3.0],
76            3 => [-SQRT_3 / 3.0, SQRT_3 / 3.0, SQRT_3 / 3.0],
77            4 => [SQRT_3 / 3.0, -SQRT_3 / 3.0, -SQRT_3 / 3.0],
78            5 => [SQRT_3 / 3.0, -SQRT_3 / 3.0, SQRT_3 / 3.0],
79            6 => [SQRT_3 / 3.0, SQRT_3 / 3.0, -SQRT_3 / 3.0],
80            7 => [SQRT_3 / 3.0, SQRT_3 / 3.0, SQRT_3 / 3.0],
81            _ => panic!(),
82        }
83    }
84    const fn integration_weight() -> Scalar {
85        1.0
86    }
87    #[cfg(test)]
88    const fn shape_functions([xi_1, xi_2, xi_3]: [Scalar; M]) -> ShapeFunctions<N> {
89        tensor_rank_1([
90            (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
91            (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
92            (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
93            (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
94            (1.0 - xi_1) * (1.0 - xi_2) * (1.0 + xi_3) / 8.0,
95            (1.0 + xi_1) * (1.0 - xi_2) * (1.0 + xi_3) / 8.0,
96            (1.0 + xi_1) * (1.0 + xi_2) * (1.0 + xi_3) / 8.0,
97            (1.0 - xi_1) * (1.0 + xi_2) * (1.0 + xi_3) / 8.0,
98        ])
99    }
100    #[cfg(test)]
101    const fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, Q> {
102        ShapeFunctionsAtIntegrationPoints::const_from([
103            Self::shape_functions(Self::integration_point(0)),
104            Self::shape_functions(Self::integration_point(1)),
105            Self::shape_functions(Self::integration_point(2)),
106            Self::shape_functions(Self::integration_point(3)),
107            Self::shape_functions(Self::integration_point(4)),
108            Self::shape_functions(Self::integration_point(5)),
109            Self::shape_functions(Self::integration_point(6)),
110            Self::shape_functions(Self::integration_point(7)),
111        ])
112    }
113    const fn shape_functions_gradients(
114        [xi_1, xi_2, xi_3]: [Scalar; M],
115    ) -> ShapeFunctionsGradients<M, N> {
116        ShapeFunctionsGradients::const_from([
117            tensor_rank_1([
118                -(1.0 - xi_2) * (1.0 - xi_3) / 8.0,
119                -(1.0 - xi_1) * (1.0 - xi_3) / 8.0,
120                -(1.0 - xi_1) * (1.0 - xi_2) / 8.0,
121            ]),
122            tensor_rank_1([
123                (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
124                -(1.0 + xi_1) * (1.0 - xi_3) / 8.0,
125                -(1.0 + xi_1) * (1.0 - xi_2) / 8.0,
126            ]),
127            tensor_rank_1([
128                (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
129                (1.0 + xi_1) * (1.0 - xi_3) / 8.0,
130                -(1.0 + xi_1) * (1.0 + xi_2) / 8.0,
131            ]),
132            tensor_rank_1([
133                -(1.0 + xi_2) * (1.0 - xi_3) / 8.0,
134                (1.0 - xi_1) * (1.0 - xi_3) / 8.0,
135                -(1.0 - xi_1) * (1.0 + xi_2) / 8.0,
136            ]),
137            tensor_rank_1([
138                -(1.0 - xi_2) * (1.0 + xi_3) / 8.0,
139                -(1.0 - xi_1) * (1.0 + xi_3) / 8.0,
140                (1.0 - xi_1) * (1.0 - xi_2) / 8.0,
141            ]),
142            tensor_rank_1([
143                (1.0 - xi_2) * (1.0 + xi_3) / 8.0,
144                -(1.0 + xi_1) * (1.0 + xi_3) / 8.0,
145                (1.0 + xi_1) * (1.0 - xi_2) / 8.0,
146            ]),
147            tensor_rank_1([
148                (1.0 + xi_2) * (1.0 + xi_3) / 8.0,
149                (1.0 + xi_1) * (1.0 + xi_3) / 8.0,
150                (1.0 + xi_1) * (1.0 + xi_2) / 8.0,
151            ]),
152            tensor_rank_1([
153                -(1.0 + xi_2) * (1.0 + xi_3) / 8.0,
154                (1.0 - xi_1) * (1.0 + xi_3) / 8.0,
155                (1.0 - xi_1) * (1.0 + xi_2) / 8.0,
156            ]),
157        ])
158    }
159    const fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
160        StandardGradientOperators::const_from([
161            Self::shape_functions_gradients(Self::integration_point(0)),
162            Self::shape_functions_gradients(Self::integration_point(1)),
163            Self::shape_functions_gradients(Self::integration_point(2)),
164            Self::shape_functions_gradients(Self::integration_point(3)),
165            Self::shape_functions_gradients(Self::integration_point(4)),
166            Self::shape_functions_gradients(Self::integration_point(5)),
167            Self::shape_functions_gradients(Self::integration_point(6)),
168            Self::shape_functions_gradients(Self::integration_point(7)),
169        ])
170    }
171}