Skip to main content

conspire/geometry/mesh/from/ntree/polyhedra/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::geometry::{
5    Coordinates,
6    mesh::{
7        Connectivity, Mesh,
8        from::ntree::facets::{Facet, Facets, corner_length, facet, leaves},
9    },
10    ntree::{
11        Octree, Orthotree, Quadtree,
12        node::{cell::Cell, slot::Slot},
13    },
14};
15
16impl<T, U, V> From<Octree<T, U, V>> for Mesh<3>
17where
18    T: Cell,
19    U: Slot,
20{
21    fn from(octree: Octree<T, U, V>) -> Self {
22        let (elements_faces, faces_nodes, mut coordinates) = polytopes(&octree);
23        octree.rescale_coordinates(&mut coordinates);
24        (
25            vec![Connectivity::Polyhedral(
26                (elements_faces, faces_nodes).into(),
27            )],
28            coordinates,
29        )
30            .into()
31    }
32}
33
34impl<T, U, V> From<Quadtree<T, U, V>> for Mesh<2>
35where
36    T: Cell,
37    U: Slot,
38{
39    fn from(quadtree: Quadtree<T, U, V>) -> Self {
40        let (elements_faces, faces_nodes, mut coordinates) = polytopes(&quadtree);
41        quadtree.rescale_coordinates(&mut coordinates);
42        (
43            vec![Connectivity::Polygonal(
44                (elements_faces, faces_nodes).into(),
45            )],
46            coordinates,
47        )
48            .into()
49    }
50}
51
52fn polytopes<const D: usize, const L: usize, const M: usize, const N: usize, T, U, V>(
53    tree: &Orthotree<D, L, M, N, T, U, V>,
54) -> (Vec<Vec<usize>>, Vec<Vec<usize>>, Coordinates<D>)
55where
56    T: Cell,
57    U: Slot,
58{
59    let (leaves, element_of) = leaves(tree);
60    let facets = Facets::<D>::new::<L, M, N, T, U, V>(tree, &leaves);
61    let mut elements_faces = vec![Vec::new(); leaves.len()];
62    let mut faces_nodes = Vec::<Vec<usize>>::new();
63    let mut emit = |corner: [usize; D],
64                    size: usize,
65                    axis: usize,
66                    plane: usize,
67                    elements: &[usize],
68                    flip: bool| {
69        let index = faces_nodes.len();
70        faces_nodes.push(facets.polygon(corner, size, axis, plane, flip));
71        elements
72            .iter()
73            .for_each(|&element| elements_faces[element].push(index));
74    };
75    for &index in &leaves {
76        let element = element_of[index];
77        let (corner, length) = corner_length(&tree.nodes[index]);
78        for f in 0..M {
79            let (axis, side) = (f >> 1, f & 1);
80            let plane = corner[axis] + side * length;
81            match facet(tree, index, f) {
82                Facet::Refined(fine) => fine.into_iter().for_each(|leaf| {
83                    let (fine_corner, fine_length) = corner_length(&tree.nodes[leaf]);
84                    let fine_element = element_of[leaf];
85                    let flip = if fine_element < element {
86                        side == 1
87                    } else {
88                        side == 0
89                    };
90                    emit(
91                        fine_corner,
92                        fine_length,
93                        axis,
94                        plane,
95                        &[fine_element, element],
96                        flip,
97                    )
98                }),
99                Facet::Neighbor(neighbor) => {
100                    if side == 0 {
101                        let neighbor_element = element_of[neighbor];
102                        emit(
103                            corner,
104                            length,
105                            axis,
106                            plane,
107                            &[neighbor_element, element],
108                            element < neighbor_element,
109                        )
110                    }
111                }
112                Facet::Boundary => emit(corner, length, axis, plane, &[element], side == 0),
113                Facet::Absent => {}
114            }
115        }
116    }
117    (elements_faces, faces_nodes, facets.into_coordinates())
118}