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

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

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