conspire/fem/block/element/linear/hexahedron/
mod.rs1#[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}