conspire/geometry/mesh/buffer/
mod.rs

1#[cfg(test)]
2mod test;
3
4mod fit;
5mod restrict;
6
7use super::{Connectivity, Mesh, PrimitiveConnectivity, Tessellation};
8use crate::{
9    geometry::Coordinates,
10    math::{Tensor, TensorVec},
11};
12use std::{
13    array::from_fn,
14    collections::{HashMap, HashSet, hash_map::Entry},
15};
16
17/// The four faces of a tetrahedron, as triples of local node indices.
18const TET_FACES: [[usize; 3]; 4] = [[0, 2, 1], [0, 1, 3], [1, 2, 3], [2, 0, 3]];
19
20/// A mesh peeled open along its boundary, with that boundary's nodes
21/// duplicated so a layer of elements can span the two copies.
22struct Peeled {
23    connectivities: Vec<Connectivity>,
24    coordinates: Coordinates<3>,
25    count: usize,
26    duplicates: HashMap<usize, usize>,
27    layer: Vec<usize>,
28}
29
30/// Splits the prism standing on an outward boundary triangle into three
31/// tetrahedra.
32///
33/// Each lateral quadrilateral is cut by the diagonal running from its
34/// lower-numbered base node to the duplicate of the higher one. That depends
35/// on the edge's two nodes alone, so the prisms either side of a boundary edge
36/// cut their shared quadrilateral the same way, and since node numbers are a
37/// total order the three diagonals can never wind around the prism, which is
38/// what would leave it untetrahedralizable without a new node.
39fn prism(face: &[usize], duplicates: &HashMap<usize, usize>) -> [[usize; 4]; 3] {
40    let first = (0..3)
41        .min_by_key(|&i| face[i])
42        .expect("empty boundary face");
43    let [p0, p1, p2]: [usize; 3] = from_fn(|i| face[(first + i) % 3]);
44    let [q0, q1, q2] = [p0, p1, p2].map(|node| duplicates[&node]);
45    if p1 < p2 {
46        [[p0, p1, p2, q2], [p0, p1, q2, q1], [p0, q1, q2, q0]]
47    } else {
48        [[p0, p1, p2, q1], [p0, p2, q2, q1], [p0, q2, q0, q1]]
49    }
50}
51
52/// Adds `cells` to the last block of their own kind, or as a new one.
53fn merge<const N: usize>(
54    connectivities: &mut Vec<Connectivity>,
55    cells: Vec<[usize; N]>,
56    of_kind: fn(&Connectivity) -> bool,
57    variant: fn(PrimitiveConnectivity<3, N>) -> Connectivity,
58) -> Result<(), &'static str>
59where
60    PrimitiveConnectivity<3, N>: TryFrom<Connectivity, Error = &'static str>,
61{
62    match connectivities.iter().rposition(of_kind) {
63        Some(index) => {
64            let block = PrimitiveConnectivity::<3, N>::try_from(connectivities.remove(index))?;
65            connectivities.insert(
66                index,
67                variant(block.into_iter().chain(cells).collect::<Vec<_>>().into()),
68            )
69        }
70        None => connectivities.push(variant(cells.into())),
71    }
72    Ok(())
73}
74
75/// Drops tetrahedra until the mesh boundary is edge-manifold: every boundary
76/// edge carried by exactly two boundary faces.
77///
78/// [`trim`](Tessellation::trim) keeps or discards a background cell by the
79/// signed distances at its nodes alone, with no topological guard, so a
80/// tetrahedral background can be left pinched along an edge or hanging by one.
81/// Such a boundary cannot be [peeled](Mesh::peel), so each tetrahedron that
82/// carries a boundary face on a non-manifold edge is removed, and the check
83/// repeated, until the boundary closes up. A hexahedral background trimmed the
84/// same way rarely pinches, and [`buffer`](Mesh::buffer) leaves it alone.
85fn manifold_boundary(mut mesh: Mesh<3>) -> Result<Mesh<3>, &'static str> {
86    for _ in 0..64 {
87        let tets: Vec<[usize; 4]> = mesh
88            .iter()
89            .flatten()
90            .map(|tet| from_fn(|i| tet[i]))
91            .collect();
92        let mut face_tets: HashMap<[usize; 3], Vec<usize>> = HashMap::new();
93        for (element, tet) in tets.iter().enumerate() {
94            for face in TET_FACES {
95                let mut key = face.map(|node| tet[node]);
96                key.sort_unstable();
97                face_tets.entry(key).or_default().push(element);
98            }
99        }
100        let mut edge_faces: HashMap<[usize; 2], u32> = HashMap::new();
101        for (face, owners) in &face_tets {
102            if owners.len() == 1 {
103                for [a, b] in [[0, 1], [1, 2], [2, 0]] {
104                    let mut edge = [face[a], face[b]];
105                    edge.sort_unstable();
106                    *edge_faces.entry(edge).or_insert(0) += 1;
107                }
108            }
109        }
110        let bad: HashSet<[usize; 2]> = edge_faces
111            .into_iter()
112            .filter_map(|(edge, count)| (count != 2).then_some(edge))
113            .collect();
114        if bad.is_empty() {
115            return Ok(mesh);
116        }
117        let mut discard: HashSet<usize> = HashSet::new();
118        for (face, owners) in &face_tets {
119            if owners.len() == 1
120                && [[0, 1], [1, 2], [2, 0]].iter().any(|&[a, b]| {
121                    let mut edge = [face[a], face[b]];
122                    edge.sort_unstable();
123                    bad.contains(&edge)
124                })
125            {
126                discard.insert(owners[0]);
127            }
128        }
129        if discard.is_empty() {
130            return Err("non-manifold boundary");
131        }
132        let mut element = 0;
133        mesh.retain_elements(|_, _, _| {
134            let keep = !discard.contains(&element);
135            element += 1;
136            keep
137        });
138    }
139    Err("non-manifold boundary")
140}
141
142/// Constraint on how the buffer layer approaches the target surface.
143#[derive(Clone, Copy, Debug)]
144pub enum Fitting {
145    /// The layer settles wherever the quality and fit energies balance.
146    Soft,
147    /// The layer settles as above, but is then projected onto the surface,
148    /// after which the interior relaxes.
149    Snap,
150}
151
152impl Mesh<3> {
153    pub fn buffer(mut self, target: &Tessellation, fitting: Fitting) -> Result<Self, &'static str> {
154        self.restrict()?;
155        let boundary = self.exterior_faces();
156        let Peeled {
157            mut connectivities,
158            coordinates,
159            count,
160            duplicates,
161            layer,
162        } = self.peel(&boundary, 4, "non-quadrilateral boundary face")?;
163        let cells = boundary
164            .iter()
165            .map(|face| {
166                [
167                    face[0],
168                    face[1],
169                    face[2],
170                    face[3],
171                    duplicates[&face[0]],
172                    duplicates[&face[1]],
173                    duplicates[&face[2]],
174                    duplicates[&face[3]],
175                ]
176            })
177            .collect::<Vec<_>>();
178        merge(
179            &mut connectivities,
180            cells,
181            |connectivity| matches!(connectivity, Connectivity::Hexahedral(_)),
182            Connectivity::Hexahedral,
183        )?;
184        let mut mesh = Self::from((connectivities, coordinates));
185        let nodes: Vec<usize> = layer.iter().copied().chain(0..count).collect();
186        mesh.fit(&nodes, target)?;
187        if let Fitting::Snap = fitting {
188            mesh.project(target, &layer)?;
189            mesh.fit(&(0..count).collect::<Vec<_>>(), target)?;
190        }
191        Ok(mesh)
192    }
193    /// Adds a buffer layer of tetrahedra to a tetrahedral mesh and fits it to
194    /// the target.
195    ///
196    /// The counterpart of [`buffer`](Self::buffer), differing in that each
197    /// boundary triangle raises a prism split into three tetrahedra rather
198    /// than one hexahedron, so the result stays a single tetrahedral block.
199    ///
200    /// It also runs no clearance pre-pass. [`restrict`](Self::restrict) is
201    /// defined on hexahedral boundary quadrilaterals and has no tetrahedral
202    /// analogue yet, so a boundary leaving some node no feasible direction is
203    /// fitted here rather than pruned first.
204    pub fn buffer_tets(
205        self,
206        target: &Tessellation,
207        fitting: Fitting,
208    ) -> Result<Self, &'static str> {
209        let cleaned = manifold_boundary(self)?;
210        let boundary = cleaned.exterior_faces();
211        let Peeled {
212            mut connectivities,
213            coordinates,
214            count,
215            duplicates,
216            layer,
217        } = cleaned.peel(&boundary, 3, "non-triangular boundary face")?;
218        let cells: Vec<[usize; 4]> = boundary
219            .iter()
220            .flat_map(|face| prism(face, &duplicates))
221            .collect();
222        merge(
223            &mut connectivities,
224            cells,
225            |connectivity| matches!(connectivity, Connectivity::Tetrahedral(_)),
226            Connectivity::Tetrahedral,
227        )?;
228        let mut mesh = Self::from((connectivities, coordinates));
229        let nodes: Vec<usize> = layer.iter().copied().chain(0..count).collect();
230        mesh.fit(&nodes, target)?;
231        if let Fitting::Snap = fitting {
232            mesh.project(target, &layer)?;
233            mesh.fit(&(0..count).collect::<Vec<_>>(), target)?;
234        }
235        Ok(mesh)
236    }
237    /// Checks the boundary is a manifold of `arity`-node faces and duplicates
238    /// its nodes.
239    fn peel(
240        self,
241        boundary: &[Vec<usize>],
242        arity: usize,
243        misshapen: &'static str,
244    ) -> Result<Peeled, &'static str> {
245        let mut edges = HashMap::new();
246        boundary.iter().try_for_each(|face| {
247            if face.len() != arity {
248                return Err(misshapen);
249            }
250            (0..arity).for_each(|i| {
251                let mut edge = [face[i], face[(i + 1) % arity]];
252                edge.sort_unstable();
253                *edges.entry(edge).or_insert(0u8) += 1;
254            });
255            Ok(())
256        })?;
257        if edges.values().any(|&count| count != 2) {
258            return Err("non-manifold boundary");
259        }
260        let (connectivities, mut coordinates) = self.into();
261        let connectivities = connectivities.into_members();
262        let count = coordinates.len();
263        let mut duplicates = HashMap::new();
264        let mut layer = Vec::new();
265        boundary.iter().flatten().for_each(|&node| {
266            if let Entry::Vacant(slot) = duplicates.entry(node) {
267                slot.insert(coordinates.len());
268                layer.push(coordinates.len());
269                let point = coordinates[node].clone();
270                coordinates.push(point);
271            }
272        });
273        Ok(Peeled {
274            connectivities,
275            coordinates,
276            count,
277            duplicates,
278            layer,
279        })
280    }
281    /// Moves the layer's nodes onto the closest point of the target.
282    fn project(&mut self, target: &Tessellation, layer: &[usize]) -> Result<(), &'static str> {
283        let surface = target.mesh();
284        let surface_coordinates = surface.coordinates();
285        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
286        let bvh = target.bvh();
287        let coordinates = self.coordinates.members_mut();
288        layer.iter().try_for_each(|&node| {
289            let (point, _) = bvh
290                .closest_point(&coordinates[node], surface_coordinates, &elements)
291                .ok_or("empty tessellation")?;
292            coordinates[node] = point;
293            Ok(())
294        })
295    }
296}
conspire/geometry/mesh/buffer/mod.rs

conspire/geometry/mesh/buffer/
mod.rs
