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

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