conspire/domain/fem/block/element/
mod.rs1#[cfg(test)]
2mod test;
3
4pub mod cohesive;
5pub mod composite;
6pub mod linear;
7pub mod quadratic;
8pub mod serendipity;
9pub mod solid;
10pub mod surface;
11pub mod thermal;
12
13use crate::{
14 defeat_message,
15 math::{Scalar, ScalarList, TensorRank1, TensorRank1List, TensorRank1List2D, TestError},
16 mechanics::{CoordinateList, CurrentCoordinates, ReferenceCoordinates, VectorList2D},
17};
18use std::fmt::{self, Debug, Display, Formatter};
19
20const A: usize = 9;
21const FRAC_1_SQRT_3: Scalar = 0.577_350_269_189_625_8; const FRAC_SQRT_3_5: Scalar = 0.774_596_669_241_483;
23
24pub type ElementNodalCoordinates<const N: usize> = CurrentCoordinates<N>;
25pub type ElementNodalVelocities<const N: usize> = CurrentCoordinates<N>;
26pub type ElementNodalEitherCoordinates<const I: usize, const N: usize> = CoordinateList<I, N>;
27pub type ElementNodalReferenceCoordinates<const N: usize> = ReferenceCoordinates<N>;
28pub type GradientVectors<const G: usize, const N: usize> = VectorList2D<0, N, G>;
29pub type ParametricCoordinate<const M: usize> = TensorRank1<M, A>;
30pub type ParametricCoordinates<const G: usize, const M: usize> = TensorRank1List<M, A, G>;
31pub type ParametricReference<const M: usize, const N: usize> = TensorRank1List<M, A, N>;
32pub type ShapeFunctions<const N: usize> = TensorRank1<N, A>;
33pub type ShapeFunctionsAtIntegrationPoints<const G: usize, const N: usize> =
34 TensorRank1List<N, A, G>;
35pub type ShapeFunctionsGradients<const M: usize, const N: usize> = TensorRank1List<M, 0, N>;
36pub type StandardGradientOperators<const M: usize, const O: usize, const P: usize> =
37 TensorRank1List2D<M, 0, O, P>;
38pub type StandardGradientOperatorsTransposed<const M: usize, const O: usize, const P: usize> =
39 TensorRank1List2D<M, 0, P, O>;
40
41pub trait FiniteElement<const G: usize, const M: usize, const N: usize, const P: usize>
42where
43 Self: Debug,
44{
45 fn integration_points() -> ParametricCoordinates<G, M>;
46 fn integration_weights(&self) -> &ScalarList<G>;
47 fn parametric_reference() -> ParametricReference<M, N>;
48 fn parametric_weights() -> ScalarList<G>;
49 fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<P>;
50 fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, P> {
51 Self::integration_points()
52 .into_iter()
53 .map(|integration_point| Self::shape_functions(integration_point))
54 .collect()
55 }
56 fn shape_functions_gradients(
57 parametric_coordinate: ParametricCoordinate<M>,
58 ) -> ShapeFunctionsGradients<M, P>;
59 fn shape_functions_gradients_at_integration_points() -> StandardGradientOperators<M, P, G> {
60 Self::integration_points()
61 .into_iter()
62 .map(|integration_point| Self::shape_functions_gradients(integration_point))
63 .collect()
64 }
65 fn volume(&self) -> Scalar {
66 self.integration_weights().into_iter().sum()
67 }
68}
69
70pub struct Element<const G: usize, const N: usize, const O: usize> {
71 gradient_vectors: GradientVectors<G, N>,
72 integration_weights: ScalarList<G>,
73}
74
75impl<const G: usize, const N: usize, const O: usize> Element<G, N, O> {
76 fn gradient_vectors(&self) -> &GradientVectors<G, N> {
77 &self.gradient_vectors
78 }
79}
80
81impl<const G: usize, const N: usize, const O: usize> Debug for Element<G, N, O> {
82 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
83 let element = match (G, N, O) {
84 (8, 8, 1) => "LinearHexahedron",
85 (8, 5, 1) => "LinearPyramid",
86 (1, 4, 1) => "LinearTetrahedron",
87 (6, 6, 1) => "LinearWedge",
88 (27, 27, 2) => "QuadraticHexahedron",
89 (4, 10, 2) => "QuadraticTetrahedron",
90 (27, 13, 2) => "QuadraticPyramid",
91 (18, 15, 2) => "QuadraticWedge",
92 (27, 20, 2) => "SerendipityHexahedron",
93 (4, 10, 0) => "CompositeTetrahedron",
94 _ => panic!(),
95 };
96 write!(f, "{element} {{ integration points: {G}, nodes: {N} }}",)
97 }
98}
99
100impl<const G: usize, const N: usize, const O: usize> Default for Element<G, N, O>
101where
102 Self: FiniteElement<G, 3, N, N> + From<ElementNodalReferenceCoordinates<N>>,
103{
104 fn default() -> Self {
105 ElementNodalReferenceCoordinates::from(Self::parametric_reference()).into()
106 }
107}
108
109fn basic_from<const G: usize, const N: usize, const O: usize>(
110 reference_nodal_coordinates: ElementNodalReferenceCoordinates<N>,
111) -> Element<G, N, O>
112where
113 Element<G, N, O>: FiniteElement<G, 3, N, N>,
114{
115 let gradient_vectors = Element::shape_functions_gradients_at_integration_points()
116 .into_iter()
117 .map(|standard_gradient_operator| {
118 (&reference_nodal_coordinates * &standard_gradient_operator).inverse_transpose()
119 * standard_gradient_operator
120 })
121 .collect();
122 let integration_weights = Element::shape_functions_gradients_at_integration_points()
123 .into_iter()
124 .zip(Element::parametric_weights())
125 .map(|(standard_gradient_operator, integration_weight)| {
126 (&reference_nodal_coordinates * standard_gradient_operator).determinant()
127 * integration_weight
128 })
129 .collect();
130 Element {
131 gradient_vectors,
132 integration_weights,
133 }
134}
135
136pub enum FiniteElementError {
137 Upstream(String, String),
138}
139
140impl From<FiniteElementError> for TestError {
141 fn from(error: FiniteElementError) -> Self {
142 Self {
143 message: error.to_string(),
144 }
145 }
146}
147
148impl Debug for FiniteElementError {
149 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
150 let error = match self {
151 Self::Upstream(error, element) => {
152 format!(
153 "{error}\x1b[0;91m\n\
154 In finite element: {element}."
155 )
156 }
157 };
158 write!(f, "\n{error}\n\x1b[0;2;31m{}\x1b[0m\n", defeat_message())
159 }
160}
161
162impl Display for FiniteElementError {
163 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
164 let error = match self {
165 Self::Upstream(error, element) => {
166 format!(
167 "{error}\x1b[0;91m\n\
168 In finite element: {element}."
169 )
170 }
171 };
172 write!(f, "{error}\x1b[0m")
173 }
174}