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