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}