1pub mod solid;
2
3use crate::{
4 fem::block::element::{
5 ElementNodalCoordinates as FemElementNodalCoordinates,
6 ElementNodalReferenceCoordinates as FemElementNodalReferenceCoordinates, FiniteElement,
7 linear::Tetrahedron,
8 },
9 math::{
10 CrossProduct, Scalar, Scalars, Style, StyledError, Tensor, TensorRank1Vec2D,
11 assert::AssertionError, styled_error,
12 },
13 mechanics::{CurrentCoordinate, CurrentCoordinatesRef, ReferenceCoordinate, Vectors2D},
14 vem::{NodalCoordinates, NodalReferenceCoordinates},
15};
16
17#[cfg(test)]
18use crate::math::assert::Assert;
19use std::{
20 collections::VecDeque,
21 fmt::{self, Debug, Formatter},
22};
23
24pub type ElementNodalCoordinates<'a> = CurrentCoordinatesRef<'a>;
25pub type ElementNodalReferenceCoordinates = TensorRank1Vec2D<3, 0>;
26pub type GradientVectors = Vectors2D<0>;
27
28pub type TetrahedraCoordinates = Vec<FemElementNodalCoordinates<4>>;
29
30pub struct Element {
31 faces_nodes: Vec<Vec<usize>>,
32 gradient_vectors: GradientVectors,
33 integration_weights: Scalars,
34 stabilization: Scalar,
35 tetrahedra: Vec<Tetrahedron>,
36 tetrahedra_nodes: Vec<[usize; 3]>,
37}
38
39pub trait VirtualElement
40where
41 for<'a> Self: From<(
42 ElementNodalReferenceCoordinates,
43 &'a [usize],
44 &'a [usize],
45 &'a [Vec<usize>],
46 )>,
47{
48 fn element_center<'a>(nodal_coordinates: &ElementNodalCoordinates<'a>) -> CurrentCoordinate;
49 fn faces_centers<'a>(
50 &'a self,
51 nodal_coordinates: &ElementNodalCoordinates<'a>,
52 ) -> NodalCoordinates;
53 fn faces_nodes(&self) -> &[Vec<usize>];
54 fn gradient_vectors(&self) -> &GradientVectors;
55 fn integration_weights(&self) -> &Scalars;
56 fn stabilization(&self) -> Scalar;
57 fn tetrahedra(&self) -> &[Tetrahedron];
58 fn tetrahedra_coordinates<'a>(
59 &'a self,
60 nodal_coordinates: &ElementNodalCoordinates<'a>,
61 ) -> TetrahedraCoordinates;
62 fn tetrahedra_nodes(&self) -> &[[usize; 3]];
63}
64
65impl VirtualElement for Element {
66 fn element_center<'a>(nodal_coordinates: &ElementNodalCoordinates<'a>) -> CurrentCoordinate {
67 nodal_coordinates
68 .iter()
69 .map(|&nodal_coordinate| nodal_coordinate.clone())
70 .sum::<CurrentCoordinate>()
71 / nodal_coordinates.len() as Scalar
72 }
73 fn faces_centers<'a>(
74 &'a self,
75 nodal_coordinates: &ElementNodalCoordinates<'a>,
76 ) -> NodalCoordinates {
77 self.faces_nodes()
78 .iter()
79 .map(|face_nodes| {
80 face_nodes
81 .iter()
82 .map(|&face_node| nodal_coordinates[face_node].clone())
83 .sum::<CurrentCoordinate>()
84 / (face_nodes.len() as Scalar)
85 })
86 .collect()
87 }
88 fn faces_nodes(&self) -> &[Vec<usize>] {
89 &self.faces_nodes
90 }
91 fn gradient_vectors(&self) -> &GradientVectors {
92 &self.gradient_vectors
93 }
94 fn integration_weights(&self) -> &Scalars {
95 &self.integration_weights
96 }
97 fn stabilization(&self) -> Scalar {
98 self.stabilization
99 }
100 fn tetrahedra(&self) -> &[Tetrahedron] {
101 &self.tetrahedra
102 }
103 fn tetrahedra_coordinates<'a>(
104 &'a self,
105 nodal_coordinates: &ElementNodalCoordinates<'a>,
106 ) -> TetrahedraCoordinates {
107 let element_center = Self::element_center(nodal_coordinates);
108 let faces_centers = self.faces_centers(nodal_coordinates);
109 self.tetrahedra_nodes()
110 .iter()
111 .map(|&[face, node_b, node_a]| {
112 [
113 faces_centers[face].clone(),
114 nodal_coordinates[node_b].clone(),
115 nodal_coordinates[node_a].clone(),
116 element_center.clone(),
117 ]
118 .into()
119 })
120 .collect()
121 }
122 fn tetrahedra_nodes(&self) -> &[[usize; 3]] {
123 &self.tetrahedra_nodes
124 }
125}
126
127impl
128 From<(
129 ElementNodalReferenceCoordinates,
130 &[usize],
131 &[usize],
132 &[Vec<usize>],
133 )> for Element
134{
135 fn from(
136 (reference_nodal_coordinates, element_faces, element_nodes, block_faces_nodes): (
137 ElementNodalReferenceCoordinates,
138 &[usize],
139 &[usize],
140 &[Vec<usize>],
141 ),
142 ) -> Self {
143 let faces_nodes = element_faces
144 .iter()
145 .map(|&element_face| {
146 block_faces_nodes[element_face]
147 .iter()
148 .map(|face_node| {
149 element_nodes
150 .iter()
151 .position(|element_node| face_node == element_node)
152 .unwrap()
153 })
154 .collect::<Vec<_>>()
155 })
156 .collect::<Vec<_>>();
157 let mut nodal_coordinates =
158 NodalReferenceCoordinates::from(vec![
159 ReferenceCoordinate::from([0.0, 0.0, 0.0]);
160 element_nodes.len()
161 ]);
162 block_faces_nodes
163 .iter()
164 .flatten()
165 .zip(reference_nodal_coordinates.iter().flatten())
166 .for_each(|(&node, coordinates)| nodal_coordinates[node] = coordinates.clone());
167 let element_center = nodal_coordinates.into_iter().sum::<ReferenceCoordinate>()
168 / (element_nodes.len() as Scalar);
169 let tetrahedra = reference_nodal_coordinates
170 .iter()
171 .flat_map(|face_coordinates| {
172 let face_center = face_coordinates
173 .iter()
174 .cloned()
175 .sum::<ReferenceCoordinate>()
176 / (face_coordinates.len() as Scalar);
177 let mut face_coordinates_one_ahead = VecDeque::from(face_coordinates.clone());
178 let first_entry = face_coordinates_one_ahead.pop_front().unwrap();
179 face_coordinates_one_ahead.push_back(first_entry);
180 face_coordinates
181 .iter()
182 .zip(face_coordinates_one_ahead)
183 .map(|(node_a_coordinates, node_b_coordinates)| {
184 Tetrahedron::from(FemElementNodalReferenceCoordinates::from([
185 face_center.clone(),
186 node_b_coordinates,
187 node_a_coordinates.clone(),
188 element_center.clone(),
189 ]))
190 })
191 .collect::<Vec<_>>()
192 })
193 .collect::<Vec<_>>();
194 let tetrahedra_nodes = faces_nodes
195 .iter()
196 .enumerate()
197 .flat_map(|(face, face_nodes)| {
198 let mut face_nodes_one_ahead = VecDeque::from(face_nodes.clone());
199 let first_entry = face_nodes_one_ahead.pop_front().unwrap();
200 face_nodes_one_ahead.push_back(first_entry);
201 face_nodes
202 .iter()
203 .zip(face_nodes_one_ahead)
204 .map(|(&node_a, node_b)| [face, node_b, node_a])
205 .collect::<Vec<_>>()
206 })
207 .collect::<Vec<_>>();
208 let element_volume = tetrahedra
209 .iter()
210 .map(|tetrahedron| tetrahedron.volume())
211 .sum();
212 let integration_weights = Scalars::from([element_volume]);
213 let gradient_vectors = vec![
214 element_nodes
215 .iter()
216 .map(|&node| {
217 element_faces
218 .iter()
219 .zip(reference_nodal_coordinates.iter())
220 .filter_map(|(&face, face_coordinates)| {
221 let face_nodes = &block_faces_nodes[face];
222 if face_nodes.contains(&node) {
223 let num_nodes_face = face_coordinates.len() as Scalar;
224 let face_center = face_coordinates
225 .iter()
226 .cloned()
227 .sum::<ReferenceCoordinate>()
228 / num_nodes_face;
229 let mut face_coordinates_one_ahead =
230 VecDeque::from(face_coordinates.clone());
231 let first_entry = face_coordinates_one_ahead.pop_front().unwrap();
232 face_coordinates_one_ahead.push_back(first_entry);
233 Some(
234 face_coordinates
235 .into_iter()
236 .zip(face_coordinates_one_ahead)
237 .zip(face_nodes.iter())
238 .map(
239 |(
240 (node_a_coordinates, node_b_coordinates),
241 &node_a,
242 )| {
243 let node_a_spot = face_nodes
244 .iter()
245 .position(|&n| n == node_a)
246 .unwrap();
247 let node_b = if node_a_spot + 1 == face_nodes.len()
248 {
249 face_nodes[0]
250 } else {
251 face_nodes[node_a_spot + 1]
252 };
253 let factor = if node == node_a || node == node_b {
254 1.0 + 1.0 / num_nodes_face
255 } else {
256 1.0 / num_nodes_face
257 };
258 let e_1 = &node_b_coordinates - node_a_coordinates;
259 let e_2 = &face_center - node_b_coordinates;
260 e_1.cross(e_2) * factor
261 },
262 )
263 .sum::<ReferenceCoordinate>(),
264 )
265 } else {
266 None
267 }
268 })
269 .sum::<ReferenceCoordinate>()
270 / (element_volume * 6.0)
271 })
272 .collect(),
273 ]
274 .into();
275 Self {
276 faces_nodes,
277 gradient_vectors,
278 integration_weights,
279 stabilization: 0.1,
280 tetrahedra,
281 tetrahedra_nodes,
282 }
283 }
284}
285
286impl Debug for Element {
287 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
288 write!(f, "VirtualElement {{ ... }}",)
289 }
290}
291
292pub enum VirtualElementError {
293 Upstream(String, String),
294}
295
296impl From<VirtualElementError> for AssertionError {
297 fn from(error: VirtualElementError) -> Self {
298 Self {
299 message: error.to_string(),
300 }
301 }
302}
303
304impl StyledError for VirtualElementError {
305 fn message(&self, style: &Style) -> String {
306 let c = style.frame;
307 match self {
308 Self::Upstream(error, element) => format!(
309 "{error}{c}\n\
310 In virtual element: {element}."
311 ),
312 }
313 }
314}
315
316styled_error!(VirtualElementError);
317
318#[test]
319fn temporary_poly_0() {
320 use crate::vem::NodalReferenceCoordinates;
321 let phi = (1.0 + 5.0_f64.sqrt()) / 2.0;
322 let coordinates = NodalReferenceCoordinates::from(vec![
323 [-1.0, -1.0, -1.0],
324 [-1.0, -1.0, 1.0],
325 [-1.0, 1.0, -1.0],
326 [-1.0, 1.0, 1.0],
327 [1.0, -1.0, -1.0],
328 [1.0, -1.0, 1.0],
329 [1.0, 1.0, -1.0],
330 [1.0, 1.0, 1.0],
331 [0.0, -phi, -1.0 / phi],
332 [0.0, -phi, 1.0 / phi],
333 [0.0, phi, -1.0 / phi],
334 [0.0, phi, 1.0 / phi],
335 [-phi, -1.0 / phi, 0.0],
336 [-phi, 1.0 / phi, 0.0],
337 [phi, -1.0 / phi, 0.0],
338 [phi, 1.0 / phi, 0.0],
339 [-1.0 / phi, 0.0, -phi],
340 [1.0 / phi, 0.0, -phi],
341 [-1.0 / phi, 0.0, phi],
342 [1.0 / phi, 0.0, phi],
343 ]);
344 let face_node_connectivity = vec![
345 vec![16, 17, 4, 8, 0],
346 vec![12, 13, 2, 16, 0],
347 vec![8, 9, 1, 12, 0],
348 vec![9, 5, 19, 18, 1],
349 vec![18, 3, 13, 12, 1],
350 vec![10, 6, 17, 16, 2],
351 vec![13, 3, 11, 10, 2],
352 vec![7, 11, 3, 18, 19],
353 vec![14, 5, 9, 8, 4],
354 vec![6, 15, 14, 4, 17],
355 vec![5, 14, 15, 7, 19],
356 vec![6, 10, 11, 7, 15],
357 ];
358 let element_face_connectivity = vec![vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]];
359 use crate::constitutive::solid::hyperelastic::NeoHookean;
360 use crate::fem::solid::elastic::ElasticElements;
361 use crate::vem::block::{Block, solid::SolidVirtualElements};
362 let block = Block::<_, Element>::from((
363 NeoHookean {
364 shear_modulus: 3.0,
365 bulk_modulus: 13.0,
366 },
367 element_face_connectivity.clone(),
368 face_node_connectivity.clone(),
369 &coordinates,
370 ));
371 use crate::fem::solid::NodalForcesSolid;
372 use crate::math::TensorArray;
373 use crate::mechanics::DeformationGradient;
374 use crate::vem::NodalCoordinates;
375 let coordinates_current = NodalCoordinates::from(coordinates.clone());
376 Assert::default()
377 .eq_within_tols(
378 DeformationGradient::identity(),
379 &block.deformation_gradients(&coordinates_current)[0][0],
380 )
381 .unwrap();
382 Assert::default()
383 .eq_within_tols(
384 NodalForcesSolid::zero(coordinates_current.len()),
385 &block.nodal_forces(&coordinates_current).unwrap(),
386 )
387 .unwrap();
388 let length = (coordinates[face_node_connectivity[0][0]].clone()
389 - coordinates[face_node_connectivity[0][1]].clone())
390 .norm();
391 let volume = (15.0 + 7.0 * 5.0_f64.sqrt()) / 4.0 * length.powi(3);
392 assert!((block.elements()[0].integration_weights()[0] / volume - 1.0).abs() < 1e-14);
393}
394
395#[test]
396fn temporary_poly_1() {
397 use crate::vem::NodalReferenceCoordinates;
398 let coordinates = NodalReferenceCoordinates::from(vec![
399 [-0.7727027, -0.65398245, -0.80050964],
400 [-0.55585269, -1.31907453, 1.32652506],
401 [-0.68068751, 0.86362469, -0.58348725],
402 [-1.2475506, 1.06566759, 1.45034587],
403 [1.47277602, -1.10640079, -0.90724596],
404 [1.10274756, -0.69153902, 1.27617253],
405 [0.64323505, 1.36639746, -1.48447683],
406 [0.91277928, 0.97322043, 0.67055],
407 [-0.19978796, -2.0201241, -0.50145446],
408 [-0.07547771, -1.54630032, 0.22127876],
409 [0.37534904, 1.50203587, -0.81372091],
410 [-0.20273152, 1.4672534, 0.27738481],
411 [-1.98854772, -0.25595864, 0.16143842],
412 [-1.80085125, 0.19913772, -0.19452172],
413 [1.3154974, -0.72436122, 0.17437191],
414 [2.09624968, 1.01585944, 0.29687302],
415 [-0.61664715, 0.18078644, -1.94806432],
416 [0.86740811, -0.38259605, -1.2754194],
417 [-1.08169702, -0.39837623, 1.63255916],
418 [0.12293689, -0.48172557, 1.4158596],
419 ]);
420 let face_node_connectivity = vec![
421 vec![16, 17, 4, 8, 0],
422 vec![12, 13, 2, 16, 0],
423 vec![8, 9, 1, 12, 0],
424 vec![9, 5, 19, 18, 1],
425 vec![18, 3, 13, 12, 1],
426 vec![10, 6, 17, 16, 2],
427 vec![13, 3, 11, 10, 2],
428 vec![7, 11, 3, 18, 19],
429 vec![14, 5, 9, 8, 4],
430 vec![6, 15, 14, 4, 17],
431 vec![5, 14, 15, 7, 19],
432 vec![6, 10, 11, 7, 15],
433 ];
434 let element_face_connectivity = vec![vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]];
435 use crate::constitutive::solid::hyperelastic::NeoHookean;
436 use crate::fem::solid::elastic::ElasticElements;
437 use crate::vem::block::{Block, solid::SolidVirtualElements};
438 let block = Block::<_, Element>::from((
439 NeoHookean {
440 shear_modulus: 3.0,
441 bulk_modulus: 13.0,
442 },
443 element_face_connectivity.clone(),
444 face_node_connectivity.clone(),
445 &coordinates,
446 ));
447 use crate::fem::solid::NodalForcesSolid;
448 use crate::math::TensorArray;
449 use crate::mechanics::DeformationGradient;
450 use crate::vem::NodalCoordinates;
451 let coordinates_current = NodalCoordinates::from(coordinates.clone());
452 Assert::default()
453 .eq_within_tols(
454 DeformationGradient::identity(),
455 &block.deformation_gradients(&coordinates_current)[0][0],
456 )
457 .unwrap();
458 Assert::default()
459 .eq_within_tols(
460 NodalForcesSolid::zero(coordinates_current.len()),
461 &block.nodal_forces(&coordinates_current).unwrap(),
462 )
463 .unwrap();
464 use crate::mechanics::test::{get_deformation_gradient, get_translation_current_configuration};
465 let coordinates_current: NodalCoordinates = coordinates
466 .iter()
467 .map(|coord| get_deformation_gradient() * coord + get_translation_current_configuration())
468 .collect();
469 Assert::default()
470 .eq_within_tols(
471 get_deformation_gradient(),
472 &block.deformation_gradients(&coordinates_current)[0][0],
473 )
474 .unwrap();
475}
476
477#[test]
478fn temporary_poly_2() {
479 use crate::vem::NodalReferenceCoordinates;
480 let phi = (1.0 + 5.0_f64.sqrt()) / 2.0;
481 let coordinates_0 = NodalReferenceCoordinates::from(vec![
482 [-1.0, -1.0, -1.0],
483 [-1.0, -1.0, 1.0],
484 [-1.0, 1.0, -1.0],
485 [-1.0, 1.0, 1.0],
486 [1.0, -1.0, -1.0],
487 [1.0, -1.0, 1.0],
488 [1.0, 1.0, -1.0],
489 [1.0, 1.0, 1.0],
490 [0.0, -phi, -1.0 / phi],
491 [0.0, -phi, 1.0 / phi],
492 [0.0, phi, -1.0 / phi],
493 [0.0, phi, 1.0 / phi],
494 [-phi, -1.0 / phi, 0.0],
495 [-phi, 1.0 / phi, 0.0],
496 [phi, -1.0 / phi, 0.0],
497 [phi, 1.0 / phi, 0.0],
498 [-1.0 / phi, 0.0, -phi],
499 [1.0 / phi, 0.0, -phi],
500 [-1.0 / phi, 0.0, phi],
501 [1.0 / phi, 0.0, phi],
502 ]);
503 let face_node_connectivity = vec![
504 vec![16, 17, 4, 8, 0],
505 vec![12, 13, 2, 16, 0],
506 vec![8, 9, 1, 12, 0],
507 vec![9, 5, 19, 18, 1],
508 vec![18, 3, 13, 12, 1],
509 vec![10, 6, 17, 16, 2],
510 vec![13, 3, 11, 10, 2],
511 vec![7, 11, 3, 18, 19],
512 vec![14, 5, 9, 8, 4],
513 vec![6, 15, 14, 4, 17],
514 vec![5, 14, 15, 7, 19],
515 vec![6, 10, 11, 7, 15],
516 ];
517 let element_face_connectivity = vec![vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]];
518 use crate::constitutive::solid::hyperelastic::NeoHookean;
519 use crate::fem::solid::elastic::ElasticElements;
520 use crate::vem::block::Block;
521 let block = Block::<_, Element>::from((
522 NeoHookean {
523 shear_modulus: 3.0,
524 bulk_modulus: 13.0,
525 },
526 element_face_connectivity.clone(),
527 face_node_connectivity.clone(),
528 &coordinates_0,
529 ));
530 use crate::vem::NodalCoordinates;
531 let coordinates = NodalCoordinates::from(vec![
532 [-0.7727027, -0.65398245, -0.80050964],
533 [-0.55585269, -1.31907453, 1.32652506],
534 [-0.68068751, 0.86362469, -0.58348725],
535 [-1.2475506, 1.06566759, 1.45034587],
536 [1.47277602, -1.10640079, -0.90724596],
537 [1.10274756, -0.69153902, 1.27617253],
538 [0.64323505, 1.36639746, -1.48447683],
539 [0.91277928, 0.97322043, 0.67055],
540 [-0.19978796, -2.0201241, -0.50145446],
541 [-0.07547771, -1.54630032, 0.22127876],
542 [0.37534904, 1.50203587, -0.81372091],
543 [-0.20273152, 1.4672534, 0.27738481],
544 [-1.98854772, -0.25595864, 0.16143842],
545 [-1.80085125, 0.19913772, -0.19452172],
546 [1.3154974, -0.72436122, 0.17437191],
547 [2.09624968, 1.01585944, 0.29687302],
548 [-0.61664715, 0.18078644, -1.94806432],
549 [0.86740811, -0.38259605, -1.2754194],
550 [-1.08169702, -0.39837623, 1.63255916],
551 [0.12293689, -0.48172557, 1.4158596],
552 ]);
553 use crate::EPSILON;
554 use crate::fem::solid::hyperelastic::HyperelasticElements;
555 let mut finite_difference = 0.0;
556 let nodal_forces_fd = (0..coordinates.len())
557 .map(|node| {
558 (0..3)
559 .map(|i| {
560 let mut nodal_coordinates = coordinates.clone();
561 nodal_coordinates[node][i] += 0.5 * EPSILON;
562 finite_difference = block.helmholtz_free_energy(&nodal_coordinates).unwrap();
563 nodal_coordinates[node][i] -= EPSILON;
564 finite_difference -= block.helmholtz_free_energy(&nodal_coordinates).unwrap();
565 finite_difference / EPSILON
566 })
567 .collect()
568 })
569 .collect();
570 Assert::default()
571 .eq_within_fd_tol(block.nodal_forces(&coordinates).unwrap(), &nodal_forces_fd)
572 .unwrap();
573 let mut finite_difference = 0.0;
574 let nodal_stiffnesses_fd = (0..coordinates.len())
575 .map(|a| {
576 (0..coordinates.len())
577 .map(|b| {
578 (0..3)
579 .map(|i| {
580 (0..3)
581 .map(|j| {
582 let mut nodal_coordinates = coordinates.clone();
583 nodal_coordinates[b][j] += 0.5 * EPSILON;
584 finite_difference =
585 block.nodal_forces(&nodal_coordinates).unwrap()[a][i];
586 nodal_coordinates[b][j] -= EPSILON;
587 finite_difference -=
588 block.nodal_forces(&nodal_coordinates).unwrap()[a][i];
589 finite_difference / EPSILON
590 })
591 .collect()
592 })
593 .collect()
594 })
595 .collect()
596 })
597 .collect();
598 Assert::default()
599 .eq_within_fd_tol(
600 block.nodal_stiffnesses(&coordinates).unwrap(),
601 &nodal_stiffnesses_fd,
602 )
603 .unwrap();
604}