conspire/geometry/ntree/from/tessellation/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    geometry::{
6        Coordinate, CoordinateList,
7        bbox::BoundingBox,
8        mesh::Tessellation,
9        ntree::{
10            Octree,
11            balance::Balancing,
12            node::{Kind, Node, cell::Cell, slot::Slot},
13            pair::Pairing,
14            rescale::Rescaling,
15            sizing::{Sizing, curvature::CurvatureSizing},
16        },
17    },
18    math::{Quantity, Scalar},
19};
20use std::array::from_fn;
21
22const D: usize = 3;
23const M: usize = 6;
24
25impl<T, U> Octree<T, U>
26where
27    T: Cell,
28    U: Slot,
29{
30    /// Builds an octree from a tessellation, refining cells where either the
31    /// local thickness or the local curvature demands a smaller size.
32    ///
33    /// `scale` controls cells-per-thickness; `curvature` controls
34    /// curvature-driven refinement independent of thickness (e.g. a sphere
35    /// has ~constant thickness everywhere but can still demand refinement
36    /// from curvature alone). `padding` adds extra empty root levels in case
37    /// the tessellation's boundary overlaps the primordial primal node.
38    pub fn from_features(
39        tessellation: &Tessellation,
40        scale: Scalar,
41        curvature: CurvatureSizing,
42        padding: u16,
43    ) -> Result<Self, &'static str> {
44        Self::refine(&Sizing::new(tessellation, scale, curvature, padding))
45    }
46    /// Refines an octree to a given size field.
47    pub fn refine(sizing: &Sizing) -> Result<Self, &'static str> {
48        let Sizing {
49            center,
50            coordinates,
51            elements,
52            levels,
53            min_length,
54            scale,
55            targets,
56        } = sizing;
57        let (center, min_length, scale) = (center, *min_length, *scale);
58        if elements.is_empty() {
59            return Ok(Self {
60                balanced: Balancing::None,
61                nodes: vec![Node {
62                    corner: from_fn(|_| T::ZERO),
63                    length: T::ONE,
64                    facets: [None; M],
65                    kind: Kind::Leaf,
66                    value: None,
67                }],
68                paired: Pairing::None,
69                rescale: Rescaling {
70                    center: Coordinate::const_from([0.0; D]),
71                    cell: Quantity::new(1.0),
72                    half: 0.0,
73                },
74            });
75        }
76        let root_length = 1usize
77            .checked_shl(*levels)
78            .and_then(T::length)
79            .ok_or("sizing field exceeds maximum octree depth")?;
80        let half = root_length.scalar() / 2.0;
81        let mut tree = Self {
82            balanced: Balancing::None,
83            rescale: Rescaling {
84                center: center.clone(),
85                cell: min_length,
86                half,
87            },
88            nodes: vec![Node {
89                corner: from_fn(|_| T::ZERO),
90                length: root_length,
91                facets: [None; M],
92                kind: Kind::Leaf,
93                value: None,
94            }],
95            paired: Pairing::None,
96        };
97        let overlaps = |bbox: &BoundingBox<3>, triangle: usize| {
98            let element = elements[triangle];
99            bbox.overlaps_triangle(
100                &coordinates[element[0]],
101                &coordinates[element[1]],
102                &coordinates[element[2]],
103            )
104        };
105        let mut stack: Vec<(usize, Vec<usize>)> = vec![(0, (0..elements.len()).collect())];
106        while let Some((index, overlapping)) = stack.pop() {
107            let cells: usize = tree.nodes[index].length.cells();
108            let extent: Scalar = tree.nodes[index].length.scalar();
109            let target = overlapping
110                .iter()
111                .map(|&triangle| targets[triangle])
112                .fold(Quantity::new(Scalar::INFINITY), Quantity::min);
113            if min_length * (extent * scale) <= target {
114                continue;
115            }
116            if cells <= 1 {
117                continue;
118            }
119            tree.subdivide(index)?;
120            let children: Vec<usize> = tree.nodes[index]
121                .orthants()
122                .unwrap()
123                .iter()
124                .map(|&child| child.slot())
125                .collect();
126            for child in children {
127                let corner = tree.nodes[child].corner;
128                let child_extent: Scalar = tree.nodes[child].length.scalar();
129                let minimum = Coordinate::<3>::from(from_fn::<_, 3, _>(|ax| {
130                    center[ax] + min_length * (corner[ax].scalar() - half)
131                }));
132                let maximum = Coordinate::<3>::from(from_fn::<_, 3, _>(|ax| {
133                    minimum[ax] + min_length * child_extent
134                }));
135                let bbox = BoundingBox::from(CoordinateList::const_from([minimum, maximum]));
136                let inside: Vec<usize> = overlapping
137                    .iter()
138                    .copied()
139                    .filter(|&triangle| overlaps(&bbox, triangle))
140                    .collect();
141                if !inside.is_empty() {
142                    stack.push((child, inside));
143                }
144            }
145        }
146        Ok(tree)
147    }
148}
conspire/geometry/ntree/from/tessellation/mod.rs

conspire/geometry/ntree/from/tessellation/
mod.rs
