conspire/geometry/mesh/tessellation/dual/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    geometry::{
6        Coordinate, Coordinates,
7        bvh::BoundingVolumeHierarchy,
8        mesh::{
9            Connectivity, Mesh,
10            tessellation::{D, Tessellation},
11        },
12        ntree::{Balance, Balancing, Dualization, Octree, Pairing},
13    },
14    math::{Scalar, Tensor, TensorVec},
15};
16use std::{
17    array::from_fn,
18    collections::{HashMap, hash_map::Entry},
19    thread::{available_parallelism, scope},
20};
21
22const GRAZING_TOLERANCE: Scalar = 1.0e-4;
23const TRIM_MARGIN: Scalar = 0.5;
24const EDGES: [[usize; 2]; 12] = [
25    [0, 1],
26    [1, 2],
27    [2, 3],
28    [3, 0],
29    [4, 5],
30    [5, 6],
31    [6, 7],
32    [7, 4],
33    [0, 4],
34    [1, 5],
35    [2, 6],
36    [3, 7],
37];
38const DIRECTIONS: [Coordinate<D>; 3] = [
39    Coordinate::const_from([1.0, 0.140_412_03, 0.092_153_88]),
40    Coordinate::const_from([0.097_153_2, 1.0, 0.131_771_4]),
41    Coordinate::const_from([0.123_456_7, 0.087_654_3, 1.0]),
42];
43impl Tessellation {
44    pub fn dualize(&self, balancing: Balancing, scale: Scalar) -> Result<Mesh<D>, &'static str> {
45        let mut octree = Octree::<u16, usize>::from_sdf(self, scale);
46        octree.equilibrate(balancing, Pairing::Regular)?;
47        let mut mesh = octree.dualize();
48        self.trim(&mut mesh, self.bvh());
49        self.buffer(mesh, self.bvh())
50    }
51    fn buffer(
52        &self,
53        mesh: Mesh<D>,
54        bvh: &BoundingVolumeHierarchy<D>,
55    ) -> Result<Mesh<D>, &'static str> {
56        let surface = self.mesh();
57        let surface_coordinates = surface.coordinates();
58        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
59        let boundary = mesh.exterior_faces();
60        let mut edges = HashMap::new();
61        boundary.iter().for_each(|face| {
62            (0..face.len()).for_each(|i| {
63                let mut edge = [face[i], face[(i + 1) % face.len()]];
64                edge.sort_unstable();
65                *edges.entry(edge).or_insert(0u8) += 1;
66            })
67        });
68        if edges.values().any(|&count| count != 2) {
69            return Err("non-manifold boundary");
70        }
71        let (connectivities, mut coordinates) = mesh.into();
72        let mut connectivity = Vec::try_from(connectivities)?;
73        let mut projection = HashMap::new();
74        boundary.iter().flatten().try_for_each(|&node| {
75            if let Entry::Vacant(slot) = projection.entry(node) {
76                let point = bvh
77                    .closest_point(&coordinates[node], surface_coordinates, &elements)
78                    .ok_or("empty tessellation")?
79                    .0;
80                slot.insert(coordinates.len());
81                coordinates.push(point);
82            }
83            Ok(())
84        })?;
85        boundary.iter().for_each(|face| {
86            let [a, b, c, d] = [face[0], face[1], face[2], face[3]];
87            connectivity.push([
88                a,
89                b,
90                c,
91                d,
92                projection[&a],
93                projection[&b],
94                projection[&c],
95                projection[&d],
96            ])
97        });
98        Ok((
99            vec![Connectivity::Hexahedral(connectivity.into())],
100            coordinates,
101        )
102            .into())
103    }
104    fn trim(&self, mesh: &mut Mesh<D>, bvh: &BoundingVolumeHierarchy<D>) {
105        let surface = self.mesh();
106        let surface_coordinates = surface.coordinates();
107        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
108        let normals: Vec<&Coordinate<D>> = self.normals().iter().flatten().collect();
109        let directions = DIRECTIONS.map(|direction| direction.normalized());
110        let coordinates = mesh.coordinates();
111        let number_of_nodes = coordinates.len();
112        let mut inside = vec![false; number_of_nodes];
113        let mut clearance = vec![0.0; number_of_nodes];
114        let threads = available_parallelism().map_or(1, |threads| threads.get());
115        let chunk_size = number_of_nodes.div_ceil(threads).max(1);
116        scope(|scope| {
117            let (elements, normals, directions) = (&elements, &normals, &directions);
118            inside
119                .chunks_mut(chunk_size)
120                .zip(clearance.chunks_mut(chunk_size))
121                .enumerate()
122                .for_each(|(chunk, (flags, distances))| {
123                    scope.spawn(move || {
124                        let offset = chunk * chunk_size;
125                        flags
126                            .iter_mut()
127                            .zip(distances.iter_mut())
128                            .enumerate()
129                            .for_each(|(local, (flag, distance))| {
130                                let point = &coordinates[offset + local];
131                                *flag = directions
132                                    .iter()
133                                    .find_map(|direction| {
134                                        let ray = (point.clone(), direction.clone()).into();
135                                        match bvh.intersect(&ray, surface_coordinates, elements) {
136                                            None => Some(false),
137                                            Some(hit) => {
138                                                let normal = normals[hit.index()];
139                                                let cosine = (direction * normal) / normal.norm();
140                                                (cosine.abs() > GRAZING_TOLERANCE)
141                                                    .then_some(cosine > 0.0)
142                                            }
143                                        }
144                                    })
145                                    .unwrap_or(false);
146                                if *flag
147                                    && let Some((closest, _)) =
148                                        bvh.closest_point(point, surface_coordinates, elements)
149                                {
150                                    *distance = (&closest - point).norm();
151                                }
152                            });
153                    });
154                });
155        });
156        let mut remap = vec![usize::MAX; inside.len()];
157        let mut coordinates = Coordinates::new();
158        let mut connectivity = Vec::new();
159        mesh.iter()
160            .flatten()
161            .filter(|element| {
162                element.iter().all(|&node| inside[node]) && {
163                    let margin = TRIM_MARGIN
164                        * EDGES
165                            .iter()
166                            .map(|&[a, b]| {
167                                (&mesh.coordinates()[element[a]] - &mesh.coordinates()[element[b]])
168                                    .norm()
169                            })
170                            .fold(Scalar::INFINITY, Scalar::min);
171                    element.iter().all(|&node| clearance[node] >= margin)
172                }
173            })
174            .for_each(|element| {
175                connectivity.push(from_fn(|i| {
176                    let node = element[i];
177                    if remap[node] == usize::MAX {
178                        remap[node] = coordinates.len();
179                        coordinates.push(mesh.coordinates()[node].clone());
180                    }
181                    remap[node]
182                }))
183            });
184        *mesh = (
185            vec![Connectivity::Hexahedral(connectivity.into())],
186            coordinates,
187        )
188            .into();
189    }
190}
conspire/geometry/mesh/tessellation/dual/mod.rs

conspire/geometry/mesh/tessellation/dual/
mod.rs
