1use crate::{
2 constitutive::solid::viscoelastic::Viscoelastic,
3 fem::block::element::{
4 Element, ElementNodalCoordinates, ElementNodalVelocities, FiniteElement,
5 FiniteElementError, GradientVectors,
6 solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidFiniteElement},
7 surface::{SurfaceElement, SurfaceFiniteElement},
8 },
9 math::{ContractSecondFourthIndicesWithFirstIndicesOf, IDENTITY, Scalar, Tensor},
10 mechanics::{FirstPiolaKirchhoffRateTangentStiffnesses, FirstPiolaKirchhoffStressList},
11};
12
13pub trait ViscoelasticFiniteElement<
14 C,
15 const G: usize,
16 const M: usize,
17 const N: usize,
18 const P: usize,
19> where
20 C: Viscoelastic,
21 Self: SolidFiniteElement<G, M, N, P>,
22{
23 fn nodal_forces(
24 &self,
25 constitutive_model: &C,
26 nodal_coordinates: &ElementNodalCoordinates<N>,
27 nodal_velocities: &ElementNodalVelocities<N>,
28 ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError>;
29 fn nodal_stiffnesses(
30 &self,
31 constitutive_model: &C,
32 nodal_coordinates: &ElementNodalCoordinates<N>,
33 nodal_velocities: &ElementNodalVelocities<N>,
34 ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError>;
35}
36
37impl<C, const G: usize, const N: usize, const O: usize, const P: usize>
38 ViscoelasticFiniteElement<C, G, 3, N, P> for Element<G, N, O>
39where
40 C: Viscoelastic,
41 Self: SolidFiniteElement<G, 3, N, P>,
42{
43 fn nodal_forces(
44 &self,
45 constitutive_model: &C,
46 nodal_coordinates: &ElementNodalCoordinates<N>,
47 nodal_velocities: &ElementNodalVelocities<N>,
48 ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
49 nodal_forces::<_, _, _, _, _, O, _>(
50 self,
51 constitutive_model,
52 self.gradient_vectors(),
53 nodal_coordinates,
54 nodal_velocities,
55 )
56 }
57 fn nodal_stiffnesses(
58 &self,
59 constitutive_model: &C,
60 nodal_coordinates: &ElementNodalCoordinates<N>,
61 nodal_velocities: &ElementNodalVelocities<N>,
62 ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
63 match self
64 .deformation_gradients(nodal_coordinates)
65 .iter()
66 .zip(
67 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
68 .iter(),
69 )
70 .map(|(deformation_gradient, deformation_gradient_rate)| {
71 constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(
72 deformation_gradient,
73 deformation_gradient_rate,
74 )
75 })
76 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
77 {
78 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => {
79 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
80 .iter()
81 .zip(
82 self.gradient_vectors().iter().zip(
83 self.gradient_vectors()
84 .iter()
85 .zip(self.integration_weights()),
86 ),
87 )
88 .map(
89 |(
90 first_piola_kirchhoff_rate_tangent_stiffness,
91 (gradient_vectors_a, (gradient_vectors_b, integration_weight)),
92 )| {
93 gradient_vectors_a
94 .iter()
95 .map(|gradient_vector_a| {
96 gradient_vectors_b
97 .iter()
98 .map(|gradient_vector_b| {
99 first_piola_kirchhoff_rate_tangent_stiffness
100 .contract_second_fourth_indices_with_first_indices_of(
101 gradient_vector_a,
102 gradient_vector_b,
103 )
104 * integration_weight
105 })
106 .collect()
107 })
108 .collect()
109 },
110 )
111 .sum())
112 }
113 Err(error) => Err(FiniteElementError::Upstream(
114 format!("{error}"),
115 format!("{self:?}"),
116 )),
117 }
118 }
119}
120
121impl<C, const G: usize, const N: usize, const O: usize> ViscoelasticFiniteElement<C, G, 2, N, N>
122 for SurfaceElement<G, N, O>
123where
124 C: Viscoelastic,
125 Self: SolidFiniteElement<G, 2, N, N>,
126{
127 fn nodal_forces(
128 &self,
129 constitutive_model: &C,
130 nodal_coordinates: &ElementNodalCoordinates<N>,
131 nodal_velocities: &ElementNodalVelocities<N>,
132 ) -> Result<ElementNodalForcesSolid<N>, FiniteElementError> {
133 nodal_forces::<_, _, _, _, _, O, _>(
134 self,
135 constitutive_model,
136 self.gradient_vectors(),
137 nodal_coordinates,
138 nodal_velocities,
139 )
140 }
141 fn nodal_stiffnesses(
142 &self,
143 constitutive_model: &C,
144 nodal_coordinates: &ElementNodalCoordinates<N>,
145 nodal_velocities: &ElementNodalVelocities<N>,
146 ) -> Result<ElementNodalStiffnessesSolid<N>, FiniteElementError> {
147 match self.deformation_gradients(nodal_coordinates).iter().zip(self.deformation_gradient_rates(nodal_coordinates, nodal_velocities).iter())
148 .map(|(deformation_gradient, deformation_gradient_rate)| {
149 constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)
150 })
151 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses<G>, _>>()
152 {
153 Ok(first_piola_kirchhoff_rate_tangent_stiffnesses) => Ok(first_piola_kirchhoff_rate_tangent_stiffnesses
154 .iter()
155 .zip(
156 self.gradient_vectors()
157 .iter()
158 .zip(self.integration_weights().iter()
159 .zip(self.reference_normals().iter()
160 .zip(Self::normal_gradients(nodal_coordinates))
161 )
162 ),
163 )
164 .map(
165 |(
166 first_piola_kirchoff_rate_tangent_stiffness_mjkl,
167 (gradient_vectors, (integration_weight, (reference_normal, normal_gradients))),
168 )| {
169 gradient_vectors.iter()
170 .map(|gradient_vector_a|
171 gradient_vectors.iter()
172 .zip(normal_gradients.iter())
173 .map(|(gradient_vector_b, normal_gradient_b)|
174 first_piola_kirchoff_rate_tangent_stiffness_mjkl.iter()
175 .map(|first_piola_kirchhoff_rate_tangent_stiffness_m|
176 IDENTITY.iter()
177 .zip(normal_gradient_b.iter())
178 .map(|(identity_n, normal_gradient_b_n)|
179 first_piola_kirchhoff_rate_tangent_stiffness_m.iter()
180 .zip(gradient_vector_a.iter())
181 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mj, gradient_vector_a_j)|
182 first_piola_kirchhoff_rate_tangent_stiffness_mj.iter()
183 .zip(identity_n.iter()
184 .zip(normal_gradient_b_n.iter()))
185 .map(|(first_piola_kirchhoff_rate_tangent_stiffness_mjk, (identity_nk, normal_gradient_b_n_k))|
186 first_piola_kirchhoff_rate_tangent_stiffness_mjk.iter()
187 .zip(gradient_vector_b.iter()
188 .zip(reference_normal.iter()))
189 .map(|(first_piola_kirchoff_rate_tangent_stiffness_mjkl, (gradient_vector_b_l, reference_normal_l))|
190 first_piola_kirchoff_rate_tangent_stiffness_mjkl * gradient_vector_a_j * (
191 identity_nk * gradient_vector_b_l + normal_gradient_b_n_k * reference_normal_l
192 ) * integration_weight
193 ).sum::<Scalar>()
194 ).sum::<Scalar>()
195 ).sum::<Scalar>()
196 ).collect()
197 ).collect()
198 ).collect()
199 ).collect()
200 }
201 )
202 .sum()),
203 Err(error) => Err(FiniteElementError::Upstream(
204 format!("{error}"),
205 format!("{self:?}"),
206 )),
207 }
208 }
209}
210
211fn nodal_forces<
212 C,
213 F,
214 const G: usize,
215 const M: usize,
216 const N: usize,
217 const O: usize,
218 const P: usize,
219>(
220 element: &F,
221 constitutive_model: &C,
222 gradient_vectors: &GradientVectors<G, N>,
223 nodal_coordinates: &ElementNodalCoordinates<N>,
224 nodal_velocities: &ElementNodalVelocities<N>,
225) -> Result<ElementNodalForcesSolid<N>, FiniteElementError>
226where
227 C: Viscoelastic,
228 F: SolidFiniteElement<G, M, N, P>,
229{
230 match element
231 .deformation_gradients(nodal_coordinates)
232 .iter()
233 .zip(
234 element
235 .deformation_gradient_rates(nodal_coordinates, nodal_velocities)
236 .iter(),
237 )
238 .map(|(deformation_gradient, deformation_gradient_rate)| {
239 constitutive_model
240 .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
241 })
242 .collect::<Result<FirstPiolaKirchhoffStressList<G>, _>>()
243 {
244 Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
245 .iter()
246 .zip(gradient_vectors.iter().zip(element.integration_weights()))
247 .map(
248 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
249 gradient_vectors
250 .iter()
251 .map(|gradient_vector| {
252 (first_piola_kirchhoff_stress * gradient_vector) * integration_weight
253 })
254 .collect()
255 },
256 )
257 .sum()),
258 Err(error) => Err(FiniteElementError::Upstream(
259 format!("{error}"),
260 format!("{element:?}"),
261 )),
262 }
263}