conspire/domain/vem/block/element/solid/viscoelastic/
mod.rs1use crate::{
2 constitutive::solid::viscoelastic::Viscoelastic,
3 domain::block::element::solid::viscoelastic::ViscoelasticElement,
4 fem::block::element::FiniteElementError,
5 math::{ContractSecondFourthWithFirst, Scalar, Tensor, TensorArray},
6 mechanics::{
7 Damping, FirstPiolaKirchhoffRateTangentStiffnesses, FirstPiolaKirchhoffStresses, Force,
8 },
9 vem::block::element::{
10 Element, ElementNodalCoordinates, ElementNodalVelocities, VirtualElement,
11 VirtualElementError,
12 solid::{
13 ElementNodalDampingsSolid, ElementNodalForcesSolid, SolidElement, SolidVirtualElement,
14 },
15 },
16};
17
18pub trait ViscoelasticVirtualElement<C>
19where
20 C: Viscoelastic,
21 Self: SolidVirtualElement
22 + ViscoelasticElement<
23 C,
24 Forces = ElementNodalForcesSolid,
25 Dampings = ElementNodalDampingsSolid,
26 Error = VirtualElementError,
27 >,
28{
29}
30
31impl<T, C> ViscoelasticVirtualElement<C> for T
32where
33 C: Viscoelastic,
34 T: SolidVirtualElement
35 + ViscoelasticElement<
36 C,
37 Forces = ElementNodalForcesSolid,
38 Dampings = ElementNodalDampingsSolid,
39 Error = VirtualElementError,
40 >,
41{
42}
43
44impl<C> ViscoelasticElement<C> for Element
45where
46 C: Viscoelastic,
47{
48 type Forces = ElementNodalForcesSolid;
49 type Dampings = ElementNodalDampingsSolid;
50 type Error = VirtualElementError;
51 fn nodal_forces(
52 &self,
53 constitutive_model: &C,
54 nodal_coordinates: &ElementNodalCoordinates,
55 nodal_velocities: &ElementNodalVelocities,
56 ) -> Result<ElementNodalForcesSolid, VirtualElementError> {
57 let stabilization = self.stabilization();
58 let inverse_num_nodes = 1.0 / nodal_coordinates.len() as Scalar;
59 let tetrahedra_coordinates = self.tetrahedra_coordinates(nodal_coordinates);
60 let tetrahedra_velocities = self.tetrahedra_coordinates(nodal_velocities);
61 let mut forces = self
62 .deformation_gradients(nodal_coordinates)
63 .iter()
64 .zip(
65 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
66 .iter(),
67 )
68 .map(|(deformation_gradient, deformation_gradient_rate)| {
69 constitutive_model
70 .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)
71 })
72 .collect::<Result<FirstPiolaKirchhoffStresses, _>>()
73 .map_err(|error| self.upstream(error))?
74 .iter()
75 .zip(
76 self.gradient_vectors()
77 .iter()
78 .zip(self.integration_weights()),
79 )
80 .map(
81 |(first_piola_kirchhoff_stress, (gradient_vectors, integration_weight))| {
82 gradient_vectors
83 .iter()
84 .map(|gradient_vector| {
85 (first_piola_kirchhoff_stress * gradient_vector)
86 * (integration_weight * (1.0 - stabilization))
87 })
88 .collect()
89 },
90 )
91 .sum::<ElementNodalForcesSolid>();
92 let mut faces_forces = vec![Force::zero(); self.faces_nodes().len()];
93 let mut center_force = Force::zero();
94 self.tetrahedra()
95 .iter()
96 .zip(
97 tetrahedra_coordinates
98 .iter()
99 .zip(tetrahedra_velocities.iter()),
100 )
101 .zip(self.tetrahedra_nodes().iter())
102 .try_for_each(
103 |(
104 (tetrahedron, (tetrahedron_coordinates, tetrahedron_velocities)),
105 &[face, node_b, node_a],
106 )| {
107 let nodal_forces = tetrahedron.nodal_forces(
108 constitutive_model,
109 tetrahedron_coordinates,
110 tetrahedron_velocities,
111 )?;
112 faces_forces[face] += &nodal_forces[0];
113 forces[node_b] += &nodal_forces[1] * stabilization;
114 forces[node_a] += &nodal_forces[2] * stabilization;
115 center_force += &nodal_forces[3];
116 Ok::<(), FiniteElementError>(())
117 },
118 )
119 .map_err(|error| self.upstream(error))?;
120 self.faces_nodes()
121 .iter()
122 .zip(faces_forces.iter())
123 .for_each(|(face_nodes, face_force)| {
124 let face_force = face_force * (stabilization / face_nodes.len() as Scalar);
125 face_nodes
126 .iter()
127 .for_each(|&face_node| forces[face_node] += &face_force)
128 });
129 center_force *= stabilization * inverse_num_nodes;
130 forces.iter_mut().for_each(|force| *force += ¢er_force);
131 Ok(forces)
132 }
133 fn nodal_stiffnesses(
134 &self,
135 constitutive_model: &C,
136 nodal_coordinates: &ElementNodalCoordinates,
137 nodal_velocities: &ElementNodalVelocities,
138 ) -> Result<ElementNodalDampingsSolid, VirtualElementError> {
139 let num_nodes = nodal_coordinates.len();
140 let stabilization = self.stabilization();
141 let inverse_num_nodes = 1.0 / num_nodes as Scalar;
142 let tetrahedra_coordinates = self.tetrahedra_coordinates(nodal_coordinates);
143 let tetrahedra_velocities = self.tetrahedra_coordinates(nodal_velocities);
144 let mut stiffnesses = self
145 .deformation_gradients(nodal_coordinates)
146 .iter()
147 .zip(
148 self.deformation_gradient_rates(nodal_coordinates, nodal_velocities)
149 .iter(),
150 )
151 .map(|(deformation_gradient, deformation_gradient_rate)| {
152 constitutive_model.first_piola_kirchhoff_rate_tangent_stiffness(
153 deformation_gradient,
154 deformation_gradient_rate,
155 )
156 })
157 .collect::<Result<FirstPiolaKirchhoffRateTangentStiffnesses, _>>()
158 .map_err(|error| self.upstream(error))?
159 .iter()
160 .zip(
161 self.gradient_vectors()
162 .iter()
163 .zip(self.integration_weights()),
164 )
165 .map(
166 |(
167 first_piola_kirchhoff_rate_tangent_stiffness,
168 (gradient_vectors, integration_weight),
169 )| {
170 let weight = integration_weight * (1.0 - stabilization);
171 gradient_vectors
172 .iter()
173 .map(|gradient_vector_a| {
174 gradient_vectors
175 .iter()
176 .map(|gradient_vector_b| {
177 first_piola_kirchhoff_rate_tangent_stiffness
178 .contract_second_fourth_with_first(
179 gradient_vector_a,
180 gradient_vector_b,
181 )
182 * weight
183 })
184 .collect()
185 })
186 .collect()
187 },
188 )
189 .sum::<ElementNodalDampingsSolid>();
190 let num_faces = self.faces_nodes().len();
191 let mut faces_stiffnesses = vec![Damping::zero(); num_faces];
192 let mut faces_rows = vec![Damping::zero(); num_faces];
193 let mut faces_columns = vec![Damping::zero(); num_faces];
194 let mut rows = vec![Damping::zero(); num_nodes];
195 let mut columns = vec![Damping::zero(); num_nodes];
196 let mut center_stiffness = Damping::zero();
197 self.tetrahedra()
198 .iter()
199 .zip(
200 tetrahedra_coordinates
201 .iter()
202 .zip(tetrahedra_velocities.iter()),
203 )
204 .zip(self.tetrahedra_nodes().iter())
205 .try_for_each(
206 |(
207 (tetrahedron, (tetrahedron_coordinates, tetrahedron_velocities)),
208 &[face, node_b, node_a],
209 )| {
210 let nodal_stiffnesses = tetrahedron.nodal_stiffnesses(
211 constitutive_model,
212 tetrahedron_coordinates,
213 tetrahedron_velocities,
214 )?;
215 let face_nodes = &self.faces_nodes()[face];
216 let weight = stabilization / face_nodes.len() as Scalar;
217 faces_stiffnesses[face] += &nodal_stiffnesses[0][0];
218 faces_rows[face] += &nodal_stiffnesses[0][3];
219 faces_columns[face] += &nodal_stiffnesses[3][0];
220 let face_node_b = &nodal_stiffnesses[0][1] * weight;
221 let face_node_a = &nodal_stiffnesses[0][2] * weight;
222 let node_b_face = &nodal_stiffnesses[1][0] * weight;
223 let node_a_face = &nodal_stiffnesses[2][0] * weight;
224 face_nodes.iter().for_each(|&face_node| {
225 stiffnesses[face_node][node_b] += &face_node_b;
226 stiffnesses[face_node][node_a] += &face_node_a;
227 stiffnesses[node_b][face_node] += &node_b_face;
228 stiffnesses[node_a][face_node] += &node_a_face;
229 });
230 stiffnesses[node_b][node_b] += &nodal_stiffnesses[1][1] * stabilization;
231 stiffnesses[node_b][node_a] += &nodal_stiffnesses[1][2] * stabilization;
232 stiffnesses[node_a][node_b] += &nodal_stiffnesses[2][1] * stabilization;
233 stiffnesses[node_a][node_a] += &nodal_stiffnesses[2][2] * stabilization;
234 rows[node_b] += &nodal_stiffnesses[1][3] * (stabilization * inverse_num_nodes);
235 rows[node_a] += &nodal_stiffnesses[2][3] * (stabilization * inverse_num_nodes);
236 columns[node_b] +=
237 &nodal_stiffnesses[3][1] * (stabilization * inverse_num_nodes);
238 columns[node_a] +=
239 &nodal_stiffnesses[3][2] * (stabilization * inverse_num_nodes);
240 center_stiffness += &nodal_stiffnesses[3][3]
241 * (stabilization * inverse_num_nodes * inverse_num_nodes);
242 Ok::<(), FiniteElementError>(())
243 },
244 )
245 .map_err(|error| self.upstream(error))?;
246 self.faces_nodes()
247 .iter()
248 .zip(
249 faces_stiffnesses
250 .iter()
251 .zip(faces_rows.iter().zip(faces_columns.iter())),
252 )
253 .for_each(|(face_nodes, (face_stiffness, (face_row, face_column)))| {
254 let inverse_num_nodes_face = 1.0 / face_nodes.len() as Scalar;
255 let face_stiffness = face_stiffness
256 * (stabilization * inverse_num_nodes_face * inverse_num_nodes_face);
257 let face_row =
258 face_row * (stabilization * inverse_num_nodes_face * inverse_num_nodes);
259 let face_column =
260 face_column * (stabilization * inverse_num_nodes_face * inverse_num_nodes);
261 face_nodes.iter().for_each(|&face_node_a| {
262 rows[face_node_a] += &face_row;
263 columns[face_node_a] += &face_column;
264 face_nodes.iter().for_each(|&face_node_b| {
265 stiffnesses[face_node_a][face_node_b] += &face_stiffness
266 })
267 })
268 });
269 rows.iter_mut().for_each(|row| *row += ¢er_stiffness);
270 stiffnesses
271 .iter_mut()
272 .zip(rows.iter())
273 .for_each(|(stiffness, row)| {
274 stiffness
275 .iter_mut()
276 .zip(columns.iter())
277 .for_each(|(entry, column)| {
278 *entry += row;
279 *entry += column
280 })
281 });
282 Ok(stiffnesses)
283 }
284}