1#[cfg(test)]
2mod test;
3
4pub mod composite;
5pub mod linear;
6
7use super::*;
8use crate::{
9 defeat_message,
10 math::{IDENTITY, LEVI_CIVITA, tensor_rank_1_zero},
11 mechanics::Scalar,
12};
13use std::{
14 any::type_name,
15 fmt::{Debug, Display},
16};
17
18pub struct Element<'a, C, const G: usize, const N: usize> {
19 constitutive_model: &'a C,
20 gradient_vectors: GradientVectors<G, N>,
21 integration_weights: Scalars<G>,
22}
23
24#[derive(Debug)]
25pub struct SurfaceElement<'a, C, const G: usize, const N: usize, const P: usize> {
26 constitutive_model: &'a C,
27 gradient_vectors: GradientVectors<G, N>,
28 integration_weights: Scalars<G>,
29 reference_normals: ReferenceNormals<P>,
30}
31
32impl<'a, C, const G: usize, const N: usize> Debug for Element<'a, C, G, N> {
33 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
34 let element = match (G, N) {
35 (1, 3) => "LinearTriangle",
36 (1, 4) => "LinearTetrahedron",
37 (8, 8) => "LinearHexahedron",
38 (4, 10) => "CompositeTetrahedron",
39 _ => panic!(),
40 };
41 write!(
42 f,
43 "{element} {{ constitutive_model: {} }}",
44 type_name::<C>()
45 .rsplit("::")
46 .next()
47 .unwrap()
48 .split("<")
49 .next()
50 .unwrap()
51 )
52 }
53}
54
55pub trait FiniteElement<'a, C, const G: usize, const N: usize>
56where
57 Self: FiniteElementMethods<G, N>,
58{
59 fn new(
60 constitutive_model: &'a C,
61 reference_nodal_coordinates: ReferenceNodalCoordinates<N>,
62 ) -> Self;
63 fn reference() -> ReferenceNodalCoordinates<N>;
64 fn reset(&mut self);
65}
66
67pub trait SurfaceFiniteElement<'a, C, const G: usize, const N: usize, const P: usize>
68where
69 Self: FiniteElementMethods<G, N>,
70{
71 fn new(
72 constitutive_model: &'a C,
73 reference_nodal_coordinates: ReferenceNodalCoordinates<N>,
74 thickness: &Scalar,
75 ) -> Self;
76}
77
78pub trait FiniteElementMethods<const G: usize, const N: usize> {
79 fn deformation_gradients(
80 &self,
81 nodal_coordinates: &NodalCoordinates<N>,
82 ) -> DeformationGradientList<G>;
83 fn deformation_gradient_rates(
84 &self,
85 nodal_coordinates: &NodalCoordinates<N>,
86 nodal_velocities: &NodalVelocities<N>,
87 ) -> DeformationGradientRateList<G>;
88 fn gradient_vectors(&self) -> &GradientVectors<G, N>;
89 fn integration_weights(&self) -> &Scalars<G>;
90}
91
92pub trait SurfaceFiniteElementMethods<
93 const G: usize,
94 const M: usize,
95 const N: usize,
96 const P: usize,
97> where
98 Self: SurfaceFiniteElementMethodsExtra<M, N, P>,
99{
100 fn bases<const I: usize>(nodal_coordinates: &Coordinates<I, N>) -> Bases<I, P>;
101 fn dual_bases<const I: usize>(nodal_coordinates: &Coordinates<I, N>) -> Bases<I, P>;
102 fn normals(nodal_coordinates: &NodalCoordinates<N>) -> Normals<P>;
103 fn normal_gradients(nodal_coordinates: &NodalCoordinates<N>) -> NormalGradients<N, P>;
104 fn normal_rates(
105 nodal_coordinates: &NodalCoordinates<N>,
106 nodal_velocities: &NodalVelocities<N>,
107 ) -> NormalRates<P>;
108 fn reference_normals(&self) -> &ReferenceNormals<P>;
109}
110
111pub trait SurfaceFiniteElementMethodsExtra<const M: usize, const N: usize, const P: usize> {
113 fn standard_gradient_operators() -> StandardGradientOperators<M, N, P>;
114}
115
116pub enum FiniteElementError {
117 Upstream(String, String),
118}
119
120impl From<FiniteElementError> for TestError {
121 fn from(error: FiniteElementError) -> Self {
122 Self {
123 message: error.to_string(),
124 }
125 }
126}
127
128impl Debug for FiniteElementError {
129 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
130 let error = match self {
131 Self::Upstream(error, element) => {
132 format!(
133 "{error}\x1b[0;91m\n\
134 In finite element: {element}."
135 )
136 }
137 };
138 write!(f, "\n{error}\n\x1b[0;2;31m{}\x1b[0m\n", defeat_message())
139 }
140}
141
142impl Display for FiniteElementError {
143 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
144 let error = match self {
145 Self::Upstream(error, element) => {
146 format!(
147 "{error}\x1b[0;91m\n\
148 In finite element: {element}."
149 )
150 }
151 };
152 write!(f, "{error}\x1b[0m")
153 }
154}
155
156impl<'a, C, const G: usize, const N: usize> FiniteElementMethods<G, N> for Element<'a, C, G, N> {
157 fn deformation_gradients(
158 &self,
159 nodal_coordinates: &NodalCoordinates<N>,
160 ) -> DeformationGradientList<G> {
161 self.gradient_vectors()
162 .iter()
163 .map(|gradient_vectors| {
164 nodal_coordinates
165 .iter()
166 .zip(gradient_vectors.iter())
167 .map(|(nodal_coordinate, gradient_vector)| {
168 DeformationGradient::dyad(nodal_coordinate, gradient_vector)
169 })
170 .sum()
171 })
172 .collect()
173 }
174 fn deformation_gradient_rates(
175 &self,
176 _: &NodalCoordinates<N>,
177 nodal_velocities: &NodalVelocities<N>,
178 ) -> DeformationGradientRateList<G> {
179 self.gradient_vectors()
180 .iter()
181 .map(|gradient_vectors| {
182 nodal_velocities
183 .iter()
184 .zip(gradient_vectors.iter())
185 .map(|(nodal_velocity, gradient_vector)| {
186 DeformationGradientRate::dyad(nodal_velocity, gradient_vector)
187 })
188 .sum()
189 })
190 .collect()
191 }
192 fn gradient_vectors(&self) -> &GradientVectors<G, N> {
193 &self.gradient_vectors
194 }
195 fn integration_weights(&self) -> &Scalars<G> {
196 &self.integration_weights
197 }
198}
199
200impl<'a, C, const G: usize, const M: usize, const N: usize, const P: usize>
201 SurfaceFiniteElementMethods<G, M, N, P> for SurfaceElement<'a, C, G, N, P>
202where
203 Self: SurfaceFiniteElementMethodsExtra<M, N, P>,
204{
205 fn bases<const I: usize>(nodal_coordinates: &Coordinates<I, N>) -> Bases<I, P> {
206 Self::standard_gradient_operators()
207 .iter()
208 .map(|standard_gradient_operator| {
209 standard_gradient_operator
210 .iter()
211 .zip(nodal_coordinates.iter())
212 .map(|(standard_gradient_operator_a, nodal_coordinate_a)| {
213 standard_gradient_operator_a
214 .iter()
215 .map(|standard_gradient_operator_a_m| {
216 nodal_coordinate_a * standard_gradient_operator_a_m
217 })
218 .collect()
219 })
220 .sum()
221 })
222 .collect()
223 }
224 fn dual_bases<const I: usize>(nodal_coordinates: &Coordinates<I, N>) -> Bases<I, P> {
225 Self::bases(nodal_coordinates)
226 .iter()
227 .map(|basis_vectors| {
228 basis_vectors
229 .iter()
230 .map(|basis_vectors_m| {
231 basis_vectors
232 .iter()
233 .map(|basis_vectors_n| basis_vectors_m * basis_vectors_n)
234 .collect()
235 })
236 .collect::<TensorRank2<2, I, I>>()
237 .inverse()
238 .iter()
239 .map(|metric_tensor_m| {
240 metric_tensor_m
241 .iter()
242 .zip(basis_vectors.iter())
243 .map(|(metric_tensor_mn, basis_vectors_n)| {
244 basis_vectors_n * metric_tensor_mn
245 })
246 .sum()
247 })
248 .collect()
249 })
250 .collect()
251 }
252 fn normals(nodal_coordinates: &NodalCoordinates<N>) -> Normals<P> {
253 Self::bases(nodal_coordinates)
254 .iter()
255 .map(|basis_vectors| basis_vectors[0].cross(&basis_vectors[1]).normalized())
256 .collect()
257 }
258 fn normal_gradients(nodal_coordinates: &NodalCoordinates<N>) -> NormalGradients<N, P> {
259 let levi_civita_symbol = LEVI_CIVITA;
260 let mut normalization: Scalar = 0.0;
261 let mut normal_vector = tensor_rank_1_zero();
262 Self::standard_gradient_operators().iter()
263 .zip(Self::bases(nodal_coordinates).iter())
264 .map(|(standard_gradient_operator, basis_vectors)|{
265 normalization = basis_vectors[0].cross(&basis_vectors[1]).norm();
266 normal_vector = basis_vectors[0].cross(&basis_vectors[1])/normalization;
267 standard_gradient_operator.iter()
268 .map(|standard_gradient_operator_a|
269 levi_civita_symbol.iter()
270 .map(|levi_civita_symbol_m|
271 IDENTITY.iter()
272 .zip(normal_vector.iter())
273 .map(|(identity_i, normal_vector_i)|
274 levi_civita_symbol_m.iter()
275 .zip(basis_vectors[0].iter()
276 .zip(basis_vectors[1].iter()))
277 .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
278 levi_civita_symbol_mn.iter()
279 .zip(identity_i.iter()
280 .zip(normal_vector.iter()))
281 .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
282 levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
283 ).sum::<Scalar>() * (
284 standard_gradient_operator_a[0] * basis_vector_1_n
285 - standard_gradient_operator_a[1] * basis_vector_0_n
286 )
287 ).sum::<Scalar>() / normalization
288 ).collect()
289 ).collect()
290 ).collect()
291 }).collect()
292 }
293 fn normal_rates(
294 nodal_coordinates: &NodalCoordinates<N>,
295 nodal_velocities: &NodalVelocities<N>,
296 ) -> NormalRates<P> {
297 let identity = IDENTITY;
298 let levi_civita_symbol = LEVI_CIVITA;
299 let mut normalization = 0.0;
300 Self::bases(nodal_coordinates)
301 .iter()
302 .zip(Self::normals(nodal_coordinates).iter()
303 .zip(Self::standard_gradient_operators().iter()))
304 .map(|(basis, (normal, standard_gradient_operator))| {
305 normalization = basis[0].cross(&basis[1]).norm();
306 identity.iter()
307 .zip(normal.iter())
308 .map(|(identity_i, normal_vector_i)|
309 nodal_velocities.iter()
310 .zip(standard_gradient_operator.iter())
311 .map(|(nodal_velocity_a, standard_gradient_operator_a)|
312 levi_civita_symbol.iter()
313 .zip(nodal_velocity_a.iter())
314 .map(|(levi_civita_symbol_m, nodal_velocity_a_m)|
315 levi_civita_symbol_m.iter()
316 .zip(basis[0].iter()
317 .zip(basis[1].iter()))
318 .map(|(levi_civita_symbol_mn, (basis_vector_0_n, basis_vector_1_n))|
319 levi_civita_symbol_mn.iter()
320 .zip(identity_i.iter()
321 .zip(normal.iter()))
322 .map(|(levi_civita_symbol_mno, (identity_io, normal_vector_o))|
323 levi_civita_symbol_mno * (identity_io - normal_vector_i * normal_vector_o)
324 ).sum::<Scalar>() * (
325 standard_gradient_operator_a[0] * basis_vector_1_n
326 - standard_gradient_operator_a[1] * basis_vector_0_n
327 )
328 ).sum::<Scalar>() * nodal_velocity_a_m
329 ).sum::<Scalar>()
330 ).sum::<Scalar>() / normalization
331 ).collect()
332 }).collect()
333 }
334 fn reference_normals(&self) -> &ReferenceNormals<P> {
335 &self.reference_normals
336 }
337}
338
339pub trait ElasticFiniteElement<C, const G: usize, const N: usize>
340where
341 C: Elastic,
342 Self: Debug + FiniteElementMethods<G, N>,
343{
344 fn nodal_forces(
345 &self,
346 nodal_coordinates: &NodalCoordinates<N>,
347 ) -> Result<NodalForces<N>, FiniteElementError>;
348 fn nodal_stiffnesses(
349 &self,
350 nodal_coordinates: &NodalCoordinates<N>,
351 ) -> Result<NodalStiffnesses<N>, FiniteElementError>;
352}
353
354pub trait HyperelasticFiniteElement<C, const G: usize, const N: usize>
355where
356 C: Hyperelastic,
357 Self: ElasticFiniteElement<C, G, N>,
358{
359 fn helmholtz_free_energy(
360 &self,
361 nodal_coordinates: &NodalCoordinates<N>,
362 ) -> Result<Scalar, FiniteElementError>;
363}
364
365pub trait ViscoelasticFiniteElement<C, const G: usize, const N: usize>
366where
367 C: Viscoelastic,
368 Self: FiniteElementMethods<G, N>,
369{
370 fn nodal_forces(
371 &self,
372 nodal_coordinates: &NodalCoordinates<N>,
373 nodal_velocities: &NodalVelocities<N>,
374 ) -> Result<NodalForces<N>, FiniteElementError>;
375 fn nodal_stiffnesses(
376 &self,
377 nodal_coordinates: &NodalCoordinates<N>,
378 nodal_velocities: &NodalVelocities<N>,
379 ) -> Result<NodalStiffnesses<N>, FiniteElementError>;
380}
381
382pub trait ElasticHyperviscousFiniteElement<C, const G: usize, const N: usize>
383where
384 C: ElasticHyperviscous,
385 Self: ViscoelasticFiniteElement<C, G, N>,
386{
387 fn viscous_dissipation(
388 &self,
389 nodal_coordinates: &NodalCoordinates<N>,
390 nodal_velocities: &NodalVelocities<N>,
391 ) -> Result<Scalar, FiniteElementError>;
392 fn dissipation_potential(
393 &self,
394 nodal_coordinates: &NodalCoordinates<N>,
395 nodal_velocities: &NodalVelocities<N>,
396 ) -> Result<Scalar, FiniteElementError>;
397}
398
399pub trait HyperviscoelasticFiniteElement<C, const G: usize, const N: usize>
400where
401 C: Hyperviscoelastic,
402 Self: ElasticHyperviscousFiniteElement<C, G, N>,
403{
404 fn helmholtz_free_energy(
405 &self,
406 nodal_coordinates: &NodalCoordinates<N>,
407 ) -> Result<Scalar, FiniteElementError>;
408}
409
410impl<'a, C, const G: usize, const N: usize> ElasticFiniteElement<C, G, N> for Element<'a, C, G, N>
411where
412 C: Elastic,
413{
414 fn nodal_forces(
415 &self,
416 nodal_coordinates: &NodalCoordinates<N>,
417 ) -> Result<NodalForces<N>, FiniteElementError> {
418 match self
419 .deformation_gradients(nodal_coordinates)
420 .iter()
421 .map(|deformation_gradient| {
422 self.constitutive_model
423 .first_piola_kirchhoff_stress(deformation_gradient)
424 })
425 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
426 {
427 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
428 .iter()
429 .zip(
430 self.gradient_vectors()
431 .iter()
432 .zip(self.integration_weights().iter()),
433 )
434 .map(
435 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
436 gradient_vectors
437 .iter()
438 .map(|gradient_vector| {
439 (first_piola_kirchhoff_stress * gradient_vector)
440 * integration_weight
441 })
442 .collect()
443 },
444 )
445 .sum()),
446 Err(error) => Err(FiniteElementError::Upstream(
447 format!("{error}"),
448 format!("{self:?}"),
449 )),
450 }
451 }
452 fn nodal_stiffnesses(
453 &self,
454 nodal_coordinates: &NodalCoordinates<N>,
455 ) -> Result<NodalStiffnesses<N>, FiniteElementError> {
456 match self
457 .deformation_gradients(nodal_coordinates)
458 .iter()
459 .map(|deformation_gradient| {
460 self.constitutive_model
461 .first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
462 })
463 .collect::<Result<FirstPiolaKirchhoffTangentStiffnesses<G>, _>>()
464 {
465 Ok(first_piola_kirchhoff_tangent_stiffnesses) => {
466 Ok(first_piola_kirchhoff_tangent_stiffnesses
467 .iter()
468 .zip(
469 self.gradient_vectors()
470 .iter()
471 .zip(self.integration_weights().iter()),
472 )
473 .map(
474 |(
475 first_piola_kirchhoff_tangent_stiffness,
476 (gradient_vectors, integration_weight),
477 )| {
478 gradient_vectors
479 .iter()
480 .map(|gradient_vector_a| {
481 gradient_vectors
482 .iter()
483 .map(|gradient_vector_b| {
484 first_piola_kirchhoff_tangent_stiffness
485 .contract_second_fourth_indices_with_first_indices_of(
486 gradient_vector_a,
487 gradient_vector_b,
488 )
489 * integration_weight
490 })
491 .collect()
492 })
493 .collect()
494 },
495 )
496 .sum())
497 }
498 Err(error) => Err(FiniteElementError::Upstream(
499 format!("{error}"),
500 format!("{self:?}"),
501 )),
502 }
503 }
504}
505
506impl<'a, C, const G: usize, const N: usize> HyperelasticFiniteElement<C, G, N>
507 for Element<'a, C, G, N>
508where
509 C: Hyperelastic,
510{
511 fn helmholtz_free_energy(
512 &self,
513 nodal_coordinates: &NodalCoordinates<N>,
514 ) -> Result<Scalar, FiniteElementError> {
515 match self
516 .deformation_gradients(nodal_coordinates)
517 .iter()
518 .zip(self.integration_weights().iter())
519 .map(|(deformation_gradient, integration_weight)| {
520 Ok::<Scalar, ConstitutiveError>(
521 self.constitutive_model
522 .helmholtz_free_energy_density(deformation_gradient)?
523 * integration_weight,
524 )
525 })
526 .sum()
527 {
528 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
529 Err(error) => Err(FiniteElementError::Upstream(
530 format!("{error}"),
531 format!("{self:?}"),
532 )),
533 }
534 }
535}
536
537impl<'a, C, const G: usize, const N: usize> ViscoelasticFiniteElement<C, G, N>
538 for Element<'a, C, G, N>
539where
540 C: Viscoelastic,
541{
542 fn nodal_forces(
543 &self,
544 nodal_coordinates: &NodalCoordinates<N>,
545 nodal_velocities: &NodalVelocities<N>,
546 ) -> Result<NodalForces<N>, FiniteElementError> {
547 match self
548 .deformation_gradients(nodal_coordinates)
549 .iter()
550 .zip(
551 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
552 .iter(),
553 )
554 .map(|(deformation_gradient, deformation_gradient_rate)| {
555 self.constitutive_model
556 .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
557 })
558 .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
559 {
560 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
561 .iter()
562 .zip(
563 self.gradient_vectors()
564 .iter()
565 .zip(self.integration_weights().iter()),
566 )
567 .map(
568 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
569 gradient_vectors
570 .iter()
571 .map(|gradient_vector| {
572 (first_piola_kirchhoff_stress * gradient_vector)
573 * integration_weight
574 })
575 .collect()
576 },
577 )
578 .sum()),
579 Err(error) => Err(FiniteElementError::Upstream(
580 format!("{error}"),
581 format!("{self:?}"),
582 )),
583 }
584 }
585 fn nodal_stiffnesses(
586 &self,
587 nodal_coordinates: &NodalCoordinates<N>,
588 nodal_velocities: &NodalVelocities<N>,
589 ) -> Result<NodalStiffnesses<N>, FiniteElementError> {
590 match self
591 .deformation_gradients(nodal_coordinates)
592 .iter()
593 .zip(
594 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
595 .iter(),
596 )
597 .map(|(deformation_gradient, deformation_gradient_rate)| {
598 self.constitutive_model
599 .first_piola_kirchhoff_rate_tangent_stiffness(
600 deformation_gradient,
601 deformation_gradient_rate,
602 )
603 })
604 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
605 {
606 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => {
607 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
608 .iter()
609 .zip(
610 self.gradient_vectors().iter().zip(
611 self.gradient_vectors()
612 .iter()
613 .zip(self.integration_weights().iter()),
614 ),
615 )
616 .map(
617 |(
618 first_piola_kirchhoff_rate_tangent_stiffness,
619 (gradient_vectors_a, (gradient_vectors_b, integration_weight)),
620 )| {
621 gradient_vectors_a
622 .iter()
623 .map(|gradient_vector_a| {
624 gradient_vectors_b
625 .iter()
626 .map(|gradient_vector_b| {
627 first_piola_kirchhoff_rate_tangent_stiffness
628 .contract_second_fourth_indices_with_first_indices_of(
629 gradient_vector_a,
630 gradient_vector_b,
631 )
632 * integration_weight
633 })
634 .collect()
635 })
636 .collect()
637 },
638 )
639 .sum())
640 }
641 Err(error) => Err(FiniteElementError::Upstream(
642 format!("{error}"),
643 format!("{self:?}"),
644 )),
645 }
646 }
647}
648
649impl<'a, C, const G: usize, const N: usize> ElasticHyperviscousFiniteElement<C, G, N>
650 for Element<'a, C, G, N>
651where
652 C: ElasticHyperviscous,
653{
654 fn viscous_dissipation(
655 &self,
656 nodal_coordinates: &NodalCoordinates<N>,
657 nodal_velocities: &NodalVelocities<N>,
658 ) -> Result<Scalar, FiniteElementError> {
659 match self
660 .deformation_gradients(nodal_coordinates)
661 .iter()
662 .zip(
663 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
664 .iter()
665 .zip(self.integration_weights().iter()),
666 )
667 .map(
668 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
669 Ok::<Scalar, ConstitutiveError>(
670 self.constitutive_model
671 .viscous_dissipation(deformation_gradient, deformation_gradient_rate)?
672 * integration_weight,
673 )
674 },
675 )
676 .sum()
677 {
678 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
679 Err(error) => Err(FiniteElementError::Upstream(
680 format!("{error}"),
681 format!("{self:?}"),
682 )),
683 }
684 }
685 fn dissipation_potential(
686 &self,
687 nodal_coordinates: &NodalCoordinates<N>,
688 nodal_velocities: &NodalVelocities<N>,
689 ) -> Result<Scalar, FiniteElementError> {
690 match self
691 .deformation_gradients(nodal_coordinates)
692 .iter()
693 .zip(
694 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
695 .iter()
696 .zip(self.integration_weights().iter()),
697 )
698 .map(
699 |(deformation_gradient, (deformation_gradient_rate, integration_weight))| {
700 Ok::<Scalar, ConstitutiveError>(
701 self.constitutive_model.dissipation_potential(
702 deformation_gradient,
703 deformation_gradient_rate,
704 )? * integration_weight,
705 )
706 },
707 )
708 .sum()
709 {
710 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
711 Err(error) => Err(FiniteElementError::Upstream(
712 format!("{error}"),
713 format!("{self:?}"),
714 )),
715 }
716 }
717}
718
719impl<'a, C, const G: usize, const N: usize> HyperviscoelasticFiniteElement<C, G, N>
720 for Element<'a, C, G, N>
721where
722 C: Hyperviscoelastic,
723{
724 fn helmholtz_free_energy(
725 &self,
726 nodal_coordinates: &NodalCoordinates<N>,
727 ) -> Result<Scalar, FiniteElementError> {
728 match self
729 .deformation_gradients(nodal_coordinates)
730 .iter()
731 .zip(self.integration_weights().iter())
732 .map(|(deformation_gradient, integration_weight)| {
733 Ok::<Scalar, ConstitutiveError>(
734 self.constitutive_model
735 .helmholtz_free_energy_density(deformation_gradient)?
736 * integration_weight,
737 )
738 })
739 .sum()
740 {
741 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
742 Err(error) => Err(FiniteElementError::Upstream(
743 format!("{error}"),
744 format!("{self:?}"),
745 )),
746 }
747 }
748}