conspire/geometry/mesh/tessellation/cut/classify/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{Class, DIRECTIONS, topology::element_faces};
5use crate::{
6    geometry::{
7        Coordinate, Coordinates, CoordinatesRef, Direction, DirectionsRef,
8        bbox::BoundingBox,
9        mesh::{Mesh, tessellation::D, tessellation::Tessellation},
10    },
11    math::Tensor,
12};
13use std::collections::{HashMap, hash_map::Entry};
14
15impl Tessellation {
16    pub fn classify(&self, mesh: &Mesh<D>) -> Vec<Class> {
17        let surface = self.mesh();
18        let surface_coordinates = surface.coordinates();
19        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
20        let normals: DirectionsRef<'_, D> = self.normals().iter().flatten().collect();
21        let directions = DIRECTIONS.map(|direction| direction.normalized());
22        let bvh = self.bvh();
23        let coordinates = mesh.coordinates();
24        let number_of_elements = mesh.number_of_elements();
25        let mut cut = vec![false; number_of_elements];
26        mesh.iter()
27            .flat_map(|block| {
28                block
29                    .iter()
30                    .map(move |element| block.element_nodes(element))
31            })
32            .zip(cut.iter_mut())
33            .for_each(|(nodes, flag)| {
34                let bbox: BoundingBox<D> = nodes
35                    .iter()
36                    .map(|&node| &coordinates[node])
37                    .collect::<CoordinatesRef<'_, D>>()
38                    .into();
39                *flag = bvh.overlapping(&bbox).into_iter().any(|triangle| {
40                    let nodes = elements[triangle];
41                    bbox.overlaps_triangle(
42                        &surface_coordinates[nodes[0]],
43                        &surface_coordinates[nodes[1]],
44                        &surface_coordinates[nodes[2]],
45                    )
46                })
47            });
48        let mut faces = HashMap::new();
49        let mut neighbors: Vec<Vec<usize>> = vec![Vec::new(); number_of_elements];
50        let mut offset = 0;
51        mesh.iter().for_each(|block| {
52            block.iter().enumerate().for_each(|(local, element)| {
53                let index = offset + local;
54                if !cut[index] {
55                    element_faces(block, element).into_iter().for_each(|face| {
56                        let mut key = face;
57                        key.sort_unstable();
58                        match faces.entry(key) {
59                            Entry::Occupied(entry) => {
60                                let other = *entry.get();
61                                neighbors[index].push(other);
62                                neighbors[other].push(index);
63                            }
64                            Entry::Vacant(slot) => {
65                                slot.insert(index);
66                            }
67                        }
68                    })
69                }
70            });
71            offset += block.number_of_elements();
72        });
73        let centroids = mesh.centroids();
74        let mut classes: Vec<Class> = cut
75            .iter()
76            .map(|&flag| if flag { Class::Cut } else { Class::Outside })
77            .collect();
78        let mut visited = cut;
79        let mut stack = Vec::new();
80        (0..number_of_elements).for_each(|seed| {
81            if !visited[seed] {
82                let class = if self.encloses(
83                    &centroids[seed],
84                    surface_coordinates,
85                    &elements,
86                    &normals,
87                    &directions,
88                ) {
89                    Class::Inside
90                } else {
91                    Class::Outside
92                };
93                visited[seed] = true;
94                stack.push(seed);
95                while let Some(index) = stack.pop() {
96                    classes[index] = class;
97                    neighbors[index].iter().for_each(|&next| {
98                        if !visited[next] {
99                            visited[next] = true;
100                            stack.push(next);
101                        }
102                    })
103                }
104            }
105        });
106        classes
107    }
108    pub(super) fn encloses(
109        &self,
110        point: &Coordinate<D>,
111        surface_coordinates: &Coordinates<D>,
112        elements: &[&[usize]],
113        normals: &DirectionsRef<'_, D>,
114        directions: &[Direction<D>; 3],
115    ) -> bool {
116        directions
117            .iter()
118            .find_map(|direction| {
119                let ray = (point.clone(), direction.clone()).into();
120                match self.bvh().intersect(&ray, surface_coordinates, elements) {
121                    None => Some(false),
122                    Some(hit) => {
123                        let normal = &normals[hit.index()];
124                        let cosine = (direction * normal) / normal.norm();
125                        (cosine.abs() > super::GRAZING_TOLERANCE).then_some(cosine > 0.0)
126                    }
127                }
128            })
129            .unwrap_or(false)
130    }
131}
conspire/geometry/mesh/tessellation/cut/classify/mod.rs

conspire/geometry/mesh/tessellation/cut/classify/
mod.rs
