Skip to main content

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

1#[cfg(test)]
2mod test;
3
4mod assemble;
5mod build;
6mod classify;
7mod cleanup;
8mod face;
9mod geometry;
10mod lattice;
11mod snap;
12mod split;
13mod tables;
14mod topology;
15
16use crate::{
17    geometry::{
18        Coordinate, Direction,
19        mesh::{
20            Dualization, Mesh,
21            tessellation::{D, Tessellation, cut::geometry::contained},
22        },
23        ntree::{Balance, Balancing, CurvatureSizing, Octree, Pairing, Sizing},
24    },
25    math::{Quantity, Scalar},
26    units::Length,
27};
28use std::{collections::HashMap, num::NonZeroU32};
29
30const COLLAPSE_FRACTION: Scalar = 0.2;
31const CROSSING_TOLERANCE: Quantity<Length> = Length::meters(1.0e-8);
32const GRAZING_TOLERANCE: Scalar = 1.0e-4;
33const PADDING: u16 = 2;
34const SLIVER_FRACTION: Scalar = 0.1;
35const SNAP_FEATURE: Scalar = 0.5;
36const SNAP_HARD: Scalar = 0.05;
37const SNAP_QUALITY: Scalar = 0.3;
38const SNAP_SOFT: Scalar = 0.2;
39const FACES: [[usize; 4]; 6] = [
40    [0, 1, 5, 4],
41    [1, 2, 6, 5],
42    [2, 3, 7, 6],
43    [3, 0, 4, 7],
44    [0, 3, 2, 1],
45    [4, 5, 6, 7],
46];
47const EDGES: [[usize; 2]; 12] = [
48    [0, 1],
49    [1, 2],
50    [2, 3],
51    [3, 0],
52    [4, 5],
53    [5, 6],
54    [6, 7],
55    [7, 4],
56    [0, 4],
57    [1, 5],
58    [2, 6],
59    [3, 7],
60];
61const DIRECTIONS: [Direction<D>; 3] = [
62    Direction::const_from([1.0, 0.140_412_03, 0.092_153_88]),
63    Direction::const_from([0.097_153_2, 1.0, 0.131_771_4]),
64    Direction::const_from([0.123_456_7, 0.087_654_3, 1.0]),
65];
66
67/// What an octree background's cells are meshed into.
68enum Cells {
69    Polyhedral,
70    Tetrahedral,
71}
72
73#[derive(Clone, Copy, Debug, PartialEq)]
74pub enum Class {
75    Inside,
76    Cut,
77    Outside,
78}
79
80#[derive(Clone, Copy, Debug, PartialEq)]
81pub enum Sign {
82    Inside,
83    On,
84    Outside,
85}
86
87/// Identifies a point used while stitching a cut face/cell.
88///
89/// Either an original mesh node, or the `ordinal`-th crossing,
90/// (in canonical ascending-node-order direction) along the sorted edge.
91#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
92pub enum Vertex {
93    Node(usize),
94    Crossing([usize; 2], usize),
95}
96
97pub struct Tables {
98    signs: HashMap<usize, Sign>,
99    crossings: HashMap<[usize; 2], Vec<Coordinate<D>>>,
100    faces: HashMap<[usize; 4], [usize; 4]>,
101    segments: HashMap<[usize; 4], Vec<[Vertex; 2]>>,
102}
103
104impl Tables {
105    pub fn signs(&self) -> &HashMap<usize, Sign> {
106        &self.signs
107    }
108    pub fn crossings(&self) -> &HashMap<[usize; 2], Vec<Coordinate<D>>> {
109        &self.crossings
110    }
111    pub fn faces(&self) -> &HashMap<[usize; 4], [usize; 4]> {
112        &self.faces
113    }
114    pub fn segments(&self) -> &HashMap<[usize; 4], Vec<[Vertex; 2]>> {
115        &self.segments
116    }
117}
118
119impl Tessellation {
120    /// Builds the dual of an octree fitted to this tessellation, with each
121    /// cell classified against the surface.
122    ///
123    /// The background for [`cut`](Self::cut). `balancing` must be `Strong(1)`
124    /// or `Weak(1)`, which is what dualization requires.
125    pub fn dual_background(
126        &self,
127        balancing: Balancing,
128        scale: Scalar,
129    ) -> Result<(Mesh<D>, Vec<Class>), &'static str> {
130        let sizing = Sizing::new(self, scale, CurvatureSizing::default(), PADDING);
131        let mesh = if sizing.fits::<u16>() {
132            let mut octree = Octree::<u16, NonZeroU32>::refine(&sizing)?;
133            octree.equilibrate(balancing, Pairing::Regular)?;
134            octree.dualize()
135        } else {
136            let mut octree = Octree::<u32, NonZeroU32>::refine(&sizing)?;
137            octree.equilibrate(balancing, Pairing::Regular)?;
138            octree.dualize()
139        };
140        let classes = self.classify(&mesh);
141        Ok((mesh, classes))
142    }
143    /// Builds a uniform lattice of cubes of the given edge length around this
144    /// tessellation, with each cell classified against the surface.
145    ///
146    /// The lattice spans the cells the surface passes through, those its
147    /// interior encloses, and a single shell of cells beyond them, so it is a
148    /// background to be [cut](Self::cut), or [trimmed](Self::trim) and
149    /// [buffered](Mesh::buffer), rather than a finished mesh.
150    ///
151    /// Unlike [`dual_background`](Self::dual_background) the cells are all
152    /// axis-aligned cubes, at the cost of the grading a tree provides, and
153    /// the classes fall out of rasterizing rather than being found again.
154    pub fn lattice_background(
155        &self,
156        spacing: Quantity<Length>,
157    ) -> Result<(Mesh<D>, Vec<Class>), &'static str> {
158        Ok(self.lattice_cells(spacing)?.mesh())
159    }
160    /// Builds a uniform lattice around this tessellation and splits every cell
161    /// into six tetrahedra, with each one classified against the surface.
162    ///
163    /// The tetrahedral counterpart of
164    /// [`lattice_background`](Self::lattice_background). The cells are still
165    /// classified by rasterizing, so the six tetrahedra of a cell all take the
166    /// class of the cell they came from.
167    pub fn lattice_tet_background(
168        &self,
169        spacing: Quantity<Length>,
170    ) -> Result<(Mesh<D>, Vec<Class>), &'static str> {
171        Ok(self.lattice_cells(spacing)?.tets())
172    }
173    /// Builds an octree fitted to this tessellation, with each cell
174    /// classified against the surface.
175    ///
176    /// The background for [`cut_polyhedral`](Self::cut_polyhedral), taking
177    /// the octree directly rather than its dual. This places no 2:1
178    /// requirement on `balancing`, since hanging nodes become extra vertices
179    /// on a face rather than something to be dualized away. `Weak(n)` and
180    /// `Strong(n)` for `n > 1` are therefore available here, permitting
181    /// coarser trees than dualization allows.
182    pub fn octree_background(
183        &self,
184        balancing: Balancing,
185        scale: Scalar,
186    ) -> Result<(Mesh<D>, Vec<Class>), &'static str> {
187        let mesh = self.octree_mesh(
188            balancing,
189            Pairing::Regular,
190            scale,
191            CurvatureSizing::default(),
192            Cells::Polyhedral,
193        )?;
194        let classes = self.classify(&mesh);
195        Ok((mesh, classes))
196    }
197    /// Builds an octree fitted to this tessellation and meshes it as
198    /// tetrahedra, with each one classified against the surface.
199    ///
200    /// The tetrahedral counterpart of
201    /// [`octree_background`](Self::octree_background), to be
202    /// [trimmed](Self::trim). `balancing` must be `Strong(1)`: the templates
203    /// filling a graded cell only span a one-level difference, and only a
204    /// balance over edges and vertices as well as faces holds them to it.
205    /// `pairing` need not be `Regular`; the tetrahedra conform under any
206    /// pairing, and `None` yields a smaller background.
207    ///
208    /// `tolerance` is the Dunyach chord-error tolerance for curvature-driven
209    /// refinement; `None` disables it.
210    pub fn octree_tet_background(
211        &self,
212        balancing: Balancing,
213        pairing: Pairing,
214        scale: Scalar,
215        tolerance: Option<Quantity<Length>>,
216    ) -> Result<(Mesh<D>, Vec<Class>), &'static str> {
217        if !matches!(balancing, Balancing::Strong(1)) {
218            return Err("tetrahedra require Strong(1) balancing");
219        }
220        let curvature = CurvatureSizing {
221            tolerance,
222            ..Default::default()
223        };
224        let mesh = self.octree_mesh(balancing, pairing, scale, curvature, Cells::Tetrahedral)?;
225        let classes = self.classify(&mesh);
226        Ok((mesh, classes))
227    }
228    fn octree_mesh(
229        &self,
230        balancing: Balancing,
231        pairing: Pairing,
232        scale: Scalar,
233        curvature: CurvatureSizing,
234        cells: Cells,
235    ) -> Result<Mesh<D>, &'static str> {
236        let sizing = Sizing::new(self, scale, curvature, PADDING);
237        if sizing.fits::<u16>() {
238            let mut octree = Octree::<u16, NonZeroU32>::refine(&sizing)?;
239            octree.equilibrate(balancing, pairing)?;
240            Ok(match cells {
241                Cells::Polyhedral => Mesh::from(octree),
242                Cells::Tetrahedral => Mesh::tetrahedra_from(octree),
243            })
244        } else {
245            let mut octree = Octree::<u32, NonZeroU32>::refine(&sizing)?;
246            octree.equilibrate(balancing, pairing)?;
247            Ok(match cells {
248                Cells::Polyhedral => Mesh::from(octree),
249                Cells::Tetrahedral => Mesh::tetrahedra_from(octree),
250            })
251        }
252    }
253    /// Cuts a classified background mesh to this tessellation, leaving
254    /// hexahedra everywhere but at the boundary.
255    ///
256    /// Snaps the nodes that nearly lie on the surface onto it, builds the
257    /// crossing tables, and assembles the cut cells into polyhedra.
258    pub fn cut(&self, mesh: Mesh<D>, classes: &[Class]) -> Result<Mesh<D>, &'static str> {
259        if !contained(&mesh, classes) {
260            return Err("tessellation is not contained within the background mesh");
261        }
262        let (mesh, snapped) = self.snap(mesh, classes)?;
263        let tables = self.tables(&mesh, classes, &snapped)?;
264        self.assemble(&mesh, classes, &tables)
265    }
266    /// Cuts a classified background mesh to this tessellation, leaving
267    /// polyhedra throughout.
268    ///
269    /// The counterpart of [`cut`](Self::cut) for a background whose cells
270    /// carry hanging nodes, such as an octree taken directly.
271    pub fn cut_polyhedral(
272        &self,
273        mesh: Mesh<D>,
274        classes: &[Class],
275    ) -> Result<Mesh<D>, &'static str> {
276        if !contained(&mesh, classes) {
277            return Err("tessellation is not contained within the background mesh");
278        }
279        let (mesh, snapped) = self.snap_generic(mesh, classes)?;
280        let tables = self.tables_generic(&mesh, classes, &snapped)?;
281        self.assemble_generic(&mesh, classes, &tables)
282    }
283}