1#[cfg(test)]
2mod test;
3
4use super::*;
5use crate::{
6 math::{IDENTITY, TensorRank1List, tensor_rank_1},
7 mechanics::Scalar,
8};
9
10const G: usize = 1;
11const M: usize = 2;
12const N: usize = 3;
13const P: usize = G;
14
15#[cfg(test)]
16const Q: usize = N;
17
18pub type Triangle<'a, C> = SurfaceElement<'a, C, G, N, P>;
19
20impl<'a, C> SurfaceFiniteElement<'a, C, G, N, P> for Triangle<'a, C> {
21 fn new(
22 constitutive_model: &'a C,
23 reference_nodal_coordinates: ReferenceNodalCoordinates<N>,
24 thickness: &Scalar,
25 ) -> Self {
26 let integration_weights = Self::bases(&reference_nodal_coordinates)
27 .iter()
28 .map(|reference_basis| {
29 reference_basis[0].cross(&reference_basis[1]).norm()
30 * Self::integration_weight()
31 * thickness
32 })
33 .collect();
34 let reference_dual_bases = Self::dual_bases(&reference_nodal_coordinates);
35 let gradient_vectors = Self::standard_gradient_operators()
36 .iter()
37 .zip(reference_dual_bases.iter())
38 .map(|(standard_gradient_operator, reference_dual_basis)| {
39 standard_gradient_operator
40 .iter()
41 .map(|standard_gradient_operator_a| {
42 standard_gradient_operator_a
43 .iter()
44 .zip(reference_dual_basis.iter())
45 .map(|(standard_gradient_operator_a_m, reference_dual_basis_m)| {
46 reference_dual_basis_m * standard_gradient_operator_a_m
47 })
48 .sum()
49 })
50 .collect()
51 })
52 .collect();
53 let reference_normals = reference_dual_bases
54 .iter()
55 .map(|reference_dual_basis| {
56 reference_dual_basis[0]
57 .cross(&reference_dual_basis[1])
58 .normalized()
59 })
60 .collect();
61 Self {
62 constitutive_model,
63 gradient_vectors,
64 integration_weights,
65 reference_normals,
66 }
67 }
68}
69
70impl<'a, C> Triangle<'a, C> {
71 const fn integration_weight() -> Scalar {
72 1.0 / 2.0
73 }
74 #[cfg(test)]
75 const fn shape_functions_at_integration_points() -> ShapeFunctionsAtIntegrationPoints<G, Q> {
76 ShapeFunctionsAtIntegrationPoints::const_from([tensor_rank_1([1.0 / 3.0; Q])])
77 }
78 const fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
79 StandardGradientOperators::const_from([TensorRank1List::const_from([
80 tensor_rank_1([-1.0, -1.0]),
81 tensor_rank_1([1.0, 0.0]),
82 tensor_rank_1([0.0, 1.0]),
83 ])])
84 }
85}
86
87impl<'a, C> SurfaceFiniteElementMethodsExtra<M, N, P> for Triangle<'a, C> {
88 fn standard_gradient_operators() -> StandardGradientOperators<M, N, P> {
89 Self::standard_gradient_operators()
90 }
91}
92
93impl<'a, C> FiniteElementMethods<G, N> for Triangle<'a, C> {
94 fn deformation_gradients(
95 &self,
96 nodal_coordinates: &NodalCoordinates<N>,
97 ) -> DeformationGradientList<G> {
98 self.gradient_vectors()
99 .iter()
100 .zip(
101 Self::normals(nodal_coordinates)
102 .iter()
103 .zip(self.reference_normals().iter()),
104 )
105 .map(|(gradient_vectors, (normal, reference_normal))| {
106 nodal_coordinates
107 .iter()
108 .zip(gradient_vectors.iter())
109 .map(|(nodal_coordinate, gradient_vector)| {
110 DeformationGradient::dyad(nodal_coordinate, gradient_vector)
111 })
112 .sum::<DeformationGradient>()
113 + DeformationGradient::dyad(normal, reference_normal)
114 })
115 .collect()
116 }
117 fn deformation_gradient_rates(
118 &self,
119 nodal_coordinates: &NodalCoordinates<N>,
120 nodal_velocities: &NodalVelocities<N>,
121 ) -> DeformationGradientRateList<G> {
122 self.gradient_vectors()
123 .iter()
124 .zip(
125 Self::normal_rates(nodal_coordinates, nodal_velocities)
126 .iter()
127 .zip(self.reference_normals().iter()),
128 )
129 .map(|(gradient_vectors, (normal_rate, reference_normal))| {
130 nodal_velocities
131 .iter()
132 .zip(gradient_vectors.iter())
133 .map(|(nodal_velocity, gradient_vector)| {
134 DeformationGradientRate::dyad(nodal_velocity, gradient_vector)
135 })
136 .sum::<DeformationGradientRate>()
137 + DeformationGradientRate::dyad(normal_rate, reference_normal)
138 })
139 .collect()
140 }
141 fn gradient_vectors(&self) -> &GradientVectors<G, N> {
142 &self.gradient_vectors
143 }
144 fn integration_weights(&self) -> &Scalars<G> {
145 &self.integration_weights
146 }
147}
148
149impl<'a, C> ElasticFiniteElement<C, G, N> for Triangle<'a, C>
150where
151 C: Elastic,
152{
153 fn nodal_forces(
154 &self,
155 nodal_coordinates: &NodalCoordinates<N>,
156 ) -> Result<NodalForces<N>, FiniteElementError> {
157 match self
158 .deformation_gradients(nodal_coordinates)
159 .iter()
160 .map(|deformation_gradient| {
161 self.constitutive_model
162 .first_piola_kirchhoff_stress(deformation_gradient)
163 })
164 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
165 {
166 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
167 .iter()
168 .zip(
169 self.gradient_vectors()
170 .iter()
171 .zip(self.integration_weights().iter()),
172 )
173 .map(
174 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
175 gradient_vectors
176 .iter()
177 .map(|gradient_vector| {
178 (first_piola_kirchhoff_stress * gradient_vector)
179 * integration_weight
180 })
181 .collect()
182 },
183 )
184 .sum()),
185 Err(error) => Err(FiniteElementError::Upstream(
186 format!("{error}"),
187 format!("{self:?}"),
188 )),
189 }
190 }
191 fn nodal_stiffnesses(
192 &self,
193 nodal_coordinates: &NodalCoordinates<N>,
194 ) -> Result<NodalStiffnesses<N>, FiniteElementError> {
195 match self.deformation_gradients(nodal_coordinates).iter()
196 .map(|deformation_gradient| {
197 self.constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
198 })
199 .collect::<Result<FirstPiolaKirchhoffTangentStiffnesses<G>, _>>()
200 {
201 Ok(first_piola_kirchhoff_tangent_stiffnesses) => Ok(first_piola_kirchhoff_tangent_stiffnesses
202 .iter()
203 .zip(
204 self.gradient_vectors()
205 .iter()
206 .zip(self.integration_weights().iter()
207 .zip(self.reference_normals().iter()
208 .zip(Self::normal_gradients(nodal_coordinates).iter())
209 )
210 ),
211 )
212 .map(
213 |(
214 first_piola_kirchhoff_tangent_stiffness,
215 (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
216 )| {
217 gradient_vectors.iter()
218 .map(|gradient_vector_a|
219 gradient_vectors.iter()
220 .zip(normal_gradients.iter())
221 .map(|(gradient_vector_b, normal_gradient_b)|
222 first_piola_kirchhoff_tangent_stiffness.iter()
223 .map(|first_piola_kirchhoff_tangent_stiffness_m|
224 IDENTITY.iter()
225 .zip(normal_gradient_b.iter())
226 .map(|(identity_n, normal_gradient_b_n)|
227 first_piola_kirchhoff_tangent_stiffness_m.iter()
228 .zip(gradient_vector_a.iter())
229 .map(|(first_piola_kirchhoff_tangent_stiffness_mj, gradient_vector_a_j)|
230 first_piola_kirchhoff_tangent_stiffness_mj.iter()
231 .zip(identity_n.iter()
232 .zip(normal_gradient_b_n.iter()))
233 .map(|(first_piola_kirchhoff_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
234 first_piola_kirchhoff_tangent_stiffness_mjk.iter()
235 .zip(gradient_vector_b.iter()
236 .zip(reference_normal.iter()))
237 .map(|(first_piola_kirchhoff_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
238 first_piola_kirchhoff_tangent_stiffness_mjkl * gradient_vector_a_j * (
239 identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
240 ) * integration_weight
241 ).sum::<Scalar>()
242 ).sum::<Scalar>()
243 ).sum::<Scalar>()
244 ).collect()
245 ).collect()
246 ).collect()
247 ).collect()
248 }
249 )
250 .sum()),
251 Err(error) => Err(FiniteElementError::Upstream(
252 format!("{error}"),
253 format!("{self:?}"),
254 )),
255 }
256 }
257}
258
259impl<'a, C> HyperelasticFiniteElement<C, G, N> for Triangle<'a, C>
260where
261 C: Hyperelastic,
262{
263 fn helmholtz_free_energy(
264 &self,
265 nodal_coordinates: &NodalCoordinates<N>,
266 ) -> Result<Scalar, FiniteElementError> {
267 match self
268 .deformation_gradients(nodal_coordinates)
269 .iter()
270 .zip(self.integration_weights().iter())
271 .map(|(deformation_gradient, integration_weight)| {
272 Ok::<Scalar, ConstitutiveError>(
273 self.constitutive_model
274 .helmholtz_free_energy_density(deformation_gradient)?
275 * integration_weight,
276 )
277 })
278 .sum()
279 {
280 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
281 Err(error) => Err(FiniteElementError::Upstream(
282 format!("{error}"),
283 format!("{self:?}"),
284 )),
285 }
286 }
287}
288
289impl<'a, C> ViscoelasticFiniteElement<C, G, N> for Triangle<'a, C>
290where
291 C: Viscoelastic,
292{
293 fn nodal_forces(
294 &self,
295 nodal_coordinates: &NodalCoordinates<N>,
296 nodal_velocities: &NodalVelocities<N>,
297 ) -> Result<NodalForces<N>, FiniteElementError> {
298 match self
299 .deformation_gradients(nodal_coordinates)
300 .iter()
301 .zip(
302 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
303 .iter(),
304 )
305 .map(|(deformation_gradient, deformation_gradient_rate)| {
306 self.constitutive_model
307 .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
308 })
309 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
310 {
311 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
312 .iter()
313 .zip(
314 self.gradient_vectors()
315 .iter()
316 .zip(self.integration_weights().iter()),
317 )
318 .map(
319 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
320 gradient_vectors
321 .iter()
322 .map(|gradient_vector| {
323 (first_piola_kirchhoff_stress * gradient_vector)
324 * integration_weight
325 })
326 .collect()
327 },
328 )
329 .sum()),
330 Err(error) => Err(FiniteElementError::Upstream(
331 format!("{error}"),
332 format!("{self:?}"),
333 )),
334 }
335 }
336 fn nodal_stiffnesses(
337 &self,
338 nodal_coordinates: &NodalCoordinates<N>,
339 nodal_velocities: &NodalVelocities<N>,
340 ) -> Result<NodalStiffnesses<N>, FiniteElementError> {
341 match self.deformation_gradients(nodal_coordinates).iter().zip(self.deformation_gradient_rates(nodal_coordinates, nodal_velocities).iter())
342 .map(|(deformation_gradient, deformation_gradient_rate)| {
343 self.constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)
344 })
345 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
346 {
347 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
348 .iter()
349 .zip(
350 self.gradient_vectors()
351 .iter()
352 .zip(self.integration_weights().iter()
353 .zip(self.reference_normals().iter()
354 .zip(Self::normal_gradients(nodal_coordinates).iter())
355 )
356 ),
357 )
358 .map(
359 |(
360 first_piola_kirchoff_rate_tangent_stiffness_mjkl,
361 (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
362 )| {
363 gradient_vectors.iter()
364 .map(|gradient_vector_a|
365 gradient_vectors.iter()
366 .zip(normal_gradients.iter())
367 .map(|(gradient_vector_b, normal_gradient_b)|
368 first_piola_kirchoff_rate_tangent_stiffness_mjkl.iter()
369 .map(|first_piola_kirchhoff_rate_tangent_stiffness_m|
370 IDENTITY.iter()
371 .zip(normal_gradient_b.iter())
372 .map(|(identity_n, normal_gradient_b_n)|
373 first_piola_kirchhoff_rate_tangent_stiffness_m.iter()
374 .zip(gradient_vector_a.iter())
375 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mj, gradient_vector_a_j)|
376 first_piola_kirchhoff_rate_tangent_stiffness_mj.iter()
377 .zip(identity_n.iter()
378 .zip(normal_gradient_b_n.iter()))
379 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
380 first_piola_kirchhoff_rate_tangent_stiffness_mjk.iter()
381 .zip(gradient_vector_b.iter()
382 .zip(reference_normal.iter()))
383 .map(|(first_piola_kirchoff_rate_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
384 first_piola_kirchoff_rate_tangent_stiffness_mjkl * gradient_vector_a_j * (
385 identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
386 ) * integration_weight
387 ).sum::<Scalar>()
388 ).sum::<Scalar>()
389 ).sum::<Scalar>()
390 ).collect()
391 ).collect()
392 ).collect()
393 ).collect()
394 }
395 )
396 .sum()),
397 Err(error) => Err(FiniteElementError::Upstream(
398 format!("{error}"),
399 format!("{self:?}"),
400 )),
401 }
402 }
403}
404
405impl<'a, C> ElasticHyperviscousFiniteElement<C, G, N> for Triangle<'a, C>
406where
407 C: ElasticHyperviscous,
408{
409 fn viscous_dissipation(
410 &self,
411 nodal_coordinates: &NodalCoordinates<N>,
412 nodal_velocities: &NodalVelocities<N>,
413 ) -> Result<Scalar, FiniteElementError> {
414 match self
415 .deformation_gradients(nodal_coordinates)
416 .iter()
417 .zip(
418 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
419 .iter()
420 .zip(self.integration_weights().iter()),
421 )
422 .map(
423 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
424 Ok::<Scalar, ConstitutiveError>(
425 self.constitutive_model
426 .viscous_dissipation(deformation_gradient, deformation_gradient_rate)?
427 * integration_weight,
428 )
429 },
430 )
431 .sum()
432 {
433 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
434 Err(error) => Err(FiniteElementError::Upstream(
435 format!("{error}"),
436 format!("{self:?}"),
437 )),
438 }
439 }
440 fn dissipation_potential(
441 &self,
442 nodal_coordinates: &NodalCoordinates<N>,
443 nodal_velocities: &NodalVelocities<N>,
444 ) -> Result<Scalar, FiniteElementError> {
445 match self
446 .deformation_gradients(nodal_coordinates)
447 .iter()
448 .zip(
449 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
450 .iter()
451 .zip(self.integration_weights().iter()),
452 )
453 .map(
454 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
455 Ok::<Scalar, ConstitutiveError>(
456 self.constitutive_model.dissipation_potential(
457 deformation_gradient,
458 deformation_gradient_rate,
459 )? * integration_weight,
460 )
461 },
462 )
463 .sum()
464 {
465 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
466 Err(error) => Err(FiniteElementError::Upstream(
467 format!("{error}"),
468 format!("{self:?}"),
469 )),
470 }
471 }
472}
473
474impl<'a, C> HyperviscoelasticFiniteElement<C, G, N> for Triangle<'a, C>
475where
476 C: Hyperviscoelastic,
477{
478 fn helmholtz_free_energy(
479 &self,
480 nodal_coordinates: &NodalCoordinates<N>,
481 ) -> Result<Scalar, FiniteElementError> {
482 match self
483 .deformation_gradients(nodal_coordinates)
484 .iter()
485 .zip(self.integration_weights().iter())
486 .map(|(deformation_gradient, integration_weight)| {
487 Ok::<Scalar, ConstitutiveError>(
488 self.constitutive_model
489 .helmholtz_free_energy_density(deformation_gradient)?
490 * integration_weight,
491 )
492 })
493 .sum()
494 {
495 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
496 Err(error) => Err(FiniteElementError::Upstream(
497 format!("{error}"),
498 format!("{self:?}"),
499 )),
500 }
501 }
502}