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