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;
honeycomb_kernels/grisubal/mod.rs

honeycomb_kernels/grisubal/
mod.rs
