conspire/domain/fem/block/element/solid/viscoelastic/
mod.rs

1use crate::{
2    constitutive::solid::viscoelastic::Viscoelastic,
3    fem::block::element::{
4        Element, ElementNodalCoordinates, ElementNodalVelocities, FiniteElementError,
5        solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidFiniteElement},
6    },
7    math::{ContractSecondFourthIndicesWithFirstIndicesOf, Tensor},
8    mechanics::{FirstPiolaKirchhoffRateTangentStiffnesses, FirstPiolaKirchhoffStresses},
9};
10
11pub trait ViscoelasticFiniteElement<C, const G: usize, const N: usize>
12where
13    C: Viscoelastic,
14    Self: SolidFiniteElement<G, N>,
15{
16    fn nodal_forces(
17        &self,
18        constitutive_model: &C,
19        nodal_coordinates: &ElementNodalCoordinates<N>,
20        nodal_velocities: &ElementNodalVelocities<N>,
21    ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError>;
22    fn nodal_stiffnesses(
23        &self,
24        constitutive_model: &C,
25        nodal_coordinates: &ElementNodalCoordinates<N>,
26        nodal_velocities: &ElementNodalVelocities<N>,
27    ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError>;
28}
29
30impl<C, const G: usize, const N: usize> ViscoelasticFiniteElement<C, G, N> for Element<G, N>
31where
32    C: Viscoelastic,
33{
34    fn nodal_forces(
35        &self,
36        constitutive_model: &C,
37        nodal_coordinates: &ElementNodalCoordinates<N>,
38        nodal_velocities: &ElementNodalVelocities<N>,
39    ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
40        match self
41            .deformation_gradients(nodal_coordinates)
42            .iter()
43            .zip(
44                self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
45                    .iter(),
46            )
47            .map(|(deformation_gradient, deformation_gradient_rate)| {
48                constitutive_model
49                    .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
50            })
51            .collect::<Result<FirstPiolaKirchhoffStresses<G>, _>>()
52        {
53            Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
54                .iter()
55                .zip(
56                    self.gradient_vectors()
57                        .iter()
58                        .zip(self.integration_weights().iter()),
59                )
60                .map(
61                    |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
62                        gradient_vectors
63                            .iter()
64                            .map(|gradient_vector| {
65                                (first_piola_kirchhoff_stress * gradient_vector)
66                                    * integration_weight
67                            })
68                            .collect()
69                    },
70                )
71                .sum()),
72            Err(error) => Err(FiniteElementError::Upstream(
73                format!("{error}"),
74                format!("{self:?}"),
75            )),
76        }
77    }
78    fn nodal_stiffnesses(
79        &self,
80        constitutive_model: &C,
81        nodal_coordinates: &ElementNodalCoordinates<N>,
82        nodal_velocities: &ElementNodalVelocities<N>,
83    ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
84        match self
85            .deformation_gradients(nodal_coordinates)
86            .iter()
87            .zip(
88                self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
89                    .iter(),
90            )
91            .map(|(deformation_gradient, deformation_gradient_rate)| {
92                constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(
93                    deformation_gradient,
94                    deformation_gradient_rate,
95                )
96            })
97            .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
98        {
99            Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => {
100                Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
101                    .iter()
102                    .zip(
103                        self.gradient_vectors().iter().zip(
104                            self.gradient_vectors()
105                                .iter()
106                                .zip(self.integration_weights().iter()),
107                        ),
108                    )
109                    .map(
110                        |(
111                            first_piola_kirchhoff_rate_tangent_stiffness,
112                            (gradient_vectors_a, (gradient_vectors_b, integration_weight)),
113                        )| {
114                            gradient_vectors_a
115                                .iter()
116                                .map(|gradient_vector_a| {
117                                    gradient_vectors_b
118                                        .iter()
119                                        .map(|gradient_vector_b| {
120                                            first_piola_kirchhoff_rate_tangent_stiffness
121                                        .contract_second_fourth_indices_with_first_indices_of(
122                                            gradient_vector_a,
123                                            gradient_vector_b,
124                                        )
125                                        * integration_weight
126                                        })
127                                        .collect()
128                                })
129                                .collect()
130                        },
131                    )
132                    .sum())
133            }
134            Err(error) => Err(FiniteElementError::Upstream(
135                format!("{error}"),
136                format!("{self:?}"),
137            )),
138        }
139    }
140}