Skip to main content

parry3d/transformation/volume_mesh/
mod.rs

1//! Filling a closed boundary with a simulation-grade simplex mesh.
2
3use crate::math::{Real, Vector, DIM};
4use alloc::vec::Vec;
5
6#[cfg(feature = "dim3")]
7mod cover_octree;
8#[cfg(feature = "dim3")]
9mod cover_smoothing;
10#[cfg(feature = "dim2")]
11mod delaunay_refinement;
12#[cfg(feature = "dim3")]
13mod isosurface_stuffing;
14
15/// A simplicial mesh filling the interior of a closed boundary: triangles in 2D, tetrahedra in 3D.
16#[derive(Clone, Debug, Default)]
17#[cfg_attr(
18    feature = "serde-serialize",
19    derive(serde::Serialize, serde::Deserialize)
20)]
21pub struct VolumeMesh {
22    /// The mesh vertices.
23    pub vertices: Vec<Vector>,
24    /// The mesh cells, positively oriented (positive signed area/volume).
25    pub cells: Vec<[u32; DIM + 1]>,
26}
27
28/// Which cover the mesh is (3D only): of the shape's volume, or of its surface alone.
29///
30/// Either way the mesh is made of whole lattice cells the shape reaches, so it *contains*
31/// what it covers, which is what an embedding wants: every point of the shape (a skin
32/// vertex, say) interpolates inside a cell instead of extrapolating outside all of them.
33#[cfg(feature = "dim3")]
34#[derive(Copy, Clone, Debug, PartialEq, Eq, Default)]
35pub enum MeshEnclosure {
36    /// Every lattice cell the shape reaches is kept whole: a solid fill containing the shape by
37    /// construction, blocky at the cell size (see [`VolumeMeshParameters::cover_smoothing`] and
38    /// [`VolumeMeshParameters::cover_subdivisions`]). Needs a closed, consistently oriented mesh.
39    #[default]
40    Cover,
41    /// The cover of the shape's surface alone: only the cells the surface crosses are kept,
42    /// the interior stays empty, and the mesh need not be closed or oriented. The result is a
43    /// hollow shell that deforms like a shell, not like a solid.
44    Crust,
45}
46
47/// Parameters of [`volume_mesh`].
48#[derive(Copy, Clone, Debug, PartialEq)]
49pub struct VolumeMeshParameters {
50    /// Target size of the generated elements.
51    pub cell_size: Real,
52    /// Smallest angle, in radians, the refinement aims for (2D only).
53    ///
54    /// Above 30 degrees the refinement is not guaranteed to terminate, so it may stop early and
55    /// leave some elements badly shaped.
56    #[cfg(feature = "dim2")]
57    pub min_angle: Real,
58    /// Whether the whole shape is covered or its surface alone (3D only); see
59    /// [`MeshEnclosure`].
60    #[cfg(feature = "dim3")]
61    pub enclosure: MeshEnclosure,
62    /// How many shrink-wrap iterations smooth the staircase of a [`MeshEnclosure::Cover`] mesh,
63    /// zero leaving it raw (3D only, cover only). Each iteration pulls the boundary toward
64    /// the shape, held off by [`Self::cover_guard`], and commits a move only if containment holds.
65    #[cfg(feature = "dim3")]
66    pub cover_smoothing: u32,
67    /// How close to the shape the smoothed cover may pull its boundary, as a fraction of
68    /// the local cell size (3D only, read by the cover's smoothing only).
69    #[cfg(feature = "dim3")]
70    pub cover_guard: Real,
71    /// How many halvings below [`Self::cell_size`] a [`MeshEnclosure::Cover`] cell crossing the
72    /// shape's boundary may be refined, zero keeping the boundary at the cell size (3D only, cover
73    /// only). The refinement is on the octree (no hanging vertex) and runs before the smoothing.
74    #[cfg(feature = "dim3")]
75    pub cover_subdivisions: u32,
76}
77
78impl VolumeMeshParameters {
79    /// Parameters generating elements of size `cell_size`: in 3D a raw cover (no
80    /// smoothing, no subdivision), in 2D a refinement aiming for 30 degree angles.
81    pub fn new(cell_size: Real) -> Self {
82        Self {
83            cell_size,
84            #[cfg(feature = "dim2")]
85            #[cfg_attr(feature = "f64", expect(clippy::unnecessary_cast))]
86            min_angle: core::f64::consts::PI as Real / 6.0,
87            #[cfg(feature = "dim3")]
88            enclosure: MeshEnclosure::Cover,
89            #[cfg(feature = "dim3")]
90            cover_smoothing: 0,
91            #[cfg(feature = "dim3")]
92            cover_guard: 0.15,
93            #[cfg(feature = "dim3")]
94            cover_subdivisions: 0,
95        }
96    }
97}
98
99/// Fills a boundary `(vertices, indices)` (2D polyline, 3D triangle mesh) with a simplex mesh
100/// for finite-element simulation: 2D refines a constrained Delaunay triangulation of the boundary,
101/// 3D builds a lattice cover containing it. Returns `None` if the boundary encloses nothing.
102pub fn volume_mesh(
103    vertices: &[Vector],
104    indices: &[[u32; DIM]],
105    params: &VolumeMeshParameters,
106) -> Option<VolumeMesh> {
107    #[cfg(feature = "dim2")]
108    {
109        delaunay_refinement::triangulate(vertices, indices, params)
110    }
111    #[cfg(feature = "dim3")]
112    {
113        isosurface_stuffing::tetrahedralize(vertices, indices, params)
114    }
115}
116
117impl VolumeMesh {
118    /// The connected component of each cell: two cells sharing a vertex are in the same one.
119    ///
120    /// Components are numbered from zero, so the number of them is the largest index plus one.
121    /// A mesh in more than one piece means the input's own surface is in several pieces: the
122    /// cover keeps every cell the surface reaches, however thin the feature.
123    pub fn connected_components(&self) -> Vec<u32> {
124        let mut parent: Vec<u32> = (0..self.vertices.len() as u32).collect();
125
126        fn root(parent: &mut [u32], mut i: u32) -> u32 {
127            while parent[i as usize] != i {
128                parent[i as usize] = parent[parent[i as usize] as usize];
129                i = parent[i as usize];
130            }
131            i
132        }
133
134        for cell in &self.cells {
135            for vid in &cell[1..] {
136                let (a, b) = (root(&mut parent, cell[0]), root(&mut parent, *vid));
137                parent[a as usize] = b;
138            }
139        }
140
141        let mut ids = alloc::vec![u32::MAX; self.vertices.len()];
142        let mut count = 0;
143
144        self.cells
145            .iter()
146            .map(|cell| {
147                let id = &mut ids[root(&mut parent, cell[0]) as usize];
148                if *id == u32::MAX {
149                    *id = count;
150                    count += 1;
151                }
152                *id
153            })
154            .collect()
155    }
156
157    /// Removes the vertices no cell references, and reindexes the cells accordingly.
158    pub fn compact(&mut self) {
159        let mut remap = alloc::vec![u32::MAX; self.vertices.len()];
160        let mut vertices = Vec::with_capacity(self.vertices.len());
161
162        for cell in &mut self.cells {
163            for vid in cell {
164                let remapped = &mut remap[*vid as usize];
165                if *remapped == u32::MAX {
166                    *remapped = vertices.len() as u32;
167                    vertices.push(self.vertices[*vid as usize]);
168                }
169                *vid = *remapped;
170            }
171        }
172
173        self.vertices = vertices;
174    }
175}