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

1use crate::{
2    constitutive::solid::elastic::Elastic,
3    fem::block::element::{FiniteElementError, solid::elastic::ElasticFiniteElement},
4    math::{
5        ContractSecondFourthIndicesWithFirstIndicesOf, Scalar, Tensor, TensorArray, TensorRank2,
6        TensorRank2Vec,
7    },
8    mechanics::{FirstPiolaKirchhoffStresses, FirstPiolaKirchhoffTangentStiffnesses},
9    vem::block::element::{
10        Element, ElementNodalCoordinates, VirtualElement, VirtualElementError,
11        solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidVirtualElement},
12    },
13};
14
15pub trait ElasticVirtualElement<C>
16where
17    C: Elastic,
18    Self: SolidVirtualElement,
19{
20    fn nodal_forces<'a>(
21        &'a self,
22        constitutive_model: &'a C,
23        nodal_coordinates: ElementNodalCoordinates<'a>,
24    ) -> Result<ElementNodalForcesSolid, VirtualElementError>;
25    fn nodal_stiffnesses<'a>(
26        &'a self,
27        constitutive_model: &'a C,
28        nodal_coordinates: ElementNodalCoordinates<'a>,
29    ) -> Result<ElementNodalStiffnessesSolid, VirtualElementError>;
30}
31
32impl<C> ElasticVirtualElement<C> for Element
33where
34    C: Elastic,
35{
36    fn nodal_forces<'a>(
37        &'a self,
38        constitutive_model: &'a C,
39        nodal_coordinates: ElementNodalCoordinates<'a>,
40    ) -> Result<ElementNodalForcesSolid, VirtualElementError> {
41        let mut tetrahedra_forces =
42            ElementNodalForcesSolid::from(vec![[0.0; 3]; nodal_coordinates.len()]);
43        let num_nodes = nodal_coordinates.len() as Scalar;
44        match self
45            .tetrahedra()
46            .iter()
47            .zip(self.tetrahedra_coordinates(&nodal_coordinates).iter())
48            .zip(self.tetrahedra_nodes.iter())
49            .try_for_each(
50                |((tetrahedron, tetrahedron_coordinates), &[face, node_b, node_a])| {
51                    let num_nodes_face = self.faces_nodes()[face].len() as Scalar;
52                    let nodal_forces =
53                        tetrahedron.nodal_forces(constitutive_model, tetrahedron_coordinates)?;
54                    self.faces_nodes()[face].iter().for_each(|&face_node| {
55                        tetrahedra_forces[face_node] += &nodal_forces[0] / num_nodes_face;
56                    });
57                    tetrahedra_forces[node_b] += &nodal_forces[1];
58                    tetrahedra_forces[node_a] += &nodal_forces[2];
59                    tetrahedra_forces.iter_mut().for_each(|entry| {
60                        *entry += &nodal_forces[3] / num_nodes;
61                    });
62                    Ok::<(), FiniteElementError>(())
63                },
64            ) {
65            Ok(()) => {
66                match self
67                    .deformation_gradients(nodal_coordinates)
68                    .iter()
69                    .map(|deformation_gradient| {
70                        constitutive_model.first_piola_kirchhoff_stress(deformation_gradient)
71                    })
72                    .collect::<Result<FirstPiolaKirchhoffStresses, _>>()
73                {
74                    Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
75                        .iter()
76                        .zip(
77                            self.gradient_vectors()
78                                .iter()
79                                .zip(self.integration_weights()),
80                        )
81                        .map(
82                            |(
83                                first_piola_kirchhoff_stress,
84                                (gradient_vectors, integration_weight),
85                            )| {
86                                gradient_vectors
87                                    .iter()
88                                    .map(|gradient_vector| {
89                                        (first_piola_kirchhoff_stress * gradient_vector)
90                                            * integration_weight
91                                    })
92                                    .collect()
93                            },
94                        )
95                        .sum::<ElementNodalForcesSolid>()
96                        * (1.0 - self.stabilization())
97                        + tetrahedra_forces * self.stabilization()),
98                    Err(error) => Err(VirtualElementError::Upstream(
99                        format!("{error}"),
100                        format!("{self:?}"),
101                    )),
102                }
103            }
104            Err(error) => Err(VirtualElementError::Upstream(
105                format!("{error}"),
106                format!("{self:?}"),
107            )),
108        }
109    }
110    fn nodal_stiffnesses<'a>(
111        &'a self,
112        constitutive_model: &'a C,
113        nodal_coordinates: ElementNodalCoordinates<'a>,
114    ) -> Result<ElementNodalStiffnessesSolid, VirtualElementError> {
115        let mut tetrahedra_stiffnesses = ElementNodalStiffnessesSolid::from(vec![
116                TensorRank2Vec::from(vec![
117                    TensorRank2::zero(); nodal_coordinates.len()
118                ]);
119                nodal_coordinates.len()
120            ]);
121        let num_nodes = nodal_coordinates.len() as Scalar;
122        match self
123            .tetrahedra()
124            .iter()
125            .zip(self.tetrahedra_coordinates(&nodal_coordinates).iter())
126            .zip(self.tetrahedra_nodes.iter())
127            .try_for_each(
128                |((tetrahedron, tetrahedron_coordinates), &[face, node_b, node_a])| {
129                    let num_nodes_face = self.faces_nodes()[face].len() as Scalar;
130                    let nodal_stiffnesses = tetrahedron
131                        .nodal_stiffnesses(constitutive_model, tetrahedron_coordinates)?;
132                    self.faces_nodes()[face].iter().for_each(|&face_node_1| {
133                        tetrahedra_stiffnesses[face_node_1][node_b] +=
134                            &nodal_stiffnesses[0][1] / num_nodes_face;
135                        tetrahedra_stiffnesses[face_node_1][node_a] +=
136                            &nodal_stiffnesses[0][2] / num_nodes_face;
137                        tetrahedra_stiffnesses[node_b][face_node_1] +=
138                            &nodal_stiffnesses[1][0] / num_nodes_face;
139                        tetrahedra_stiffnesses[node_a][face_node_1] +=
140                            &nodal_stiffnesses[2][0] / num_nodes_face;
141                        tetrahedra_stiffnesses[face_node_1]
142                            .iter_mut()
143                            .for_each(|entry| {
144                                *entry += &nodal_stiffnesses[0][3] / num_nodes_face / num_nodes;
145                            });
146                        self.faces_nodes()[face].iter().for_each(|&face_node_2| {
147                            tetrahedra_stiffnesses[face_node_1][face_node_2] +=
148                                &nodal_stiffnesses[0][0] / num_nodes_face.powi(2);
149                        })
150                    });
151                    tetrahedra_stiffnesses[node_b][node_b] += &nodal_stiffnesses[1][1];
152                    tetrahedra_stiffnesses[node_b][node_a] += &nodal_stiffnesses[1][2];
153                    tetrahedra_stiffnesses[node_a][node_b] += &nodal_stiffnesses[2][1];
154                    tetrahedra_stiffnesses[node_a][node_a] += &nodal_stiffnesses[2][2];
155                    tetrahedra_stiffnesses
156                        .iter_mut()
157                        .for_each(|tetrahedra_stiffness| {
158                            tetrahedra_stiffness[node_b] += &nodal_stiffnesses[3][1] / num_nodes;
159                            tetrahedra_stiffness[node_a] += &nodal_stiffnesses[3][2] / num_nodes;
160                            self.faces_nodes()[face].iter().for_each(|&face_node| {
161                                tetrahedra_stiffness[face_node] +=
162                                    &nodal_stiffnesses[3][0] / num_nodes_face / num_nodes;
163                            });
164                            tetrahedra_stiffness.iter_mut().for_each(|entry| {
165                                *entry += &nodal_stiffnesses[3][3] / num_nodes.powi(2);
166                            })
167                        });
168                    tetrahedra_stiffnesses[node_b].iter_mut().for_each(|entry| {
169                        *entry += &nodal_stiffnesses[1][3] / num_nodes;
170                    });
171                    tetrahedra_stiffnesses[node_a].iter_mut().for_each(|entry| {
172                        *entry += &nodal_stiffnesses[2][3] / num_nodes;
173                    });
174                    Ok::<(), FiniteElementError>(())
175                },
176            ) {
177            Ok(()) => {
178                match self
179                    .deformation_gradients(nodal_coordinates)
180                    .iter()
181                    .map(|deformation_gradient| {
182                        constitutive_model.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)
183                    })
184                    .collect::<Result<FirstPiolaKirchhoffTangentStiffnesses, _>>()
185                {
186                    Ok(first_piola_kirchhoff_tangent_stiffnesses) => {
187                        Ok(first_piola_kirchhoff_tangent_stiffnesses
188                            .iter()
189                            .zip(
190                                self.gradient_vectors()
191                                    .iter()
192                                    .zip(self.integration_weights()),
193                            )
194                            .map(
195                                |(
196                                    first_piola_kirchhoff_tangent_stiffness,
197                                    (gradient_vectors, integration_weight),
198                                )| {
199                                    gradient_vectors
200                                        .iter()
201                                        .map(|gradient_vector_a| {
202                                            gradient_vectors
203                                                .iter()
204                                                .map(|gradient_vector_b| {
205                                                    first_piola_kirchhoff_tangent_stiffness
206                                                    .contract_second_fourth_indices_with_first_indices_of(
207                                                        gradient_vector_a,
208                                                        gradient_vector_b,
209                                                    )
210                                                    * integration_weight
211                                                })
212                                                .collect()
213                                        })
214                                        .collect()
215                                },
216                            )
217                            .sum::<ElementNodalStiffnessesSolid>()
218                            * (1.0 - self.stabilization())
219                            + tetrahedra_stiffnesses * self.stabilization())
220                    }
221                    Err(error) => Err(VirtualElementError::Upstream(
222                        format!("{error}"),
223                        format!("{self:?}"),
224                    )),
225                }
226            }
227            Err(error) => Err(VirtualElementError::Upstream(
228                format!("{error}"),
229                format!("{self:?}"),
230            )),
231        }
232    }
233}
conspire/domain/vem/block/element/solid/elastic/mod.rs

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