Skip to main content

conspire/geometry/ntree/into/
mod.rs

1use crate::geometry::ntree::node::cell::Cell;
2use crate::{
3    geometry::{Coordinates, grid::Grid, mesh::Connectivity, ntree::Orthotree},
4    math::{Tensor, TensorVec},
5};
6use std::{array::from_fn, collections::HashMap};
7
8impl<const D: usize, const L: usize, const M: usize, const N: usize, U>
9    From<Orthotree<D, L, M, N, u16, U>> for (Vec<[usize; N]>, Coordinates<D>)
10{
11    fn from(orthotree: Orthotree<D, L, M, N, u16, U>) -> Self {
12        let mut coord_map: HashMap<u64, usize> = HashMap::new();
13        let mut coords = Coordinates::<D>::new();
14        let face_mask: usize = if D <= 2 { (1 << D) - 1 } else { 3 };
15        let connectivity: Vec<[usize; N]> = orthotree
16            .nodes
17            .iter()
18            .filter(|node| node.is_leaf())
19            .map(|node| {
20                from_fn(|i| {
21                    let face = i & face_mask;
22                    let vertex_i = (i & !face_mask) | (face ^ (face >> 1));
23                    let vertex: [u16; D] = from_fn(|ax| {
24                        if (vertex_i >> ax) & 1 == 1 {
25                            node.corner[ax] + node.length
26                        } else {
27                            node.corner[ax]
28                        }
29                    });
30                    let key: u64 =
31                        (0..D).fold(0u64, |acc, ax| acc | ((vertex[ax] as u64) << (16 * ax)));
32                    if let Some(&idx) = coord_map.get(&key) {
33                        idx
34                    } else {
35                        let idx = coords.len();
36                        coords.push(from_fn(|ax| vertex[ax] as f64).into());
37                        coord_map.insert(key, idx);
38                        idx
39                    }
40                })
41            })
42            .collect();
43        (connectivity, coords)
44    }
45}
46
47impl<const L: usize, const M: usize, U> From<Orthotree<2, L, M, 4, u16, U>>
48    for (Connectivity, Coordinates<2>)
49{
50    fn from(orthotree: Orthotree<2, L, M, 4, u16, U>) -> Self {
51        let (connectivity, coordinates): (Vec<[usize; 4]>, _) = orthotree.into();
52        (
53            Connectivity::Quadrilateral(connectivity.into()),
54            coordinates,
55        )
56    }
57}
58
59impl<const L: usize, const M: usize, U> From<Orthotree<3, L, M, 8, u16, U>>
60    for (Connectivity, Coordinates<3>)
61{
62    fn from(orthotree: Orthotree<3, L, M, 8, u16, U>) -> Self {
63        let (connectivity, coordinates): (Vec<[usize; 8]>, _) = orthotree.into();
64        (Connectivity::Hexahedral(connectivity.into()), coordinates)
65    }
66}
67
68impl<const D: usize, const L: usize, const M: usize, const N: usize, T, U, V>
69    From<&Orthotree<D, L, M, N, T, U, V>> for Grid<D, V>
70where
71    T: Cell,
72    V: Copy,
73{
74    fn from(orthotree: &Orthotree<D, L, M, N, T, U, V>) -> Self {
75        let leaves: Vec<([usize; D], usize, V)> = orthotree
76            .nodes
77            .iter()
78            .filter_map(|node| {
79                node.value.map(|value| {
80                    (
81                        from_fn(|ax| node.corner[ax].cells()),
82                        node.length.cells(),
83                        value,
84                    )
85                })
86            })
87            .collect();
88        let nel: [usize; D] = from_fn(|ax| {
89            leaves
90                .iter()
91                .map(|(corner, length, _)| corner[ax] + length)
92                .max()
93                .unwrap_or(0)
94        });
95        let count: usize = nel.iter().product();
96        if count == 0 {
97            return Grid::new(Vec::new(), nel);
98        }
99        let mut data = vec![leaves[0].2; count];
100        for (corner, length, value) in leaves {
101            let extent: [usize; D] = from_fn(|ax| (corner[ax] + length).min(nel[ax]) - corner[ax]);
102            for cell in 0..extent.iter().product() {
103                let mut rem = cell;
104                let mut flat = 0;
105                let mut stride = 1;
106                for ax in 0..D {
107                    flat += (corner[ax] + rem % extent[ax]) * stride;
108                    rem /= extent[ax];
109                    stride *= nel[ax];
110                }
111                data[flat] = value;
112            }
113        }
114        Grid::new(data, nel)
115    }
116}