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, false)?;
181    let [cx, cy] = grid_cell_sizes;
182    let ogrid = GridDescriptor::default()
183        .n_cells([nx, ny])
184        .len_per_cell([cx, cy])
185        .origin([origin.0, origin.1]);
186    unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid);
187    //----/
188
189    // --- REMOVE REDUNDANT PoIs
190    remove_redundant_poi(&mut geometry, grid_cell_sizes, origin);
191    unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi);
192    //----/
193
194    // ------ START MAIN KERNEL TIMER
195    start_timer!(kernel);
196
197    // --- BUILD THE GRID
198    let mut cmap = CMapBuilder::from_grid_descriptor(ogrid)
199        .add_attribute::<Boundary>() // will be used for clipping
200        .build()
201        .expect("E: unreachable"); // unreachable because grid dims are valid
202    unsafe_time_section!(instant, timers::Section::BuildMeshInit);
203    //----/
204
205    // process the geometry
206
207    // --- STEP 1 & 2
208    // (1)
209    let (new_segments, intersection_metadata) =
210        generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin);
211    // (2)
212    let n_intersec = intersection_metadata.len();
213    let (edge_intersec, dart_slices) =
214        group_intersections_per_edge(&mut cmap, intersection_metadata);
215    let intersection_darts = compute_intersection_ids(n_intersec, &edge_intersec, &dart_slices);
216    unsafe_time_section!(instant, timers::Section::BuildMeshIntersecData);
217    //----/
218
219    // --- STEP 3
220    insert_intersections(&cmap, &edge_intersec, &dart_slices);
221    unsafe_time_section!(instant, timers::Section::BuildMeshInsertIntersec);
222    //----/
223
224    // --- STEP 4
225    let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts);
226    unsafe_time_section!(instant, timers::Section::BuildMeshEdgeData);
227    //----/
228
229    // --- STEP 5
230    insert_edges_in_map(&mut cmap, &edges);
231    unsafe_time_section!(instant, timers::Section::BuildMeshInsertEdge);
232    //----/
233
234    unsafe_time_section!(kernel, timers::Section::BuildMeshTot);
235    //-------/
236
237    // --- CLIP
238    match clip {
239        Clip::Left => clip_left(&mut cmap)?,
240        Clip::Right => clip_right(&mut cmap)?,
241        Clip::None => {}
242    }
243    unsafe_time_section!(instant, timers::Section::Clip);
244    //----/
245
246    // CLEANUP
247    cmap.remove_attribute_storage::<Boundary>();
248    finish!(instant);
249    //-/
250
251    Ok(cmap)
252}
253
254#[cfg(test)]
255mod tests;
honeycomb_kernels/grisubal/mod.rs

honeycomb_kernels/grisubal/
mod.rs
