conspire/domain/fem/block/element/linear/hexahedron/
mod.rs1#[cfg(test)]
2mod test;
3
4use crate::{
5 fem::block::element::{
6 ElementNodalEitherCoordinates, FRAC_1_SQRT_3, FiniteElement, ParametricCoordinate,
7 ParametricCoordinates, ParametricReference, ShapeFunctions, ShapeFunctionsGradients,
8 linear::{LinearElement, LinearFiniteElement, M},
9 },
10 math::{ScalarList, Tensor, TensorArray},
11 mechanics::Coordinate,
12};
13
14const G: usize = 8;
15const N: usize = 8;
16const P: usize = N;
17
18const CORNERS: [[usize; 3]; N] = [
19 [1, 3, 4],
20 [2, 0, 5],
21 [3, 1, 6],
22 [0, 2, 7],
23 [7, 5, 0],
24 [4, 6, 1],
25 [5, 7, 2],
26 [6, 4, 3],
27];
28
29pub type Hexahedron = LinearElement<G, N>;
30
31impl FiniteElement<G, M, N, P> for Hexahedron {
32 fn integration_points() -> ParametricCoordinates<G, M> {
33 [
34 [-FRAC_1_SQRT_3, -FRAC_1_SQRT_3, -FRAC_1_SQRT_3],
35 [-FRAC_1_SQRT_3, -FRAC_1_SQRT_3, FRAC_1_SQRT_3],
36 [-FRAC_1_SQRT_3, FRAC_1_SQRT_3, -FRAC_1_SQRT_3],
37 [-FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3],
38 [FRAC_1_SQRT_3, -FRAC_1_SQRT_3, -FRAC_1_SQRT_3],
39 [FRAC_1_SQRT_3, -FRAC_1_SQRT_3, FRAC_1_SQRT_3],
40 [FRAC_1_SQRT_3, FRAC_1_SQRT_3, -FRAC_1_SQRT_3],
41 [FRAC_1_SQRT_3, FRAC_1_SQRT_3, FRAC_1_SQRT_3],
42 ]
43 .into()
44 }
45 fn integration_weights(&self) -> &ScalarList<G> {
46 &self.integration_weights
47 }
48 fn parametric_reference() -> ParametricReference<M, N> {
49 [
50 [-1.0, -1.0, -1.0],
51 [1.0, -1.0, -1.0],
52 [1.0, 1.0, -1.0],
53 [-1.0, 1.0, -1.0],
54 [-1.0, -1.0, 1.0],
55 [1.0, -1.0, 1.0],
56 [1.0, 1.0, 1.0],
57 [-1.0, 1.0, 1.0],
58 ]
59 .into()
60 }
61 fn parametric_weights() -> ScalarList<G> {
62 [1.0; G].into()
63 }
64 fn scaled_jacobians<const I: usize>(
65 nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
66 ) -> ScalarList<P> {
67 let mut u = Coordinate::zero();
68 let mut v = Coordinate::zero();
69 let mut w = Coordinate::zero();
70 let mut n = Coordinate::zero();
71 CORNERS
72 .into_iter()
73 .enumerate()
74 .map(|(node, [node_a, node_b, node_c])| {
75 u = &nodal_coordinates[node_a] - &nodal_coordinates[node];
76 v = &nodal_coordinates[node_b] - &nodal_coordinates[node];
77 w = &nodal_coordinates[node_c] - &nodal_coordinates[node];
78 n = u.cross(&v);
79 (&n * &w) / u.norm() / v.norm() / w.norm()
80 })
81 .collect()
82 }
83 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<N> {
84 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
85 [
86 (1.0 - xi_1) * (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
87 (1.0 + xi_1) * (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
88 (1.0 + xi_1) * (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
89 (1.0 - xi_1) * (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
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 ]
95 .into()
96 }
97 fn shape_functions_gradients(
98 parametric_coordinate: ParametricCoordinate<M>,
99 ) -> ShapeFunctionsGradients<M, N> {
100 let [xi_1, xi_2, xi_3] = parametric_coordinate.into();
101 [
102 [
103 -(1.0 - xi_2) * (1.0 - xi_3) / 8.0,
104 -(1.0 - xi_1) * (1.0 - xi_3) / 8.0,
105 -(1.0 - xi_1) * (1.0 - xi_2) / 8.0,
106 ],
107 [
108 (1.0 - xi_2) * (1.0 - xi_3) / 8.0,
109 -(1.0 + xi_1) * (1.0 - xi_3) / 8.0,
110 -(1.0 + xi_1) * (1.0 - xi_2) / 8.0,
111 ],
112 [
113 (1.0 + xi_2) * (1.0 - xi_3) / 8.0,
114 (1.0 + xi_1) * (1.0 - xi_3) / 8.0,
115 -(1.0 + xi_1) * (1.0 + xi_2) / 8.0,
116 ],
117 [
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 [
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 [
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 [
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 [
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 ]
143 .into()
144 }
145}
146
147impl LinearFiniteElement<G, N> for Hexahedron {}