Skip to main content

conspire/geometry/ntree/into/
mod.rs

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