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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    geometry::{
6        Direction, DirectionsRef,
7        mesh::{
8            Mesh,
9            tessellation::{D, Tessellation},
10        },
11    },
12    math::{Scalar, Tensor},
13};
14use std::thread::{available_parallelism, scope};
15
16const GRAZING_TOLERANCE: Scalar = 1.0e-4;
17const TRIM_RATIO: Scalar = 0.1;
18const DIRECTIONS: [Direction<D>; 3] = [
19    Direction::const_from([1.0, 0.140_412_03, 0.092_153_88]),
20    Direction::const_from([0.097_153_2, 1.0, 0.131_771_4]),
21    Direction::const_from([0.123_456_7, 0.087_654_3, 1.0]),
22];
23
24impl Tessellation {
25    /// Discards the cells of a background mesh lying outside this
26    /// tessellation, leaving a mesh that covers the volume it encloses.
27    ///
28    /// A cell survives when the signed distances at its nodes satisfy
29    /// `minimum + 0.1 * maximum >= 0`, so the cells straddling the surface
30    /// are kept for [`buffer`](Mesh::buffer) to fit onto it.
31    pub fn trim(&self, mesh: &mut Mesh<D>) -> Result<(), &'static str> {
32        let bvh = self.bvh();
33        let surface = self.mesh();
34        let surface_coordinates = surface.coordinates();
35        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
36        let normals: DirectionsRef<'_, D> = self.normals().iter().flatten().collect();
37        let directions = DIRECTIONS.map(|direction| direction.normalized());
38        let coordinates = mesh.coordinates();
39        let number_of_nodes = coordinates.len();
40        let mut signed = vec![Scalar::NEG_INFINITY; number_of_nodes];
41        let threads = available_parallelism().map_or(1, |threads| threads.get());
42        let chunk_size = number_of_nodes.div_ceil(threads).max(1);
43        scope(|scope| {
44            let (elements, normals, directions) = (&elements, &normals, &directions);
45            signed
46                .chunks_mut(chunk_size)
47                .enumerate()
48                .for_each(|(chunk, distances)| {
49                    scope.spawn(move || {
50                        let offset = chunk * chunk_size;
51                        distances
52                            .iter_mut()
53                            .enumerate()
54                            .for_each(|(local, distance)| {
55                                let point = &coordinates[offset + local];
56                                let inside = directions
57                                    .iter()
58                                    .find_map(|direction| {
59                                        let ray = (point.clone(), direction.clone()).into();
60                                        match bvh.intersect(&ray, surface_coordinates, elements) {
61                                            None => Some(false),
62                                            Some(hit) => {
63                                                let normal = &normals[hit.index()];
64                                                let cosine = (direction * normal) / normal.norm();
65                                                (cosine.abs() > GRAZING_TOLERANCE)
66                                                    .then_some(cosine > 0.0)
67                                            }
68                                        }
69                                    })
70                                    .unwrap_or(false);
71                                if let Some((closest, _)) =
72                                    bvh.closest_point(point, surface_coordinates, elements)
73                                {
74                                    let magnitude = (&closest - point).norm().value();
75                                    *distance = if inside { magnitude } else { -magnitude };
76                                }
77                            });
78                    });
79                });
80        });
81        mesh.retain_elements(|_, element, _| {
82            let (minimum, maximum) = element.iter().fold(
83                (Scalar::INFINITY, Scalar::NEG_INFINITY),
84                |(minimum, maximum), &node| (minimum.min(signed[node]), maximum.max(signed[node])),
85            );
86            minimum + TRIM_RATIO * maximum >= 0.0
87        });
88        Ok(())
89    }
90}
conspire/geometry/mesh/tessellation/trim/mod.rs

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