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

1pub mod linear;
2
3use crate::{
4    fem::block::element::{
5        ElementNodalCoordinates, ElementNodalEitherCoordinates, ElementNodalReferenceCoordinates,
6        ElementNodalVelocities, FiniteElement, GradientVectors,
7    },
8    math::{IDENTITY, LEVI_CIVITA, Scalar, ScalarList, Tensor, TensorArray, TensorRank2},
9    mechanics::{Normal, NormalGradients, NormalRates, Normals, ReferenceNormals, SurfaceBases},
10};
11use std::fmt::{self, Debug, Formatter};
12
13const M: usize = 2;
14
15pub struct SurfaceElement<const G: usize, const N: usize, const O: usize> {
16    gradient_vectors: GradientVectors<G, N>,
17    integration_weights: ScalarList<G>,
18    reference_normals: ReferenceNormals<G>,
19}
20
21impl<const G: usize, const N: usize, const O: usize> SurfaceElement<G, N, O> {
22    pub fn gradient_vectors(&self) -> &GradientVectors<G, N> {
23        &self.gradient_vectors
24    }
25    pub fn reference_normals(&self) -> &ReferenceNormals<G> {
26        &self.reference_normals
27    }
28}
29
30impl<const G: usize, const N: usize, const O: usize> Debug for SurfaceElement<G, N, O> {
31    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
32        let element = match (G, N, O) {
33            (1, 3, 1) => "LinearTriangle",
34            (4, 4, 1) => "LinearQuadrilateral",
35            _ => panic!(),
36        };
37        write!(f, "{element} {{ G: {G}, N: {N} }}",)
38    }
39}
40
41pub trait SurfaceFiniteElement<const G: usize, const N: usize, const P: usize>
42where
43    Self: FiniteElement<G, M, N, P>,
44{
45    fn bases<const I: usize>(
46        nodal_coordinates: &ElementNodalEitherCoordinates<I, P>,
47    ) -> SurfaceBases<I, G> {
48        Self::shape_functions_gradients_at_integration_points()
49            .iter()
50            .map(|shape_functions_gradients| {
51                shape_functions_gradients
52                    .iter()
53                    .zip(nodal_coordinates)
54                    .map(|(shape_functions_gradient, nodal_coordinate)| {
55                        shape_functions_gradient
56                            .iter()
57                            .map(|standard_gradient_operator_m| {
58                                nodal_coordinate * standard_gradient_operator_m
59                            })
60                            .collect()
61                    })
62                    .sum()
63            })
64            .collect()
65    }
66    fn dual_bases<const I: usize>(
67        nodal_coordinates: &ElementNodalEitherCoordinates<I, P>,
68    ) -> SurfaceBases<I, G> {
69        Self::bases(nodal_coordinates)
70            .into_iter()
71            .map(|basis_vectors| {
72                basis_vectors
73                    .iter()
74                    .map(|basis_vector_m| {
75                        basis_vectors
76                            .iter()
77                            .map(|basis_vector_n| basis_vector_m * basis_vector_n)
78                            .collect()
79                    })
80                    .collect::<TensorRank2<2, I, I>>()
81                    .inverse()
82                    .iter()
83                    .map(|metric_tensor_m| {
84                        metric_tensor_m
85                            .iter()
86                            .zip(basis_vectors.iter())
87                            .map(|(metric_tensor_mn, basis_vectors_n)| {
88                                basis_vectors_n * metric_tensor_mn
89                            })
90                            .sum()
91                    })
92                    .collect()
93            })
94            .collect()
95    }
96    fn normals(nodal_coordinates: &ElementNodalCoordinates<P>) -> Normals<G> {
97        Self::bases(nodal_coordinates)
98            .into_iter()
99            .map(|basis_vectors| basis_vectors[0].cross(&basis_vectors[1]).normalized())
100            .collect()
101    }
102    fn normal_gradients(nodal_coordinates: &ElementNodalCoordinates<P>) -> NormalGradients<P, G> {
103        let levi_civita_symbol = LEVI_CIVITA;
104        let mut normalization: Scalar = 0.0;
105        let mut normal_vector = Normal::zero();
106        Self::shape_functions_gradients_at_integration_points().iter()
107        .zip(Self::bases(nodal_coordinates))
108        .map(|(standard_gradient_operator, basis_vectors)|{
109            normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
110            normal_vector = basis_vectors[0].cross(&basis_vectors[1])/normalization;
111            standard_gradient_operator.iter()
112            .map(|standard_gradient_operator_a|
113                levi_civita_symbol.iter()
114                .map(|levi_civita_symbol_m|
115                    IDENTITY.iter()
116                    .zip(normal_vector.iter())
117                    .map(|(identity_i, normal_vector_i)|
118                        levi_civita_symbol_m.iter()
119                        .zip(basis_vectors[0].iter()
120                        .zip(basis_vectors[1].iter()))
121                        .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
122                            levi_civita_symbol_mn.iter()
123                            .zip(identity_i.iter()
124                            .zip(normal_vector.iter()))
125                            .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
126                                levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
127                            ).sum::<Scalar>() * (
128                                standard_gradient_operator_a[0] * basis_vector_1_n
129                              - standard_gradient_operator_a[1] * basis_vector_0_n
130                            )
131                        ).sum::<Scalar>() / normalization
132                    ).collect()
133                ).collect()
134            ).collect()
135        }).collect()
136    }
137    fn normal_rates(
138        nodal_coordinates: &ElementNodalCoordinates<P>,
139        nodal_velocities: &ElementNodalVelocities<P>,
140    ) -> NormalRates<G> {
141        let identity = IDENTITY;
142        let levi_civita_symbol = LEVI_CIVITA;
143        let mut normalization = 0.0;
144        Self::bases(nodal_coordinates)
145            .iter()
146            .zip(Self::normals(nodal_coordinates).iter()
147            .zip(Self::shape_functions_gradients_at_integration_points()))
148            .map(|(basis, (normal, standard_gradient_operator))| {
149                normalization = basis[0].cross(&basis[1]).norm();
150                identity.iter()
151                .zip(normal.iter())
152                .map(|(identity_i, normal_vector_i)|
153                    nodal_velocities.iter()
154                    .zip(standard_gradient_operator.iter())
155                    .map(|(nodal_velocity_a, standard_gradient_operator_a)|
156                        levi_civita_symbol.iter()
157                        .zip(nodal_velocity_a.iter())
158                        .map(|(levi_civita_symbol_m, nodal_velocity_a_m)|
159                            levi_civita_symbol_m.iter()
160                            .zip(basis[0].iter()
161                            .zip(basis[1].iter()))
162                            .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
163                                levi_civita_symbol_mn.iter()
164                                .zip(identity_i.iter()
165                                .zip(normal.iter()))
166                                .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
167                                    levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
168                                ).sum::<Scalar>() * (
169                                    standard_gradient_operator_a[0] * basis_vector_1_n
170                                - standard_gradient_operator_a[1] * basis_vector_0_n
171                                )
172                            ).sum::<Scalar>() * nodal_velocity_a_m
173                        ).sum::<Scalar>()
174                    ).sum::<Scalar>() / normalization
175                ).collect()
176        }).collect()
177    }
178}
179
180impl<const G: usize, const N: usize, const O: usize, const P: usize> SurfaceFiniteElement<G, N, P>
181    for SurfaceElement<G, N, O>
182where
183    Self: FiniteElement<G, M, N, P>,
184{
185}
186
187impl<const G: usize, const N: usize, const O: usize>
188    From<(ElementNodalReferenceCoordinates<N>, Scalar)> for SurfaceElement<G, N, O>
189where
190    Self: SurfaceFiniteElement<G, N, N>,
191{
192    fn from(
193        (reference_nodal_coordinates, thickness): (ElementNodalReferenceCoordinates<N>, Scalar),
194    ) -> Self {
195        let integration_weights = Self::bases(&reference_nodal_coordinates)
196            .into_iter()
197            .zip(Self::parametric_weights())
198            .map(|(reference_basis, parametric_weight)| {
199                reference_basis[0].cross(&reference_basis[1]).norm() * parametric_weight * thickness
200            })
201            .collect();
202        let reference_dual_bases = Self::dual_bases(&reference_nodal_coordinates);
203        let gradient_vectors = Self::shape_functions_gradients_at_integration_points()
204            .into_iter()
205            .zip(reference_dual_bases.iter())
206            .map(|(standard_gradient_operator, reference_dual_basis)| {
207                standard_gradient_operator
208                    .iter()
209                    .map(|standard_gradient_operator_a| {
210                        standard_gradient_operator_a
211                            .iter()
212                            .zip(reference_dual_basis.iter())
213                            .map(|(standard_gradient_operator_a_m, reference_dual_basis_m)| {
214                                reference_dual_basis_m * standard_gradient_operator_a_m
215                            })
216                            .sum()
217                    })
218                    .collect()
219            })
220            .collect();
221        let reference_normals = reference_dual_bases
222            .into_iter()
223            .map(|reference_dual_basis| {
224                reference_dual_basis[0]
225                    .cross(&reference_dual_basis[1])
226                    .normalized()
227            })
228            .collect();
229        Self {
230            gradient_vectors,
231            integration_weights,
232            reference_normals,
233        }
234    }
235}
conspire/domain/fem/block/element/surface/mod.rs

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