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

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