honeycomb_kernels/grisubal/mod.rs
1//! *grisubal* algorithm description & implementation
2//!
3//! This module contain all code used to implement our grid submersion algorithm, or *grisubal*
4//! for short.
5//!
6//! This algorithm builds the mesh of a geometry by overlapping a grid over it and intersecting
7//! the grid with the geometry. It is inspired by the approach described in
8//! [this](https://internationalmeshingroundtable.com/assets/research-notes/imr32/2011.pdf)
9//! research note.
10//!
11//! # Assumptions / Hypotheses
12//!
13//! Boundaries are consistently oriented, i.e.:
14//! - normals of segments making up a boundary all point outward / inward, no mix,
15//! - boundaries are closed,
16//! - if there are nested boundaries, their orientation are consistent one with the other; this is
17//! an extension of the first condition.
18//!
19//! # Algorithm
20//!
21//! The steps followed by the algorithm are detailed in the user guide. The following is a summary.
22//!
23//! ### Pre-processing
24//!
25//! 1. Compute characteristics of a grid covering the entire geometry, avoiding exact intersection
26//! between the grid's segments and the geometry's vertices.
27//! 2. Remove "redundant" Points of Interest to avoid duplicated vertices.
28//! 3. Check for obvious orientation issues (open geometry & orientation per boundary).
29//!
30//! ### Main kernel
31//!
32//! 1. Compute intersection vertices between the geometry's segments and the grid.
33//! 2. Insert given intersections into the grid.
34//! 3. Build new edge data by searching through the original segments.
35//! 4. Insert new edges into the map. Mark darts on each side of the edge with the `Boundary`
36//! attribute.
37//!
38//! ### Post-processing clip
39//!
40//! Depending on the specified argument, one side (or the other) of the boundary can be clipped.
41//! This is specified using the [`Clip`] enum; The following steps describe the operation for
42//! [`Clip::Left`]:
43//!
44//! 1. Fetch all darts marked as `Boundary::Left` during the last step of the main kernel.
45//! 2. Use these darts' faces as starting point for a coloring algorithm. The search is done using
46//! a BFS and only consider adjacent faces if the adjacent dart isn't marked as a boundary.
47//! This step is also used to check for orientation inconsistencies, most importantly orientation
48//! across distinct boundaries.
49//! 3. Delete all darts making up the marked faces.
50//!
51//! The `Boundary` attribute is then removed from the map before return.
52
53pub(crate) mod model;
54pub(crate) mod routines;
55pub(crate) mod timers;
56
57use honeycomb_core::{
58 cmap::{CMap2, CMapBuilder, GridDescriptor},
59 geometry::CoordsFloat,
60};
61use thiserror::Error;
62use vtkio::Vtk;
63
64use crate::grisubal::{
65 model::{Boundary, Geometry2},
66 routines::{
67 clip_left, clip_right, compute_intersection_ids, compute_overlapping_grid,
68 detect_orientation_issue, generate_edge_data, generate_intersection_data,
69 group_intersections_per_edge, insert_edges_in_map, insert_intersections,
70 remove_redundant_poi,
71 },
72 timers::{finish, start_timer, unsafe_time_section},
73};
74
75/// Post-processing clip operation.
76///
77/// Note that the part of the map that is clipped depends on the orientation of the original geometry provided as
78/// input.
79#[derive(Default)]
80pub enum Clip {
81 /// Clip elements located on the left side of the oriented boundary.
82 Left,
83 /// Clip elements located on the right side of the oriented boundary.
84 Right,
85 /// Keep all elements. Default value.
86 #[default]
87 None,
88}
89
90#[derive(Error, Debug)]
91/// Enum used to model potential errors of the `grisubal` kernel.
92///
93/// Each variant has an associated message that details more precisely what was detected.
94pub enum GrisubalError {
95 /// An orientation issue has been detected in the input geometry.
96 #[error("boundary isn't consistently oriented - {0}")]
97 InconsistentOrientation(&'static str),
98 /// The specified geometry does not match one (or more) requirements of the algorithm.
99 #[error("input shape isn't conform to requirements - {0}")]
100 InvalidShape(&'static str),
101 /// The VTK file used to try to build a `Geometry2` object contains invalid data
102 /// (per VTK's specification).
103 #[error("invalid/corrupted data in the vtk file - {0}")]
104 BadVtkData(&'static str),
105 /// The VTK file used to try to build a `Geometry2` object contains valid but unsupported data.
106 #[error("unsupported data in the vtk file - {0}")]
107 UnsupportedVtkData(&'static str),
108}
109
110#[allow(clippy::missing_errors_doc)]
111/// Main algorithm call function.
112///
113/// # Arguments
114///
115/// - `file_path: impl AsRef<Path>` -- Path to a VTK file describing input geometry.
116/// - `grid_cell_sizes: [T; 2],` -- Desired grid cell size along the X/Y axes.
117/// - `clip: Option<Clip>` -- Indicates which part of the map should be clipped, if any, in
118/// the post-processing phase.
119///
120///
121/// At the moment, the input geometry should be specified via a file under the VTK Legacy format.
122/// Just like the `io` feature provided in the core crate, there are a few additional requirements
123/// for the geometry to be loaded correctly:
124/// - The geometry should have a consistent orientation, i.e. the order in which the points are
125/// given should form normals with a consistent direction (either pointing inward or outward the
126/// geometry).
127/// - The geometry should be described using in an `UnstructuredGrid` data set, with supported
128/// cell types (`Vertex`, `Line`). Lines will be interpreted as the boundary to intersect while
129/// vertices will be considered as points of interests.
130///
131/// # Return / Errors
132///
133/// This function returns a `Result` taking the following values:
134/// - `Ok(CMap2)` -- Algorithm ran successfully.
135/// - `Err(GrisubalError)` -- Algorithm encountered an issue. See [`GrisubalError`] for all
136/// possible errors.
137///
138/// # Panics
139///
140/// This function may panic if the specified file cannot be opened.
141///
142/// # Example
143///
144/// ```no_run
145/// # use honeycomb_core::cmap::CMap2;
146/// # use honeycomb_kernels::grisubal::*;
147/// # fn main() -> Result<(), GrisubalError>{
148/// let cmap: CMap2<f64> = grisubal("some/path/to/geometry.vtk", [1., 1.], Clip::default())?;
149/// # Ok(())
150/// # }
151/// ```
152#[allow(clippy::needless_pass_by_value)]
153pub fn grisubal<T: CoordsFloat>(
154 file_path: impl AsRef<std::path::Path>,
155 grid_cell_sizes: [T; 2],
156 clip: Clip,
157) -> Result<CMap2<T>, GrisubalError> {
158 // INIT TIMER
159 start_timer!(instant);
160
161 // --- IMPORT VTK INPUT
162 let geometry_vtk = match Vtk::import(file_path) {
163 Ok(vtk) => vtk,
164 Err(e) => panic!("E: could not open specified vtk file - {e}"),
165 };
166 unsafe_time_section!(instant, timers::Section::ImportVTK);
167 //----/
168
169 // --- BUILD OUR MODEL FROM THE VTK IMPORT
170 let mut geometry = Geometry2::try_from(geometry_vtk)?;
171 unsafe_time_section!(instant, timers::Section::BuildGeometry);
172 //----/
173
174 // --- FIRST DETECTION OF ORIENTATION ISSUES
175 detect_orientation_issue(&geometry)?;
176 unsafe_time_section!(instant, timers::Section::DetectOrientation);
177 //----/
178
179 // --- FIND AN OVERLAPPING GRID
180 let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?;
181 let [cx, cy] = grid_cell_sizes;
182 let ogrid = GridDescriptor::default()
183 .n_cells_x(nx)
184 .n_cells_y(ny)
185 .len_per_cell_x(cx)
186 .len_per_cell_y(cy)
187 .origin(origin);
188 unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid);
189 //----/
190
191 // --- REMOVE REDUNDANT PoIs
192 remove_redundant_poi(&mut geometry, grid_cell_sizes, origin);
193 unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi);
194 //----/
195
196 // ------ START MAIN KERNEL TIMER
197 start_timer!(kernel);
198
199 // --- BUILD THE GRID
200 let mut cmap = CMapBuilder::default()
201 .grid_descriptor(ogrid)
202 .add_attribute::<Boundary>() // will be used for clipping
203 .build()
204 .expect("E: unreachable"); // unreachable because grid dims are valid
205 unsafe_time_section!(instant, timers::Section::BuildMeshInit);
206 //----/
207
208 // process the geometry
209
210 // --- STEP 1 & 2
211 // (1)
212 let (new_segments, intersection_metadata) =
213 generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin);
214 // (2)
215 let n_intersec = intersection_metadata.len();
216 let (edge_intersec, dart_slices) =
217 group_intersections_per_edge(&mut cmap, intersection_metadata);
218 let intersection_darts = compute_intersection_ids(n_intersec, &edge_intersec, &dart_slices);
219 unsafe_time_section!(instant, timers::Section::BuildMeshIntersecData);
220 //----/
221
222 // --- STEP 3
223 insert_intersections(&cmap, &edge_intersec, &dart_slices);
224 unsafe_time_section!(instant, timers::Section::BuildMeshInsertIntersec);
225 //----/
226
227 // --- STEP 4
228 let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts);
229 unsafe_time_section!(instant, timers::Section::BuildMeshEdgeData);
230 //----/
231
232 // --- STEP 5
233 insert_edges_in_map(&mut cmap, &edges);
234 unsafe_time_section!(instant, timers::Section::BuildMeshInsertEdge);
235 //----/
236
237 unsafe_time_section!(kernel, timers::Section::BuildMeshTot);
238 //-------/
239
240 // --- CLIP
241 match clip {
242 Clip::Left => clip_left(&mut cmap)?,
243 Clip::Right => clip_right(&mut cmap)?,
244 Clip::None => {}
245 }
246 unsafe_time_section!(instant, timers::Section::Clip);
247 //----/
248
249 // CLEANUP
250 cmap.remove_attribute_storage::<Boundary>();
251 finish!(instant);
252 //-/
253
254 Ok(cmap)
255}
256
257#[cfg(test)]
258mod tests;