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

1use crate::math::Reference;
2pub mod solid;
3#[cfg(test)]
4mod test;
5
6use crate::{
7    domain::block::element::{ElementError, ElementKind},
8    fem::block::element::{
9        ElementNodalReferenceCoordinates as FemElementNodalReferenceCoordinates, FiniteElement,
10        linear::Tetrahedron,
11    },
12    math::{
13        CrossProduct, Current, Quantity, Scalar, Tensor, TensorArray, TensorRank1, TensorRank1List,
14        TensorRank1Vec, TensorRank1Vec2D, TensorVector,
15    },
16    mechanics::ReferenceCoordinate,
17    units::{Area, Length, ReciprocalLength, Velocity, Volume},
18    vem::{NodalCoordinates, NodalReferenceCoordinates, NodalVelocities},
19};
20
21use std::fmt::{self, Debug, Display, Formatter};
22
23pub type ElementNodalCoordinates = NodalCoordinates;
24pub type ElementNodalVelocities = NodalVelocities;
25pub type ElementNodalReferenceCoordinates = TensorRank1Vec2D<3, Reference, Length>;
26pub type GradientVectors = TensorRank1Vec2D<3, Reference, ReciprocalLength>;
27pub type IntegrationWeights = TensorVector<Quantity<Volume>>;
28
29pub type TetrahedraQuantities<U> = Vec<TensorRank1List<3, Current, 4, U>>;
30pub type TetrahedraCoordinates = TetrahedraQuantities<Length>;
31pub type TetrahedraVelocities = TetrahedraQuantities<Velocity>;
32
33pub struct Element {
34    faces_nodes: Vec<Vec<usize>>,
35    gradient_vectors: GradientVectors,
36    integration_weights: IntegrationWeights,
37    stabilization: Scalar,
38    tetrahedra: Vec<Tetrahedron>,
39    tetrahedra_nodes: Vec<[usize; 3]>,
40}
41
42impl Element {
43    pub(crate) fn upstream(&self, error: impl Display) -> VirtualElementError {
44        VirtualElementError::upstream(error, self)
45    }
46}
47
48pub trait VirtualElement
49where
50    for<'a> Self: From<(
51        ElementNodalReferenceCoordinates,
52        &'a [usize],
53        &'a [usize],
54        &'a [Vec<usize>],
55    )>,
56{
57    fn element_center<U>(
58        nodal_quantities: &TensorRank1Vec<3, Current, U>,
59    ) -> TensorRank1<3, Current, U>;
60    fn faces_centers<U>(
61        &self,
62        nodal_quantities: &TensorRank1Vec<3, Current, U>,
63    ) -> TensorRank1Vec<3, Current, U>;
64    fn faces_nodes(&self) -> &[Vec<usize>];
65    fn gradient_vectors(&self) -> &GradientVectors;
66    fn integration_weights(&self) -> &IntegrationWeights;
67    fn stabilization(&self) -> Scalar;
68    fn tetrahedra(&self) -> &[Tetrahedron];
69    fn tetrahedra_coordinates<U>(
70        &self,
71        nodal_quantities: &TensorRank1Vec<3, Current, U>,
72    ) -> TetrahedraQuantities<U>;
73    fn tetrahedra_nodes(&self) -> &[[usize; 3]];
74}
75
76impl VirtualElement for Element {
77    fn element_center<U>(
78        nodal_quantities: &TensorRank1Vec<3, Current, U>,
79    ) -> TensorRank1<3, Current, U> {
80        nodal_quantities
81            .iter()
82            .cloned()
83            .sum::<TensorRank1<3, Current, U>>()
84            / nodal_quantities.len() as Scalar
85    }
86    fn faces_centers<U>(
87        &self,
88        nodal_quantities: &TensorRank1Vec<3, Current, U>,
89    ) -> TensorRank1Vec<3, Current, U> {
90        self.faces_nodes()
91            .iter()
92            .map(|face_nodes| {
93                face_nodes
94                    .iter()
95                    .map(|&face_node| nodal_quantities[face_node].clone())
96                    .sum::<TensorRank1<3, Current, U>>()
97                    / (face_nodes.len() as Scalar)
98            })
99            .collect()
100    }
101    fn faces_nodes(&self) -> &[Vec<usize>] {
102        &self.faces_nodes
103    }
104    fn gradient_vectors(&self) -> &GradientVectors {
105        &self.gradient_vectors
106    }
107    fn integration_weights(&self) -> &IntegrationWeights {
108        &self.integration_weights
109    }
110    fn stabilization(&self) -> Scalar {
111        self.stabilization
112    }
113    fn tetrahedra(&self) -> &[Tetrahedron] {
114        &self.tetrahedra
115    }
116    fn tetrahedra_coordinates<U>(
117        &self,
118        nodal_quantities: &TensorRank1Vec<3, Current, U>,
119    ) -> TetrahedraQuantities<U> {
120        let element_center = Self::element_center(nodal_quantities);
121        let faces_centers = self.faces_centers(nodal_quantities);
122        self.tetrahedra_nodes()
123            .iter()
124            .map(|&[face, node_b, node_a]| {
125                [
126                    faces_centers[face].clone(),
127                    nodal_quantities[node_b].clone(),
128                    nodal_quantities[node_a].clone(),
129                    element_center.clone(),
130                ]
131                .into()
132            })
133            .collect()
134    }
135    fn tetrahedra_nodes(&self) -> &[[usize; 3]] {
136        &self.tetrahedra_nodes
137    }
138}
139
140impl
141    From<(
142        ElementNodalReferenceCoordinates,
143        &[usize],
144        &[usize],
145        &[Vec<usize>],
146    )> for Element
147{
148    fn from(
149        (reference_nodal_coordinates, element_faces, element_nodes, block_faces_nodes): (
150            ElementNodalReferenceCoordinates,
151            &[usize],
152            &[usize],
153            &[Vec<usize>],
154        ),
155    ) -> Self {
156        let faces_nodes = element_faces
157            .iter()
158            .map(|&element_face| {
159                block_faces_nodes[element_face]
160                    .iter()
161                    .map(|face_node| {
162                        element_nodes
163                            .iter()
164                            .position(|element_node| face_node == element_node)
165                            .unwrap()
166                    })
167                    .collect::<Vec<_>>()
168            })
169            .collect::<Vec<_>>();
170        let mut nodal_coordinates =
171            NodalReferenceCoordinates::from(vec![
172                ReferenceCoordinate::from([0.0, 0.0, 0.0]);
173                element_nodes.len()
174            ]);
175        faces_nodes
176            .iter()
177            .zip(reference_nodal_coordinates.iter())
178            .for_each(|(face_nodes, face_coordinates)| {
179                face_nodes
180                    .iter()
181                    .zip(face_coordinates.iter())
182                    .for_each(|(&node, coordinates)| nodal_coordinates[node] = coordinates.clone())
183            });
184        let element_center = nodal_coordinates.into_iter().sum::<ReferenceCoordinate>()
185            / (element_nodes.len() as Scalar);
186        let mut area_vectors = vec![TensorRank1::<3, Reference, Area>::zero(); element_nodes.len()];
187        let tetrahedra_nodes = faces_nodes
188            .iter()
189            .enumerate()
190            .flat_map(|(face, face_nodes)| {
191                (0..face_nodes.len())
192                    .map(|spot| {
193                        [
194                            face,
195                            face_nodes[(spot + 1) % face_nodes.len()],
196                            face_nodes[spot],
197                        ]
198                    })
199                    .collect::<Vec<_>>()
200            })
201            .collect::<Vec<_>>();
202        let tetrahedra = faces_nodes
203            .iter()
204            .zip(reference_nodal_coordinates.iter())
205            .flat_map(|(face_nodes, face_coordinates)| {
206                let num_nodes_face = face_coordinates.len();
207                let face_center = face_coordinates
208                    .iter()
209                    .cloned()
210                    .sum::<ReferenceCoordinate>()
211                    / (num_nodes_face as Scalar);
212                let mut face_area_vector = TensorRank1::<3, Reference, Area>::zero();
213                let face_tetrahedra = (0..num_nodes_face)
214                    .map(|spot| {
215                        let next = (spot + 1) % num_nodes_face;
216                        let e_1 = &face_coordinates[next] - &face_coordinates[spot];
217                        let e_2 = &face_center - &face_coordinates[next];
218                        let cross = e_1.cross(&e_2);
219                        face_area_vector += &cross;
220                        area_vectors[face_nodes[spot]] += &cross;
221                        area_vectors[face_nodes[next]] += &cross;
222                        Tetrahedron::from(FemElementNodalReferenceCoordinates::from([
223                            face_center.clone(),
224                            face_coordinates[next].clone(),
225                            face_coordinates[spot].clone(),
226                            element_center.clone(),
227                        ]))
228                    })
229                    .collect::<Vec<_>>();
230                let shared = &face_area_vector / (num_nodes_face as Scalar);
231                face_nodes
232                    .iter()
233                    .for_each(|&node| area_vectors[node] += &shared);
234                face_tetrahedra
235            })
236            .collect::<Vec<_>>();
237        let element_volume = tetrahedra
238            .iter()
239            .map(|tetrahedron| tetrahedron.volume())
240            .sum::<Quantity<Volume>>();
241        let gradient_vectors = GradientVectors::from(vec![
242            area_vectors
243                .into_iter()
244                .map(|area_vector| area_vector / (element_volume * 6.0))
245                .collect::<TensorRank1Vec<3, Reference, ReciprocalLength>>(),
246        ]);
247        let integration_weights = IntegrationWeights::from([element_volume]);
248        Self {
249            faces_nodes,
250            gradient_vectors,
251            integration_weights,
252            stabilization: 0.1,
253            tetrahedra,
254            tetrahedra_nodes,
255        }
256    }
257}
258
259impl Debug for Element {
260    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
261        write!(f, "VirtualElement {{ ... }}",)
262    }
263}
264
265pub struct VirtualElementKind;
266
267impl ElementKind for VirtualElementKind {
268    const NAME: &'static str = "virtual element";
269}
270
271pub type VirtualElementError = ElementError<VirtualElementKind>;
conspire/domain/vem/block/element/mod.rs

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