runmat_analysis_fea/assembly/elements/solid/
tetrahedron10.rs1use thiserror::Error;
2
3use runmat_meshing_core::TETRAHEDRON_MIDSIDE_EDGE_CORNERS;
4
5use super::{
6 material::{elasticity_matrix, SolidMaterialError},
7 SolidMaterial,
8};
9
10pub const TETRAHEDRON10_ELEMENT_NODE_COUNT: usize = 10;
11pub const TETRAHEDRON10_NODE_DOF_COUNT: usize = 3;
12pub const TETRAHEDRON10_ELEMENT_DOF_COUNT: usize =
13 TETRAHEDRON10_ELEMENT_NODE_COUNT * TETRAHEDRON10_NODE_DOF_COUNT;
14
15pub type Tetrahedron10Matrix30 =
16 [[f64; TETRAHEDRON10_ELEMENT_DOF_COUNT]; TETRAHEDRON10_ELEMENT_DOF_COUNT];
17
18const BARYCENTRIC_DERIVATIVES: [[f64; 3]; 4] = [
19 [-1.0, -1.0, -1.0],
20 [1.0, 0.0, 0.0],
21 [0.0, 1.0, 0.0],
22 [0.0, 0.0, 1.0],
23];
24const QUADRATURE_A: f64 = 0.585_410_196_624_968_5;
25const QUADRATURE_B: f64 = 0.138_196_601_125_010_5;
26const QUADRATURE_WEIGHT: f64 = 1.0 / 24.0;
27const QUADRATURE: [[f64; 4]; 4] = [
28 [QUADRATURE_A, QUADRATURE_B, QUADRATURE_B, QUADRATURE_B],
29 [QUADRATURE_B, QUADRATURE_A, QUADRATURE_B, QUADRATURE_B],
30 [QUADRATURE_B, QUADRATURE_B, QUADRATURE_A, QUADRATURE_B],
31 [QUADRATURE_B, QUADRATURE_B, QUADRATURE_B, QUADRATURE_A],
32];
33
34#[derive(Debug, Clone, Copy, PartialEq)]
35pub struct Tetrahedron10ElementGeometry {
36 pub nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
38}
39
40#[derive(Debug, Error, Clone, PartialEq)]
41pub enum Tetrahedron10ElementError {
42 #[error("Tetrahedron10 node coordinates must be finite")]
43 NonFiniteCoordinate,
44 #[error("Tetrahedron10 quadrature Jacobian must be positive and finite")]
45 DegenerateOrInverted,
46 #[error("Tetrahedron10 material is invalid: {0}")]
47 InvalidMaterial(#[from] SolidMaterialError),
48}
49
50pub fn global_stiffness_matrix(
51 material: SolidMaterial,
52 geometry: Tetrahedron10ElementGeometry,
53) -> Result<Tetrahedron10Matrix30, Tetrahedron10ElementError> {
54 validate_nodes(geometry.nodes_m)?;
55 let elasticity = elasticity_matrix(material)?;
56 let mut stiffness = [[0.0; TETRAHEDRON10_ELEMENT_DOF_COUNT]; TETRAHEDRON10_ELEMENT_DOF_COUNT];
57 for barycentric in QUADRATURE {
58 let gradients = physical_shape_gradients(geometry.nodes_m, barycentric)?;
59 let b = strain_displacement_matrix(gradients);
60 let determinant = jacobian_determinant(geometry.nodes_m, barycentric)?;
61 let scale = QUADRATURE_WEIGHT * determinant;
62 let mut elasticity_b = [[0.0; TETRAHEDRON10_ELEMENT_DOF_COUNT]; 6];
63 for row in 0..6 {
64 for column in 0..TETRAHEDRON10_ELEMENT_DOF_COUNT {
65 elasticity_b[row][column] = (0..6)
66 .map(|inner| elasticity[row][inner] * b[inner][column])
67 .sum();
68 }
69 }
70 for row in 0..TETRAHEDRON10_ELEMENT_DOF_COUNT {
71 for column in row..TETRAHEDRON10_ELEMENT_DOF_COUNT {
72 let contribution = scale
73 * (0..6)
74 .map(|inner| b[inner][row] * elasticity_b[inner][column])
75 .sum::<f64>();
76 stiffness[row][column] += contribution;
77 if row != column {
78 stiffness[column][row] += contribution;
79 }
80 }
81 }
82 }
83 if stiffness.iter().flatten().all(|value| value.is_finite()) {
84 Ok(stiffness)
85 } else {
86 Err(Tetrahedron10ElementError::DegenerateOrInverted)
87 }
88}
89
90fn physical_shape_gradients(
91 nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
92 barycentric: [f64; 4],
93) -> Result<[[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT], Tetrahedron10ElementError> {
94 let reference = reference_shape_gradients(barycentric);
95 let inverse = inverse_jacobian(nodes_m, reference)?;
96 Ok(reference.map(|gradient| {
97 std::array::from_fn(|axis| {
98 inverse[0][axis] * gradient[0]
99 + inverse[1][axis] * gradient[1]
100 + inverse[2][axis] * gradient[2]
101 })
102 }))
103}
104
105fn reference_shape_gradients(
106 barycentric: [f64; 4],
107) -> [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT] {
108 let mut gradients = [[0.0; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT];
109 for node in 0..4 {
110 let factor = 4.0 * barycentric[node] - 1.0;
111 gradients[node] = BARYCENTRIC_DERIVATIVES[node].map(|derivative| factor * derivative);
112 }
113 for (local_edge, [left, right]) in TETRAHEDRON_MIDSIDE_EDGE_CORNERS.into_iter().enumerate() {
114 gradients[4 + local_edge] = std::array::from_fn(|axis| {
115 4.0 * (barycentric[right] * BARYCENTRIC_DERIVATIVES[left][axis]
116 + barycentric[left] * BARYCENTRIC_DERIVATIVES[right][axis])
117 });
118 }
119 gradients
120}
121
122fn strain_displacement_matrix(
123 gradients: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
124) -> [[f64; TETRAHEDRON10_ELEMENT_DOF_COUNT]; 6] {
125 let mut b = [[0.0; TETRAHEDRON10_ELEMENT_DOF_COUNT]; 6];
126 for (node, [dn_dx, dn_dy, dn_dz]) in gradients.into_iter().enumerate() {
127 let column = node * TETRAHEDRON10_NODE_DOF_COUNT;
128 b[0][column] = dn_dx;
129 b[1][column + 1] = dn_dy;
130 b[2][column + 2] = dn_dz;
131 b[3][column + 1] = dn_dz;
132 b[3][column + 2] = dn_dy;
133 b[4][column] = dn_dz;
134 b[4][column + 2] = dn_dx;
135 b[5][column] = dn_dy;
136 b[5][column + 1] = dn_dx;
137 }
138 b
139}
140
141fn jacobian_determinant(
142 nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
143 barycentric: [f64; 4],
144) -> Result<f64, Tetrahedron10ElementError> {
145 let reference = reference_shape_gradients(barycentric);
146 let jacobian = jacobian(nodes_m, reference);
147 let determinant = dot(jacobian[0], cross(jacobian[1], jacobian[2]));
148 if determinant.is_finite() && determinant > 0.0 {
149 Ok(determinant)
150 } else {
151 Err(Tetrahedron10ElementError::DegenerateOrInverted)
152 }
153}
154
155fn inverse_jacobian(
156 nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
157 reference: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
158) -> Result<[[f64; 3]; 3], Tetrahedron10ElementError> {
159 let matrix = jacobian(nodes_m, reference);
160 let determinant = dot(matrix[0], cross(matrix[1], matrix[2]));
161 if !determinant.is_finite() || determinant <= 0.0 {
162 return Err(Tetrahedron10ElementError::DegenerateOrInverted);
163 }
164 let inverse_determinant = 1.0 / determinant;
165 Ok([
166 [
167 (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) * inverse_determinant,
168 (matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2]) * inverse_determinant,
169 (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1]) * inverse_determinant,
170 ],
171 [
172 (matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2]) * inverse_determinant,
173 (matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0]) * inverse_determinant,
174 (matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2]) * inverse_determinant,
175 ],
176 [
177 (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * inverse_determinant,
178 (matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1]) * inverse_determinant,
179 (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]) * inverse_determinant,
180 ],
181 ])
182}
183
184fn jacobian(
185 nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
186 reference: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
187) -> [[f64; 3]; 3] {
188 std::array::from_fn(|reference_axis| {
189 std::array::from_fn(|physical_axis| {
190 (0..TETRAHEDRON10_ELEMENT_NODE_COUNT)
191 .map(|node| nodes_m[node][physical_axis] * reference[node][reference_axis])
192 .sum()
193 })
194 })
195}
196
197fn validate_nodes(
198 nodes_m: [[f64; 3]; TETRAHEDRON10_ELEMENT_NODE_COUNT],
199) -> Result<(), Tetrahedron10ElementError> {
200 if nodes_m.iter().flatten().all(|value| value.is_finite()) {
201 Ok(())
202 } else {
203 Err(Tetrahedron10ElementError::NonFiniteCoordinate)
204 }
205}
206
207fn dot(left: [f64; 3], right: [f64; 3]) -> f64 {
208 left[0] * right[0] + left[1] * right[1] + left[2] * right[2]
209}
210
211fn cross(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
212 [
213 left[1] * right[2] - left[2] * right[1],
214 left[2] * right[0] - left[0] * right[2],
215 left[0] * right[1] - left[1] * right[0],
216 ]
217}
218
219#[cfg(test)]
220mod tests {
221 use super::*;
222
223 fn straight_unit_tetrahedron() -> Tetrahedron10ElementGeometry {
224 let corners = [
225 [0.0, 0.0, 0.0],
226 [1.0, 0.0, 0.0],
227 [0.0, 1.0, 0.0],
228 [0.0, 0.0, 1.0],
229 ];
230 Tetrahedron10ElementGeometry {
231 nodes_m: [
232 corners[0],
233 corners[1],
234 corners[2],
235 corners[3],
236 midpoint(corners[0], corners[1]),
237 midpoint(corners[1], corners[2]),
238 midpoint(corners[2], corners[0]),
239 midpoint(corners[0], corners[3]),
240 midpoint(corners[1], corners[3]),
241 midpoint(corners[2], corners[3]),
242 ],
243 }
244 }
245
246 fn steel() -> SolidMaterial {
247 SolidMaterial {
248 youngs_modulus_pa: 200.0e9,
249 poisson_ratio: 0.3,
250 }
251 }
252
253 #[test]
254 fn quadratic_shape_gradients_preserve_partition_of_unity() {
255 for barycentric in QUADRATURE {
256 let gradients = reference_shape_gradients(barycentric);
257 for axis in 0..3 {
258 assert!(
259 gradients
260 .iter()
261 .map(|gradient| gradient[axis])
262 .sum::<f64>()
263 .abs()
264 < 1.0e-14
265 );
266 }
267 }
268 }
269
270 #[test]
271 fn straight_tetrahedron10_stiffness_is_symmetric_and_rejects_rigid_translation() {
272 let stiffness = global_stiffness_matrix(steel(), straight_unit_tetrahedron()).unwrap();
273 for row in 0..TETRAHEDRON10_ELEMENT_DOF_COUNT {
274 for column in 0..TETRAHEDRON10_ELEMENT_DOF_COUNT {
275 assert!((stiffness[row][column] - stiffness[column][row]).abs() < 1.0e-4);
276 }
277 }
278 let mut translation = [0.0; TETRAHEDRON10_ELEMENT_DOF_COUNT];
279 for displacement in translation.chunks_exact_mut(3) {
280 displacement.copy_from_slice(&[1.0, -2.0, 3.0]);
281 }
282 let residual = std::array::from_fn::<_, TETRAHEDRON10_ELEMENT_DOF_COUNT, _>(|row| {
283 stiffness[row]
284 .iter()
285 .zip(translation)
286 .map(|(left, right)| left * right)
287 .sum::<f64>()
288 });
289 assert!(residual.iter().all(|value| value.abs() < 2.0e-3));
290 }
291
292 #[test]
293 fn tetrahedron10_rejects_nonfinite_and_inverted_geometry() {
294 let mut nonfinite = straight_unit_tetrahedron();
295 nonfinite.nodes_m[4][0] = f64::NAN;
296 assert_eq!(
297 global_stiffness_matrix(steel(), nonfinite).unwrap_err(),
298 Tetrahedron10ElementError::NonFiniteCoordinate
299 );
300
301 let mut inverted = straight_unit_tetrahedron();
302 inverted.nodes_m.swap(1, 2);
303 inverted.nodes_m.swap(4, 6);
304 inverted.nodes_m.swap(8, 9);
305 assert_eq!(
306 global_stiffness_matrix(steel(), inverted).unwrap_err(),
307 Tetrahedron10ElementError::DegenerateOrInverted
308 );
309 }
310
311 fn midpoint(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
312 std::array::from_fn(|axis| (left[axis] + right[axis]) * 0.5)
313 }
314}