conspire/domain/vem/block/element/solid/elastic_viscoplastic/
mod.rs

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

conspire/domain/vem/block/element/solid/elastic_viscoplastic/
mod.rs
