conspire/domain/fem/block/element/surface/
mod.rs1pub mod linear;
2
3use crate::{
4 fem::block::element::{
5 ElementNodalCoordinates, ElementNodalEitherCoordinates, ElementNodalReferenceCoordinates,
6 ElementNodalVelocities, FiniteElement, GradientVectors, IntegrationWeights,
7 },
8 math::{CrossProduct, IDENTITY, LEVI_CIVITA, Quantity, Tensor, TensorArray, TensorRank2},
9 mechanics::{
10 Normal, NormalGradients, NormalRates, Normals, ReferenceNormals, SurfaceBases,
11 SurfaceDualBases,
12 },
13 units::{Area, Length, Rate, Volume},
14};
15use std::fmt::{self, Debug, Formatter};
16
17const M: usize = 2;
18
19#[derive(Clone)]
20pub struct SurfaceElement<const G: usize, const N: usize, const O: usize> {
21 gradient_vectors: GradientVectors<3, G, N>,
22 integration_weights: IntegrationWeights<G, Volume>,
23 reference_normals: ReferenceNormals<G>,
24}
25
26impl<const G: usize, const N: usize, const O: usize> SurfaceElement<G, N, O> {
27 pub fn gradient_vectors(&self) -> &GradientVectors<3, G, N> {
28 &self.gradient_vectors
29 }
30 pub fn reference_normals(&self) -> &ReferenceNormals<G> {
31 &self.reference_normals
32 }
33}
34
35impl<const G: usize, const N: usize, const O: usize> Debug for SurfaceElement<G, N, O> {
36 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
37 let element = match (G, N, O) {
38 (1, 3, 1) => "LinearTriangle",
39 (4, 4, 1) => "LinearQuadrilateral",
40 _ => panic!(),
41 };
42 write!(f, "{element} {{ G: {G}, N: {N} }}",)
43 }
44}
45
46pub trait SurfaceFiniteElement<const G: usize, const N: usize, const P: usize, W = Volume>
47where
48 Self: FiniteElement<G, M, N, P, W>,
49{
50 fn bases<I>(nodal_coordinates: &ElementNodalEitherCoordinates<I, P>) -> 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<I>(
70 nodal_coordinates: &ElementNodalEitherCoordinates<I, P>,
71 ) -> SurfaceDualBases<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, Area>>()
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| {
103 let normal = basis_vectors[0].cross(&basis_vectors[1]);
104 let area = normal.norm();
105 normal / area
106 })
107 .collect()
108 }
109 fn normal_gradients(nodal_coordinates: &ElementNodalCoordinates<P>) -> NormalGradients<P, G> {
110 let levi_civita_symbol = LEVI_CIVITA;
111 let mut normalization = Quantity::<Area>::new(0.0);
112 let mut normal_vector = Normal::zero();
113 Self::shape_functions_gradients_at_integration_points().iter()
114 .zip(Self::bases(nodal_coordinates))
115 .map(|(standard_gradient_operator, basis_vectors)|{
116 normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
117 normal_vector = basis_vectors[0].cross(&basis_vectors[1]) / normalization;
118 standard_gradient_operator.iter()
119 .map(|standard_gradient_operator_a|
120 levi_civita_symbol.iter()
121 .map(|levi_civita_symbol_m|
122 IDENTITY.iter()
123 .zip(normal_vector.iter())
124 .map(|(identity_i, normal_vector_i)|
125 levi_civita_symbol_m.iter()
126 .zip(basis_vectors[0].iter()
127 .zip(basis_vectors[1].iter()))
128 .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
129 levi_civita_symbol_mn.iter()
130 .zip(identity_i.iter()
131 .zip(normal_vector.iter()))
132 .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
133 levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
134 ).sum::<Quantity>() * (
135 standard_gradient_operator_a[0] * basis_vector_1_n
136 - standard_gradient_operator_a[1] * basis_vector_0_n
137 )
138 ).sum::<Quantity<Length>>()
139 / normalization
140 ).collect()
141 ).collect()
142 ).collect()
143 }).collect()
144 }
145 fn normal_rates(
146 nodal_coordinates: &ElementNodalCoordinates<P>,
147 nodal_velocities: &ElementNodalVelocities<P>,
148 ) -> NormalRates<G> {
149 let identity = IDENTITY;
150 let levi_civita_symbol = LEVI_CIVITA;
151 let mut normalization = Quantity::<Area>::new(0.0);
152 Self::bases(nodal_coordinates)
153 .iter()
154 .zip(Self::normals(nodal_coordinates).iter()
155 .zip(Self::shape_functions_gradients_at_integration_points()))
156 .map(|(basis, (normal, standard_gradient_operator))| {
157 normalization = basis[0].cross(&basis[1]).norm();
158 identity.iter()
159 .zip(normal.iter())
160 .map(|(identity_i, normal_vector_i)|
161 nodal_velocities.iter()
162 .zip(standard_gradient_operator.iter())
163 .map(|(nodal_velocity_a, standard_gradient_operator_a)|
164 levi_civita_symbol.iter()
165 .zip(nodal_velocity_a.iter())
166 .map(|(levi_civita_symbol_m, nodal_velocity_a_m)|
167 levi_civita_symbol_m.iter()
168 .zip(basis[0].iter()
169 .zip(basis[1].iter()))
170 .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
171 levi_civita_symbol_mn.iter()
172 .zip(identity_i.iter()
173 .zip(normal.iter()))
174 .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
175 levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
176 ).sum::<Quantity>() * (
177 standard_gradient_operator_a[0] * basis_vector_1_n
178 - standard_gradient_operator_a[1] * basis_vector_0_n
179 )
180 ).sum::<Quantity<Length>>()
181 / normalization
182 * nodal_velocity_a_m
183 ).sum::<Quantity<Rate>>()
184 ).sum::<Quantity<Rate>>()
185 ).collect()
186 }).collect()
187 }
188}
189
190impl<const G: usize, const N: usize, const O: usize, const P: usize> SurfaceFiniteElement<G, N, P>
191 for SurfaceElement<G, N, O>
192where
193 Self: FiniteElement<G, M, N, P>,
194{
195}
196
197impl<const G: usize, const N: usize, const O: usize>
198 From<(ElementNodalReferenceCoordinates<N>, Quantity<Length>)> for SurfaceElement<G, N, O>
199where
200 Self: SurfaceFiniteElement<G, N, N>,
201{
202 fn from(
203 (reference_nodal_coordinates, thickness): (
204 ElementNodalReferenceCoordinates<N>,
205 Quantity<Length>,
206 ),
207 ) -> Self {
208 let integration_weights = Self::bases(&reference_nodal_coordinates)
209 .into_iter()
210 .zip(Self::parametric_weights())
211 .map(|(reference_basis, parametric_weight)| {
212 reference_basis[0].cross(&reference_basis[1]).norm() * parametric_weight * thickness
213 })
214 .collect();
215 let reference_dual_bases = Self::dual_bases(&reference_nodal_coordinates);
216 let gradient_vectors = Self::shape_functions_gradients_at_integration_points()
217 .into_iter()
218 .zip(reference_dual_bases.iter())
219 .map(|(standard_gradient_operator, reference_dual_basis)| {
220 standard_gradient_operator
221 .iter()
222 .map(|standard_gradient_operator_a| {
223 standard_gradient_operator_a
224 .iter()
225 .zip(reference_dual_basis.iter())
226 .map(|(standard_gradient_operator_a_m, reference_dual_basis_m)| {
227 reference_dual_basis_m * standard_gradient_operator_a_m
228 })
229 .sum()
230 })
231 .collect()
232 })
233 .collect();
234 let reference_normals = reference_dual_bases
235 .into_iter()
236 .map(|reference_dual_basis| {
237 let normal = reference_dual_basis[0].cross(&reference_dual_basis[1]);
238 let reciprocal_area = normal.norm();
239 normal / reciprocal_area
240 })
241 .collect();
242 Self {
243 gradient_vectors,
244 integration_weights,
245 reference_normals,
246 }
247 }
248}