conspire/geometry/mesh/quality/improve/untangle/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::Incidence;
5use crate::{
6    geometry::{
7        Coordinate,
8        mesh::{Mesh, Tessellation},
9    },
10    math::{Scalar, Tensor},
11};
12use std::mem::transmute_copy;
13
14const PROBES: usize = 32;
15
16impl<const D: usize> Mesh<D> {
17    pub fn untangle(&mut self, iterations: usize, margin: Scalar, surface: Option<&Tessellation>) {
18        let number_of_nodes = self.number_of_nodes();
19        let neighbors = self.node_node_connectivity().to_vec();
20        let incidence = Incidence::of(self);
21        let constrained = surface.filter(|_| D == 3);
22        let elements: Vec<&[usize]> = constrained
23            .map(|surface| surface.mesh().connectivities().iter().flatten().collect())
24            .unwrap_or_default();
25        let mut boundary = vec![false; number_of_nodes];
26        if constrained.is_some() {
27            self.exterior_faces()
28                .iter()
29                .flatten()
30                .for_each(|&node| boundary[node] = true);
31        }
32        let coordinates = self.coordinates.members_mut();
33        for _ in 0..iterations {
34            for node in 0..number_of_nodes {
35                if neighbors[node].is_empty() {
36                    continue;
37                }
38                let mut current = incidence.inversion(node, coordinates, margin);
39                if current <= 0.0 {
40                    continue;
41                }
42                let mut step = 0.5
43                    * neighbors[node]
44                        .iter()
45                        .map(|&neighbor| (&coordinates[node] - &coordinates[neighbor]).norm())
46                        .sum::<Scalar>()
47                    / (neighbors[node].len() as Scalar);
48                for _ in 0..PROBES {
49                    let mut improved = false;
50                    for axis in 0..D {
51                        for sign in [-1.0, 1.0] {
52                            let original = coordinates[node].clone();
53                            coordinates[node][axis] += sign * step;
54                            if boundary[node] {
55                                coordinates[node] = project_to_surface(
56                                    constrained.unwrap(),
57                                    &elements,
58                                    &coordinates[node],
59                                );
60                            }
61                            let trial = incidence.inversion(node, coordinates, margin);
62                            if trial < current {
63                                current = trial;
64                                improved = true;
65                            } else {
66                                coordinates[node] = original;
67                            }
68                        }
69                    }
70                    if !improved {
71                        step *= 0.5;
72                    }
73                }
74            }
75        }
76    }
77}
78
79fn project_to_surface<const D: usize>(
80    surface: &Tessellation,
81    elements: &[&[usize]],
82    point: &Coordinate<D>,
83) -> Coordinate<D> {
84    let query: &Coordinate<3> = unsafe { &*(point as *const Coordinate<D>).cast() };
85    let projected = surface
86        .bvh()
87        .closest_point(query, surface.mesh().coordinates(), elements)
88        .map_or_else(|| query.clone(), |(projected, _)| projected);
89    unsafe { transmute_copy(&projected) }
90}
conspire/geometry/mesh/quality/improve/untangle/mod.rs

conspire/geometry/mesh/quality/improve/untangle/
mod.rs
