conspire/domain/fem/block/element/
mod.rs

1#[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; // nightly feature
22const 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: Clone + Debug,
44{
45    fn integration_points() -> ParametricCoordinates<G, M>;
46    fn integration_weights(&self) -> &ScalarList<G>;
47    fn minimum_scaled_jacobian<const I: usize>(
48        nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
49    ) -> Scalar {
50        Self::scaled_jacobians(nodal_coordinates)
51            .into_iter()
52            .reduce(Scalar::min)
53            .unwrap()
54    }
55    fn parametric_reference() -> ParametricReference<M, N>;
56    fn parametric_weights() -> ScalarList<G>;
57    fn scaled_jacobians<const I: usize>(
58        nodal_coordinates: ElementNodalEitherCoordinates<I, N>,
59    ) -> ScalarList<P>;
60    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<P>;
61    fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, P> {
62        Self::integration_points()
63            .into_iter()
64            .map(|integration_point| Self::shape_functions(integration_point))
65            .collect()
66    }
67    fn shape_functions_gradients(
68        parametric_coordinate: ParametricCoordinate<M>,
69    ) -> ShapeFunctionsGradients<M, P>;
70    fn shape_functions_gradients_at_integration_points() -> StandardGradientOperators<M, P, G> {
71        Self::integration_points()
72            .into_iter()
73            .map(|integration_point| Self::shape_functions_gradients(integration_point))
74            .collect()
75    }
76    fn volume(&self) -> Scalar {
77        self.integration_weights().into_iter().sum()
78    }
79}
80
81#[derive(Clone)]
82pub struct Element<const G: usize, const N: usize, const O: usize> {
83    gradient_vectors: GradientVectors<G, N>,
84    integration_weights: ScalarList<G>,
85}
86
87impl<const G: usize, const N: usize, const O: usize> Element<G, N, O> {
88    fn gradient_vectors(&self) -> &GradientVectors<G, N> {
89        &self.gradient_vectors
90    }
91}
92
93impl<const G: usize, const N: usize, const O: usize> Debug for Element<G, N, O> {
94    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
95        let element = match (G, N, O) {
96            (8, 8, 1) => "LinearHexahedron",
97            (8, 5, 1) => "LinearPyramid",
98            (1, 4, 1) => "LinearTetrahedron",
99            (6, 6, 1) => "LinearWedge",
100            (27, 27, 2) => "QuadraticHexahedron",
101            (4, 10, 2) => "QuadraticTetrahedron",
102            (27, 13, 2) => "QuadraticPyramid",
103            (18, 15, 2) => "QuadraticWedge",
104            (27, 20, 2) => "SerendipityHexahedron",
105            (4, 10, 0) => "CompositeTetrahedron",
106            _ => panic!(),
107        };
108        write!(f, "{element} {{ integration points: {G}, nodes: {N} }}",)
109    }
110}
111
112impl<const G: usize, const N: usize, const O: usize> Default for Element<G, N, O>
113where
114    Self: FiniteElement<G, 3, N, N> + From<ElementNodalReferenceCoordinates<N>>,
115{
116    fn default() -> Self {
117        ElementNodalReferenceCoordinates::from(Self::parametric_reference()).into()
118    }
119}
120
121fn basic_from<const G: usize, const N: usize, const O: usize>(
122    reference_nodal_coordinates: ElementNodalReferenceCoordinates<N>,
123) -> Element<G, N, O>
124where
125    Element<G, N, O>: FiniteElement<G, 3, N, N>,
126{
127    let gradient_vectors = Element::shape_functions_gradients_at_integration_points()
128        .into_iter()
129        .map(|standard_gradient_operator| {
130            (&reference_nodal_coordinates * &standard_gradient_operator).inverse_transpose()
131                * standard_gradient_operator
132        })
133        .collect();
134    let integration_weights = Element::shape_functions_gradients_at_integration_points()
135        .into_iter()
136        .zip(Element::parametric_weights())
137        .map(|(standard_gradient_operator, integration_weight)| {
138            (&reference_nodal_coordinates * standard_gradient_operator).determinant()
139                * integration_weight
140        })
141        .collect();
142    Element {
143        gradient_vectors,
144        integration_weights,
145    }
146}
147
148pub enum FiniteElementError {
149    Upstream(String, String),
150}
151
152impl From<FiniteElementError> for TestError {
153    fn from(error: FiniteElementError) -> Self {
154        Self {
155            message: error.to_string(),
156        }
157    }
158}
159
160impl Debug for FiniteElementError {
161    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
162        let error = match self {
163            Self::Upstream(error, element) => {
164                format!(
165                    "{error}\x1b[0;91m\n\
166                    In finite element: {element}."
167                )
168            }
169        };
170        write!(f, "\n{error}\n\x1b[0;2;31m{}\x1b[0m\n", defeat_message())
171    }
172}
173
174impl Display for FiniteElementError {
175    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
176        let error = match self {
177            Self::Upstream(error, element) => {
178                format!(
179                    "{error}\x1b[0;91m\n\
180                    In finite element: {element}."
181                )
182            }
183        };
184        write!(f, "{error}\x1b[0m")
185    }
186}
conspire/domain/fem/block/element/mod.rs

conspire/domain/fem/block/element/
mod.rs
