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        ContractSecondFourthWithFirst, Scalar, Tensor, TensorArray, TensorRank2, TensorRank2Vec,
6    },
7    mechanics::{FirstPiolaKirchhoffStresses, FirstPiolaKirchhoffTangentStiffnesses},
8    vem::block::element::{
9        Element, ElementNodalCoordinates, VirtualElement, VirtualElementError,
10        solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid, SolidVirtualElement},
11    },
12};
13
14pub trait ElasticVirtualElement<C>
15where
16    C: Elastic,
17    Self: SolidVirtualElement,
18{
19    fn nodal_forces<'a>(
20        &'a self,
21        constitutive_model: &'a C,
22        nodal_coordinates: ElementNodalCoordinates<'a>,
23    ) -> Result<ElementNodalForcesSolid, VirtualElementError>;
24    fn nodal_stiffnesses<'a>(
25        &'a self,
26        constitutive_model: &'a C,
27        nodal_coordinates: ElementNodalCoordinates<'a>,
28    ) -> Result<ElementNodalStiffnessesSolid, VirtualElementError>;
29}
30
31impl<C> ElasticVirtualElement<C> for Element
32where
33    C: Elastic,
34{
35    fn nodal_forces<'a>(
36        &'a self,
37        constitutive_model: &'a C,
38        nodal_coordinates: ElementNodalCoordinates<'a>,
39    ) -> Result<ElementNodalForcesSolid, VirtualElementError> {
40        let mut tetrahedra_forces =
41            ElementNodalForcesSolid::from(vec![[0.0; 3]; nodal_coordinates.len()]);
42        let num_nodes = nodal_coordinates.len() as Scalar;
43        match self
44            .tetrahedra()
45            .iter()
46            .zip(self.tetrahedra_coordinates(&nodal_coordinates).iter())
47            .zip(self.tetrahedra_nodes.iter())
48            .try_for_each(
49                |((tetrahedron, tetrahedron_coordinates), &[face, node_b, node_a])| {
50                    let num_nodes_face = self.faces_nodes()[face].len() as Scalar;
51                    let nodal_forces =
52                        tetrahedron.nodal_forces(constitutive_model, tetrahedron_coordinates)?;
53                    self.faces_nodes()[face].iter().for_each(|&face_node| {
54                        tetrahedra_forces[face_node] += &nodal_forces[0] / num_nodes_face;
55                    });
56                    tetrahedra_forces[node_b] += &nodal_forces[1];
57                    tetrahedra_forces[node_a] += &nodal_forces[2];
58                    tetrahedra_forces.iter_mut().for_each(|entry| {
59                        *entry += &nodal_forces[3] / num_nodes;
60                    });
61                    Ok::<(), FiniteElementError>(())
62                },
63            ) {
64            Ok(()) => {
65                match self
66                    .deformation_gradients(nodal_coordinates)
67                    .iter()
68                    .map(|deformation_gradient| {
69                        constitutive_model.first_piola_kirchhoff_stress(deformation_gradient)
70                    })
71                    .collect::<Result<FirstPiolaKirchhoffStresses, _>>()
72                {
73                    Ok(first_piola_kirchhoff_stresses) => Ok(first_piola_kirchhoff_stresses
74                        .iter()
75                        .zip(
76                            self.gradient_vectors()
77                                .iter()
78                                .zip(self.integration_weights()),
79                        )
80                        .map(
81                            |(
82                                first_piola_kirchhoff_stress,
83                                (gradient_vectors, integration_weight),
84                            )| {
85                                gradient_vectors
86                                    .iter()
87                                    .map(|gradient_vector| {
88                                        (first_piola_kirchhoff_stress * gradient_vector)
89                                            * integration_weight
90                                    })
91                                    .collect()
92                            },
93                        )
94                        .sum::<ElementNodalForcesSolid>()
95                        * (1.0 - self.stabilization())
96                        + tetrahedra_forces * self.stabilization()),
97                    Err(error) => Err(VirtualElementError::Upstream(
98                        format!("{error}"),
99                        format!("{self:?}"),
100                    )),
101                }
102            }
103            Err(error) => Err(VirtualElementError::Upstream(
104                format!("{error}"),
105                format!("{self:?}"),
106            )),
107        }
108    }
109    fn nodal_stiffnesses<'a>(
110        &'a self,
111        constitutive_model: &'a C,
112        nodal_coordinates: ElementNodalCoordinates<'a>,
113    ) -> Result<ElementNodalStiffnessesSolid, VirtualElementError> {
114        let mut tetrahedra_stiffnesses = ElementNodalStiffnessesSolid::from(vec![
115                TensorRank2Vec::from(vec![
116                    TensorRank2::zero(); nodal_coordinates.len()
117                ]);
118                nodal_coordinates.len()
119            ]);
120        let num_nodes = nodal_coordinates.len() as Scalar;
121        match self
122            .tetrahedra()
123            .iter()
124            .zip(self.tetrahedra_coordinates(&nodal_coordinates).iter())
125            .zip(self.tetrahedra_nodes.iter())
126            .try_for_each(
127                |((tetrahedron, tetrahedron_coordinates), &[face, node_b, node_a])| {
128                    let num_nodes_face = self.faces_nodes()[face].len() as Scalar;
129                    let nodal_stiffnesses = tetrahedron
130                        .nodal_stiffnesses(constitutive_model, tetrahedron_coordinates)?;
131                    self.faces_nodes()[face].iter().for_each(|&face_node_1| {
132                        tetrahedra_stiffnesses[face_node_1][node_b] +=
133                            &nodal_stiffnesses[0][1] / num_nodes_face;
134                        tetrahedra_stiffnesses[face_node_1][node_a] +=
135                            &nodal_stiffnesses[0][2] / num_nodes_face;
136                        tetrahedra_stiffnesses[node_b][face_node_1] +=
137                            &nodal_stiffnesses[1][0] / num_nodes_face;
138                        tetrahedra_stiffnesses[node_a][face_node_1] +=
139                            &nodal_stiffnesses[2][0] / num_nodes_face;
140                        tetrahedra_stiffnesses[face_node_1]
141                            .iter_mut()
142                            .for_each(|entry| {
143                                *entry += &nodal_stiffnesses[0][3] / num_nodes_face / num_nodes;
144                            });
145                        self.faces_nodes()[face].iter().for_each(|&face_node_2| {
146                            tetrahedra_stiffnesses[face_node_1][face_node_2] +=
147                                &nodal_stiffnesses[0][0] / num_nodes_face.powi(2);
148                        })
149                    });
150                    tetrahedra_stiffnesses[node_b][node_b] += &nodal_stiffnesses[1][1];
151                    tetrahedra_stiffnesses[node_b][node_a] += &nodal_stiffnesses[1][2];
152                    tetrahedra_stiffnesses[node_a][node_b] += &nodal_stiffnesses[2][1];
153                    tetrahedra_stiffnesses[node_a][node_a] += &nodal_stiffnesses[2][2];
154                    tetrahedra_stiffnesses
155                        .iter_mut()
156                        .for_each(|tetrahedra_stiffness| {
157                            tetrahedra_stiffness[node_b] += &nodal_stiffnesses[3][1] / num_nodes;
158                            tetrahedra_stiffness[node_a] += &nodal_stiffnesses[3][2] / num_nodes;
159                            self.faces_nodes()[face].iter().for_each(|&face_node| {
160                                tetrahedra_stiffness[face_node] +=
161                                    &nodal_stiffnesses[3][0] / num_nodes_face / num_nodes;
162                            });
163                            tetrahedra_stiffness.iter_mut().for_each(|entry| {
164                                *entry += &nodal_stiffnesses[3][3] / num_nodes.powi(2);
165                            })
166                        });
167                    tetrahedra_stiffnesses[node_b].iter_mut().for_each(|entry| {
168                        *entry += &nodal_stiffnesses[1][3] / num_nodes;
169                    });
170                    tetrahedra_stiffnesses[node_a].iter_mut().for_each(|entry| {
171                        *entry += &nodal_stiffnesses[2][3] / num_nodes;
172                    });
173                    Ok::<(), FiniteElementError>(())
174                },
175            ) {
176            Ok(()) => {
177                match self
178                    .deformation_gradients(nodal_coordinates)
179                    .iter()
180                    .map(|deformation_gradient| {
181                        constitutive_model
182                            .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_with_first(
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
