conspire/geometry/ntree/dual/
mod.rs

1pub mod octree;
2pub mod quadtree;
3
4use crate::{
5    geometry::{
6        Coordinate, Coordinates,
7        mesh::Mesh,
8        ntree::{
9            Orthotree,
10            balance::Balancing,
11            node::{Kind, split::Split},
12            pair::Pairing,
13        },
14    },
15    math::{Scalar, TensorVec},
16};
17use std::{array::from_fn, collections::HashMap, ops::Add};
18
19type NodeMap<const D: usize> = HashMap<[usize; D], usize>;
20
21fn get_or_add<const D: usize>(
22    coordinate: Coordinate<D>,
23    coordinates: &mut Coordinates<D>,
24    nodes_map: &mut NodeMap<D>,
25    node_index: &mut usize,
26) -> usize {
27    let key = from_fn(|i| (2.0 * coordinate[i]) as usize);
28    if let Some(&node) = nodes_map.get(&key) {
29        node
30    } else {
31        let node = *node_index;
32        coordinates.push(coordinate);
33        nodes_map.insert(key, node);
34        *node_index += 1;
35        node
36    }
37}
38
39pub trait Dualization<const D: usize> {
40    fn dualize(&mut self) -> Mesh<D>;
41}
42
43pub trait Star<const D: usize, const N: usize> {
44    fn star(&self, center_nodes: &[usize], connectivity: &mut Vec<[usize; N]>);
45}
46
47impl<const D: usize, const L: usize, const M: usize, const N: usize, T, U> Star<D, N>
48    for Orthotree<D, L, M, N, T, U>
49where
50    T: Add<Output = T> + Copy + PartialOrd + Split + Into<usize>,
51    U: Copy + Into<usize>,
52{
53    fn star(&self, center_nodes: &[usize], connectivity: &mut Vec<[usize; N]>) {
54        let face_mask: usize = if D <= 2 { (1 << D) - 1 } else { 3 };
55        let root = &self.nodes[0];
56        let lo = root.corner;
57        let hi: [T; D] = from_fn(|a| root.corner[a] + root.length);
58        for node in self.iter().filter(|node| node.is_leaf()) {
59            let vertex: [T; D] = from_fn(|a| node.corner[a] + node.length);
60            if (0..D).all(|a| lo[a] < vertex[a] && vertex[a] < hi[a]) {
61                let cells: [usize; N] = from_fn(|d| incident_leaf(self, &vertex, d));
62                let mut distinct = cells.to_vec();
63                distinct.sort_unstable();
64                distinct.dedup();
65                if distinct.len() != N {
66                    continue;
67                }
68                let lengths: [usize; N] = from_fn(|o| self.nodes[cells[o]].length.into());
69                let shortest = *lengths.iter().min().unwrap();
70                let longest = *lengths.iter().max().unwrap();
71                let coordinate: [usize; D] = from_fn(|a| vertex[a].into());
72                if longest == shortest || (0..D).all(|a| coordinate[a].is_multiple_of(2 * longest))
73                {
74                    connectivity.push(from_fn(|i| {
75                        let bits = i & face_mask;
76                        center_nodes[cells[(i & !face_mask) | (bits ^ (bits >> 1))]]
77                    }));
78                }
79            }
80        }
81    }
82}
83
84pub(crate) fn incident_leaf<const D: usize, const L: usize, const M: usize, const N: usize, T, U>(
85    tree: &Orthotree<D, L, M, N, T, U>,
86    vertex: &[T; D],
87    direction: usize,
88) -> usize
89where
90    T: Add<Output = T> + Copy + PartialOrd + Split,
91    U: Copy + Into<usize>,
92{
93    let mut index = 0;
94    loop {
95        match &tree.nodes[index].kind {
96            Kind::Leaf => return index,
97            Kind::Tree(orthants) => {
98                let corner = tree.nodes[index].corner;
99                let half = tree.nodes[index].length.split();
100                let child = (0..D).fold(0, |acc, a| {
101                    let mid = corner[a] + half;
102                    let bit = if vertex[a] > mid {
103                        1
104                    } else if vertex[a] < mid {
105                        0
106                    } else {
107                        (direction >> a) & 1
108                    };
109                    acc | (bit << a)
110                });
111                index = orthants[child].into();
112            }
113        }
114    }
115}
116
117pub trait Initialize<const D: usize, const N: usize> {
118    fn initialize(&self) -> (Vec<usize>, Coordinates<D>, usize, Vec<[usize; N]>);
119}
120
121impl<const D: usize, const L: usize, const M: usize, const N: usize, T, U> Initialize<D, N>
122    for Orthotree<D, L, M, N, T, U>
123where
124    T: Copy + Into<Scalar> + Into<usize>,
125    U: Copy + Into<usize>,
126{
127    fn initialize(&self) -> (Vec<usize>, Coordinates<D>, usize, Vec<[usize; N]>) {
128        assert!(!matches!(self.balanced, Balancing::None));
129        assert!(!matches!(self.paired, Pairing::None));
130        let num = self.len();
131        let mut center_nodes = vec![0; num];
132        let mut coordinates = Coordinates::with_capacity(num);
133        let mut node_index = 0;
134        self.iter()
135            .enumerate()
136            .filter(|(_, node)| node.is_leaf())
137            .for_each(|(index, leaf)| {
138                center_nodes[index] = node_index;
139                let length: Scalar = leaf.length.into();
140                let center = from_fn(|i| {
141                    let c: Scalar = leaf.corner[i].into();
142                    c + length * 0.5
143                });
144                coordinates.push(center.into());
145                node_index += 1;
146            });
147        (
148            center_nodes,
149            coordinates,
150            node_index,
151            Vec::with_capacity(num),
152        )
153    }
154}
conspire/geometry/ntree/dual/mod.rs

conspire/geometry/ntree/dual/
mod.rs
