Skip to main content

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

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