Skip to main content

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

1#[cfg(test)]
2mod test;
3
4pub mod cohesive;
5pub mod linear;
6pub mod planar;
7pub mod quadratic;
8pub mod serendipity;
9pub mod solid;
10pub mod surface;
11pub mod thermal;
12
13use crate::{
14    domain::block::element::{ElementError, ElementKind},
15    math::{
16        Projection, Quantity, Reference, Scalar, ScalarList, TensorList, TensorRank1,
17        TensorRank1List, TensorRank1List2D,
18    },
19    mechanics::{CoordinateList, CurrentCoordinates, CurrentVelocities, ReferenceCoordinates},
20    units::{Length, ReciprocalLength, Volume},
21};
22use std::fmt::{self, Debug, Formatter};
23
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> = CurrentVelocities<N>;
29pub type ElementNodalEitherCoordinates<I, 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, Reference, N, G, ReciprocalLength>;
33pub type IntegrationWeights<const G: usize, U> = TensorList<Quantity<U>, G>;
34pub type ParametricCoordinate<const M: usize> = TensorRank1<M, Projection>;
35pub type ParametricCoordinates<const G: usize, const M: usize> = TensorRank1List<M, Projection, G>;
36pub type ParametricReference<const M: usize, const N: usize> = TensorRank1List<M, Projection, N>;
37pub type ShapeFunctions<const N: usize> = TensorRank1<N, Projection>;
38pub type ShapeFunctionsAtIntegrationPoints<const G: usize, const N: usize> =
39    TensorRank1List<N, Projection, G>;
40pub type ShapeFunctionsGradients<const M: usize, const N: usize> = TensorRank1List<M, Reference, N>;
41pub type StandardGradientOperators<const M: usize, const O: usize, const P: usize> =
42    TensorRank1List2D<M, Reference, O, P>;
43pub type StandardGradientOperatorsTransposed<const M: usize, const O: usize, const P: usize> =
44    TensorRank1List2D<M, Reference, P, O>;
45
46pub trait FiniteElement<const G: usize, const M: usize, const N: usize, const P: usize, W = Volume>
47where
48    Self: Clone + Debug,
49{
50    fn integration_points() -> ParametricCoordinates<G, M>;
51    fn integration_weights(&self) -> &IntegrationWeights<G, W>;
52    fn parametric_reference() -> ParametricReference<M, N>;
53    fn parametric_weights() -> ScalarList<G>;
54    fn shape_functions(parametric_coordinate: ParametricCoordinate<M>) -> ShapeFunctions<P>;
55    fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, P> {
56        Self::integration_points()
57            .into_iter()
58            .map(|integration_point| Self::shape_functions(integration_point))
59            .collect()
60    }
61    fn shape_functions_gradients(
62        parametric_coordinate: ParametricCoordinate<M>,
63    ) -> ShapeFunctionsGradients<M, P>;
64    fn shape_functions_gradients_at_integration_points() -> StandardGradientOperators<M, P, G> {
65        Self::integration_points()
66            .into_iter()
67            .map(|integration_point| Self::shape_functions_gradients(integration_point))
68            .collect()
69    }
70    fn volume(&self) -> Quantity<W> {
71        self.integration_weights().into_iter().sum()
72    }
73}
74
75#[derive(Clone)]
76pub struct Element<const D: usize, const G: usize, const N: usize, const O: usize> {
77    gradient_vectors: GradientVectors<D, G, N>,
78    integration_weights: IntegrationWeights<G, Volume>,
79}
80
81impl<const D: usize, const G: usize, const N: usize, const O: usize> Element<D, G, N, O> {
82    pub(crate) fn gradient_vectors(&self) -> &GradientVectors<D, G, N> {
83        &self.gradient_vectors
84    }
85}
86
87impl<const D: usize, const G: usize, const N: usize, const O: usize> Debug for Element<D, G, N, O> {
88    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
89        let element = match (D, G, N, O) {
90            (2, 1, 3, 1) => "LinearTriangle",
91            (2, 4, 4, 1) => "LinearQuadrilateral",
92            (3, 8, 8, 1) => "LinearHexahedron",
93            (3, 8, 5, 1) => "LinearPyramid",
94            (3, 1, 4, 1) => "LinearTetrahedron",
95            (3, 6, 6, 1) => "LinearWedge",
96            (3, 27, 27, 2) => "QuadraticHexahedron",
97            (3, 4, 10, 2) => "QuadraticTetrahedron",
98            (3, 27, 13, 2) => "QuadraticPyramid",
99            (3, 18, 15, 2) => "QuadraticWedge",
100            (3, 27, 20, 2) => "SerendipityHexahedron",
101            (3, 4, 10, 0) => "CompositeTetrahedron",
102            _ => panic!(),
103        };
104        write!(f, "{element} {{ integration points: {G}, nodes: {N} }}",)
105    }
106}
107
108fn basic_from<const D: usize, const G: usize, const N: usize, const O: usize>(
109    reference_nodal_coordinates: TensorRank1List<D, Reference, N, Length>,
110) -> Element<D, G, N, O>
111where
112    Element<D, G, N, O>: FiniteElement<G, D, N, N>,
113{
114    let gradient_vectors = Element::shape_functions_gradients_at_integration_points()
115        .into_iter()
116        .map(|standard_gradient_operator| {
117            (&reference_nodal_coordinates * &standard_gradient_operator).inverse_transpose()
118                * standard_gradient_operator
119        })
120        .collect();
121    let integration_weights = Element::shape_functions_gradients_at_integration_points()
122        .into_iter()
123        .zip(Element::parametric_weights())
124        .map(|(standard_gradient_operator, integration_weight)| {
125            Quantity::new(
126                (&reference_nodal_coordinates * standard_gradient_operator).determinant()
127                    * integration_weight,
128            )
129        })
130        .collect();
131    Element {
132        gradient_vectors,
133        integration_weights,
134    }
135}
136
137pub struct FiniteElementKind;
138
139impl ElementKind for FiniteElementKind {
140    const NAME: &'static str = "finite element";
141}
142
143pub type FiniteElementError = ElementError<FiniteElementKind>;