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; const 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>;