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