diff --git a/honeycomb-kernels/src/grisubal/grid.rs b/honeycomb-kernels/src/grisubal/grid.rs deleted file mode 100644 index b9ff73003..000000000 --- a/honeycomb-kernels/src/grisubal/grid.rs +++ /dev/null @@ -1,30 +0,0 @@ -//! Grid related code -//! -//! This module contains all code related to the usage of an overlapping - -// ------ IMPORTS - -// ------ CONTENT - -/// Structure used to index the overlapping grid's cases. -/// -/// Cells `(X, Y)` take value in range `(0, 0)` to `(N, M)`, -/// from left to right (X), from bottom to top (Y). -#[derive(PartialEq, Clone, Copy)] -pub struct GridCellId(pub usize, pub usize); - -impl GridCellId { - /// Compute the [Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry) between - /// two cells. - pub fn man_dist(lhs: &Self, rhs: &Self) -> usize { - lhs.0.abs_diff(rhs.0) + lhs.1.abs_diff(rhs.1) - } - - #[allow(clippy::cast_possible_wrap)] - pub fn diff(lhs: &Self, rhs: &Self) -> (isize, isize) { - ( - rhs.0 as isize - lhs.0 as isize, - rhs.1 as isize - lhs.1 as isize, - ) - } -} diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 26f2f6434..ea4036188 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -52,27 +52,46 @@ // ------ MODULE DECLARATIONS -pub(crate) mod clip; -pub(crate) mod grid; -pub(crate) mod kernel; pub(crate) mod model; - -// ------ RE-EXPORTS - -pub use model::Clip; +pub(crate) mod routines; +pub(crate) mod timers; // ------ IMPORTS -use crate::grisubal::clip::{clip_left, clip_right}; -use crate::grisubal::model::{ - compute_overlapping_grid, detect_orientation_issue, remove_redundant_poi, Boundary, Geometry2, +use crate::grisubal::{ + model::{Boundary, Geometry2}, + routines::{ + clip_left, clip_right, compute_intersection_ids, compute_overlapping_grid, + detect_orientation_issue, generate_edge_data, generate_intersection_data, + group_intersections_per_edge, insert_edges_in_map, insert_intersections, + remove_redundant_poi, + }, + timers::{finish, start_timer, unsafe_time_section}, +}; +use honeycomb_core::{ + cmap::{CMapBuilder, GridDescriptor}, + prelude::{CMap2, CoordsFloat}, }; -use honeycomb_core::prelude::{CMap2, CoordsFloat}; use thiserror::Error; use vtkio::Vtk; // ------ CONTENT +/// Post-processing clip operation. +/// +/// Note that the part of the map that is clipped depends on the orientation of the original geometry provided as +/// input. +#[derive(Default)] +pub enum Clip { + /// Clip elements located on the left side of the oriented boundary. + Left, + /// Clip elements located on the right side of the oriented boundary. + Right, + /// Keep all elements. Default value. + #[default] + None, +} + #[derive(Error, Debug)] /// Enum used to model potential errors of the `grisubal` kernel. /// @@ -93,38 +112,6 @@ pub enum GrisubalError { UnsupportedVtkData(&'static str), } -/// Global timers for execution times per-section. -#[cfg(feature = "profiling")] -static mut TIMERS: [Option; 13] = [None; 13]; - -/// Kernel section. -#[cfg(feature = "profiling")] -enum Section { - ImportVTK = 0, - BuildGeometry, - DetectOrientation, - ComputeOverlappingGrid, - RemoveRedundantPoi, - BuildMeshTot, - BuildMeshInit, - BuildMeshIntersecData, - BuildMeshInsertIntersec, - BuildMeshEdgeData, - BuildMeshInsertEdge, - Clip, - Cleanup, -} - -#[cfg(feature = "profiling")] -macro_rules! unsafe_time_section { - ($inst: ident, $sec: expr) => { - unsafe { - TIMERS[$sec as usize] = Some($inst.elapsed()); - $inst = std::time::Instant::now(); - } - }; -} - #[allow(clippy::missing_errors_doc)] /// Main algorithm call function. /// @@ -175,79 +162,101 @@ pub fn grisubal( grid_cell_sizes: [T; 2], clip: Clip, ) -> Result, GrisubalError> { - #[cfg(feature = "profiling")] - let mut instant = std::time::Instant::now(); + // INIT TIMER + start_timer!(instant); - // load geometry from file + // --- IMPORT VTK INPUT let geometry_vtk = match Vtk::import(file_path) { Ok(vtk) => vtk, Err(e) => panic!("E: could not open specified vtk file - {e}"), }; + unsafe_time_section!(instant, timers::Section::ImportVTK); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::ImportVTK); - - // pre-processing + // --- BUILD OUR MODEL FROM THE VTK IMPORT let mut geometry = Geometry2::try_from(geometry_vtk)?; + unsafe_time_section!(instant, timers::Section::BuildGeometry); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildGeometry); - + // --- FIRST DETECTION OF ORIENTATION ISSUES detect_orientation_issue(&geometry)?; + unsafe_time_section!(instant, timers::Section::DetectOrientation); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::DetectOrientation); + // --- FIND AN OVERLAPPING GRID + let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; + let [cx, cy] = grid_cell_sizes; + let ogrid = GridDescriptor::default() + .n_cells_x(nx) + .n_cells_y(ny) + .len_per_cell_x(cx) + .len_per_cell_y(cy) + .origin(origin); + unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid); + //----/ - // compute an overlapping grid & remove redundant PoIs - let (grid_n_cells, origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; + // --- REMOVE REDUNDANT PoIs + remove_redundant_poi(&mut geometry, grid_cell_sizes, origin); + unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::ComputeOverlappingGrid); + // ------ START MAIN KERNEL TIMER + start_timer!(kernel); - remove_redundant_poi(&mut geometry, grid_cell_sizes, origin); + // --- BUILD THE GRID + let mut cmap = CMapBuilder::default() + .grid_descriptor(ogrid) + .add_attribute::() // will be used for clipping + .build() + .expect("E: unreachable"); // unreachable because grid dims are valid + unsafe_time_section!(instant, timers::Section::BuildMeshInit); + //----/ + + // process the geometry - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::RemoveRedundantPoi); + // --- STEP 1 & 2 + // (1) + let (new_segments, intersection_metadata) = + generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin); + // (2) + let n_intersec = intersection_metadata.len(); + let (edge_intersec, dart_slices) = + group_intersections_per_edge(&mut cmap, intersection_metadata); + let intersection_darts = compute_intersection_ids(n_intersec, &edge_intersec, &dart_slices); + unsafe_time_section!(instant, timers::Section::BuildMeshIntersecData); + //----/ - // build the map - let mut cmap = kernel::build_mesh(&geometry, grid_cell_sizes, grid_n_cells, origin); + // --- STEP 3 + insert_intersections(&mut cmap, &edge_intersec, &dart_slices); + unsafe_time_section!(instant, timers::Section::BuildMeshInsertIntersec); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshTot); + // --- STEP 4 + let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts); + unsafe_time_section!(instant, timers::Section::BuildMeshEdgeData); + //----/ - // optional post-processing + // --- STEP 5 + insert_edges_in_map(&mut cmap, &edges); + unsafe_time_section!(instant, timers::Section::BuildMeshInsertEdge); + //----/ + + unsafe_time_section!(kernel, timers::Section::BuildMeshTot); + //-------/ + + // --- CLIP match clip { Clip::Left => clip_left(&mut cmap)?, Clip::Right => clip_right(&mut cmap)?, Clip::None => {} } + unsafe_time_section!(instant, timers::Section::Clip); + //----/ - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::Clip); - - // remove attribute used for clipping + // CLEANUP cmap.remove_attribute_storage::(); - - #[cfg(feature = "profiling")] - unsafe { - TIMERS[Section::Cleanup as usize] = Some(instant.elapsed()); - println!( - "{},{},{},{},{},{},{},{},{},{},{},{},{}", - TIMERS[0].unwrap().as_nanos(), - TIMERS[1].unwrap().as_nanos(), - TIMERS[2].unwrap().as_nanos(), - TIMERS[3].unwrap().as_nanos(), - TIMERS[4].unwrap().as_nanos(), - TIMERS[5].unwrap().as_nanos(), - TIMERS[6].unwrap().as_nanos(), - TIMERS[7].unwrap().as_nanos(), - TIMERS[8].unwrap().as_nanos(), - TIMERS[9].unwrap().as_nanos(), - TIMERS[10].unwrap().as_nanos(), - TIMERS[11].unwrap().as_nanos(), - TIMERS[12].unwrap().as_nanos(), - ); - } + finish!(instant); + //-/ Ok(cmap) } diff --git a/honeycomb-kernels/src/grisubal/model.rs b/honeycomb-kernels/src/grisubal/model.rs index c9fbb89e4..2ed10c95b 100644 --- a/honeycomb-kernels/src/grisubal/model.rs +++ b/honeycomb-kernels/src/grisubal/model.rs @@ -6,38 +6,41 @@ // ------ IMPORTS -use std::collections::HashSet; - +use crate::grisubal::GrisubalError; use honeycomb_core::attributes::AttrSparseVec; use honeycomb_core::prelude::{ AttributeBind, AttributeUpdate, CoordsFloat, DartIdentifier, OrbitPolicy, Vertex2, }; - use vtkio::{ model::{CellType, DataSet, VertexNumbers}, IOBuffer, Vtk, }; -use crate::grisubal::grid::GridCellId; -use crate::grisubal::GrisubalError; -#[cfg(doc)] -use honeycomb_core::prelude::CMap2; - // ------ CONTENT -/// Post-processing clip operation. +/// Structure used to index the overlapping grid's cases. /// -/// Note that the part of the map that is clipped depends on the orientation of the original geometry provided as -/// input. -#[derive(Default)] -pub enum Clip { - /// Clip elements located on the left side of the oriented boundary. - Left, - /// Clip elements located on the right side of the oriented boundary. - Right, - /// Keep all elements. Default value. - #[default] - None, +/// Cells `(X, Y)` take value in range `(0, 0)` to `(N, M)`, +/// from left to right (X), from bottom to top (Y). +#[derive(PartialEq, Clone, Copy)] +pub struct GridCellId(pub usize, pub usize); + +impl GridCellId { + /// Compute the [L1 / Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry) between + /// two cells. + pub fn l1_dist(lhs: &Self, rhs: &Self) -> usize { + lhs.0.abs_diff(rhs.0) + lhs.1.abs_diff(rhs.1) + } + + /// Compute the substraction between cell indices. This corresponds to an offset / movement over + /// the grid **from `lhs` to `rhs`**. + #[allow(clippy::cast_possible_wrap)] + pub fn offset(lhs: &Self, rhs: &Self) -> (isize, isize) { + ( + rhs.0 as isize - lhs.0 as isize, + rhs.1 as isize - lhs.1 as isize, + ) + } } /// Geometry representation structure. @@ -202,193 +205,6 @@ impl TryFrom for Geometry2 { } } -/// Check for orientation issue **per boundary**. -/// -/// This function check for the most obvious orientation issue; given a boundary, are all segments making it up -/// oriented consistently. If it is not the case, then there is at least one of: -/// -/// - a vertex being the origin of two segment -/// - a vertex being the end-point of two segment -/// -/// This does not cover consistent orientation across distinct boundaries (e.g. a geometry with a hole in it). -pub fn detect_orientation_issue( - geometry: &Geometry2, -) -> Result<(), GrisubalError> { - let mut origins = HashSet::new(); - let mut endpoints = HashSet::new(); - - for (orig, endp) in &geometry.segments { - if !origins.insert(orig) || !endpoints.insert(endp) { - return Err(GrisubalError::InconsistentOrientation( - "in-boundary inconsistency", - )); - } - } - - Ok(()) -} - -#[allow(clippy::cast_precision_loss)] -pub fn compute_overlapping_grid( - geometry: &Geometry2, - [len_cell_x, len_cell_y]: [T; 2], -) -> Result<([usize; 2], Vertex2), GrisubalError> { - // compute the minimum bounding box - let (mut min_x, mut max_x, mut min_y, mut max_y): (T, T, T, T) = { - let Some(tmp) = geometry.vertices.first() else { - return Err(GrisubalError::InvalidShape("no vertex in shape")); - }; - - (tmp.x(), tmp.x(), tmp.y(), tmp.y()) - }; - - geometry.vertices.iter().for_each(|v| { - min_x = min_x.min(v.x()); - max_x = max_x.max(v.x()); // may not be optimal - min_y = min_y.min(v.y()); // don't care - max_y = max_y.max(v.y()); - }); - - if max_x <= min_x { - return Err(GrisubalError::InvalidShape( - "bounding values along X axis are equal", - )); - } - if max_y <= min_y { - return Err(GrisubalError::InvalidShape( - "bounding values along Y axis are equal", - )); - } - - // compute characteristics of the overlapping Cartesian grid - - // create a ~one-and-a-half cell buffer to contain the geometry - // this, along with the `+1` below, guarantees that - // dart at the boundary of the grid are not intersected by the geometry - let mut og_x = min_x - len_cell_x * T::from(1.5).unwrap(); - let mut og_y = min_y - len_cell_y * T::from(1.5).unwrap(); - // we check for some extremely annoying cases here - // if some are detected, the origin is incrementally shifted - let (mut on_corner, mut reflect) = - detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y)); - let mut i = 1; - - while on_corner | reflect { - eprintln!( - "W: land on corner: {on_corner} - reflect on an axis: {reflect}, shifting origin" - ); - og_x += len_cell_x * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap(); - og_y += len_cell_y * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap(); - (on_corner, reflect) = - detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y)); - i += 1; - } - - let n_cells_x = ((max_x - og_x) / len_cell_x).ceil().to_usize().unwrap() + 1; - let n_cells_y = ((max_y - og_y) / len_cell_y).ceil().to_usize().unwrap() + 1; - - Ok(([n_cells_x, n_cells_y], Vertex2(og_x, og_y))) -} - -/// Remove from their geometry points of interest that intersect with a grid of specified dimension. -/// -/// This function works under the assumption that the grid is Cartesian & has its origin on `(0.0, 0.0)`. -pub fn remove_redundant_poi( - geometry: &mut Geometry2, - [cx, cy]: [T; 2], - origin: Vertex2, -) { - // PoI that land on the grid create a number of issues; removing them is ok since we're intersecting the grid - // at their coordinates, so the shape will be captured via intersection anyway - geometry.poi.retain(|idx| { - let v = geometry.vertices[*idx]; - // origin is assumed to be (0.0, 0.0) - let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); - let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); - !(on_x_axis | on_y_axis) - }); -} - -pub fn detect_overlaps( - geometry: &Geometry2, - [cx, cy]: [T; 2], - origin: Vertex2, -) -> (bool, bool) { - let on_corner = geometry - .vertices - .iter() - .map(|v| { - let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); - let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); - on_x_axis && on_y_axis - }) - .any(|a| a); - - let bad_reflection = geometry - .vertices - .iter() - .enumerate() - .filter_map(|(id, v)| { - let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); - let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); - if on_x_axis | on_y_axis { - return Some(id); - } - None - }) - // skip vertices that do not belong to the boundary - .filter(|id| { - geometry - .segments - .iter() - .any(|(v1, v2)| (id == v1) || (id == v2)) - }) - .map(|id| { - // if a vertex appear in the boundary, there should be both a segment landing and a - // segment starting on the vertex; hence `.expect()` - let vid_in = geometry - .segments - .iter() - .find_map(|(vin, ref_id)| { - if id == *ref_id { - return Some(*vin); - } - None - }) - .expect("E: found a vertex with no incident segment - is the geometry open?"); - // same - let vid_out = geometry - .segments - .iter() - .find_map(|(ref_id, vout)| { - if id == *ref_id { - return Some(*vout); - } - None - }) - .expect("E: found a vertex with no incident segment - is the geometry open?"); - let v_in = geometry.vertices[vid_in]; - let v_out = geometry.vertices[vid_out]; - let Vertex2(ox, oy) = origin; - let (c_in, c_out) = ( - GridCellId( - ((v_in.x() - ox) / cx).floor().to_usize().unwrap(), - ((v_in.y() - oy) / cy).floor().to_usize().unwrap(), - ), - GridCellId( - ((v_out.x() - ox) / cx).floor().to_usize().unwrap(), - ((v_out.y() - oy) / cy).floor().to_usize().unwrap(), - ), - ); - // if v_in and v_out belong to the same grid cell, there was a "reflection" on one - // of the grid's axis - c_in == c_out - }) - .any(|a| a); - - (on_corner, bad_reflection) -} - #[derive(Clone, Debug, PartialEq, Eq, Hash)] pub enum GeometryVertex { /// Regular vertex. Inner `usize` indicates the vertex ID in-geometry. diff --git a/honeycomb-kernels/src/grisubal/clip.rs b/honeycomb-kernels/src/grisubal/routines/clip.rs similarity index 99% rename from honeycomb-kernels/src/grisubal/clip.rs rename to honeycomb-kernels/src/grisubal/routines/clip.rs index 251751455..6b5c80ee9 100644 --- a/honeycomb-kernels/src/grisubal/clip.rs +++ b/honeycomb-kernels/src/grisubal/routines/clip.rs @@ -8,6 +8,7 @@ use honeycomb_core::prelude::{ CMap2, CoordsFloat, DartIdentifier, FaceIdentifier, Orbit2, OrbitPolicy, Vertex2, NULL_DART_ID, }; use std::collections::{HashSet, VecDeque}; + // ------ CONTENT /// Clip content on the left side of the boundary. diff --git a/honeycomb-kernels/src/grisubal/kernel.rs b/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs similarity index 56% rename from honeycomb-kernels/src/grisubal/kernel.rs rename to honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs index fbcbbe480..07496be96 100644 --- a/honeycomb-kernels/src/grisubal/kernel.rs +++ b/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs @@ -1,27 +1,19 @@ -//! Module short description +//! Step 1 implementation //! -//! Should you interact with this module directly? -//! -//! Content description if needed +//! The aim of this step is to build an exhaustive list of the segments making up +//! the geometry intersected with the grid: For each segment, if both vertices +//! do not belong to the same cell, we break it into sub-segments until it is the case. // ------ IMPORTS +use super::Segments; +use crate::grisubal::model::{Geometry2, GeometryVertex, GridCellId}; +use honeycomb_core::prelude::{CMap2, CoordsFloat, DartIdentifier, Vertex2, NULL_DART_ID}; use std::{ cmp::{max, min}, - collections::{HashMap, VecDeque}, -}; - -#[cfg(feature = "profiling")] -use super::{Section, TIMERS}; -use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex, MapEdge}; -use crate::{grisubal::grid::GridCellId, splits::splitn_edge_no_alloc}; -use honeycomb_core::prelude::{ - CMap2, CMapBuilder, CoordsFloat, DartIdentifier, EdgeIdentifier, GridDescriptor, Vertex2, - NULL_DART_ID, + collections::VecDeque, }; -// ------ CONTENT - macro_rules! make_geometry_vertex { ($g: ident, $vid: ident) => { if $g.poi.contains(&$vid) { @@ -60,125 +52,7 @@ macro_rules! up_intersec { }}; } -#[cfg(feature = "profiling")] -macro_rules! unsafe_time_section { - ($inst: ident, $sec: expr) => { - unsafe { - TIMERS[$sec as usize] = Some($inst.elapsed()); - $inst = std::time::Instant::now(); - } - }; -} - -/// Inner building routine. -/// -/// This function builds a combinatorial map from the described geometry. The returned -/// map is an adjusted grid that can be clipped in order to keep only part of the mesh. -/// See [`grisubal::Clip`] for more information. -/// -/// # Arguments -/// -/// - `geometry: &Geometry2` -- Description of the input geometry. -/// -/// ## Generics -/// -/// - `T: CoordsFloat` -- Floating point type used for coordinate representation. -pub fn build_mesh( - geometry: &Geometry2, - [cx, cy]: [T; 2], - [nx, ny]: [usize; 2], - origin: Vertex2, -) -> CMap2 { - #[cfg(feature = "profiling")] - let mut instant = std::time::Instant::now(); - - // compute grid characteristics - // build grid descriptor - let ogrid = GridDescriptor::default() - .n_cells_x(nx) - .n_cells_y(ny) - .len_per_cell_x(cx) - .len_per_cell_y(cy) - .origin(origin); - - // build initial map - let mut cmap = CMapBuilder::default() - .grid_descriptor(ogrid) - .add_attribute::() // will be used for clipping - .build() - .expect("E: unreachable"); // urneachable because grid dims are valid - - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshInit); - - // process the geometry - - // STEP 1 - // the aim of this step is to build an exhaustive list of the segments making up - // the GEOMETRY INTERSECTED WITH THE GRID, i.e. for each segment, if both vertices - // do not belong to the same cell, we break it into sub-segments until it is the case. - - let (new_segments, intersection_metadata) = - generate_intersection_data(&cmap, geometry, [nx, ny], [cx, cy], origin); - - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshIntersecData); - - // STEP 1.5 - // precompute stuff to - // - parallelize step 2 - // - make step 2 and step 3 independent from each other - - let n_intersec = intersection_metadata.len(); - let (edge_intersec, dart_slices) = - group_intersections_per_edge(&mut cmap, intersection_metadata); - let intersection_darts = compute_intersection_ids(n_intersec, &edge_intersec, &dart_slices); - - // STEP 2 - // insert the intersection vertices into the map & recover their encoding dart. The output Vec has consistent - // indexing with the input Vec, meaning that indices in GeometryVertex::Intersec instances are still valid. - - insert_intersections(&mut cmap, &edge_intersec, &dart_slices); - - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshInsertIntersec); - - // STEP 3 - // now that we have a list of "atomic" (non-dividable) segments, we can use it to build the actual segments that - // will be inserted into the map. Intersections serve as anchor points for the new segments while PoI make up - // "intermediate" points of segments. - - let edges = generate_edge_data(&cmap, geometry, &new_segments, &intersection_darts); - - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshEdgeData); - - // STEP 4 - // now that we have some segments that are directly defined between intersections, we can use some N-maps' - // properties to easily build the geometry into the map. - // This part relies heavily on "conventions"; the most important thing to note is that the darts in `MapEdge` - // instances are very precisely set, and can therefore be used to create all the new connectivities. - - insert_edges_in_map(&mut cmap, &edges); - - #[cfg(feature = "profiling")] - unsafe { - TIMERS[Section::BuildMeshInsertEdge as usize] = Some(instant.elapsed()); - } - - // return result - cmap -} - -// --- type aliases for clarity - -pub type Segments = HashMap; - -pub type IntersectionsPerEdge = HashMap>; - -pub type DartSlices = Vec>; - -// --- main kernels steps +// ------ CONTENT #[allow( clippy::too_many_lines, @@ -186,7 +60,7 @@ pub type DartSlices = Vec>; clippy::cast_possible_wrap, clippy::cast_sign_loss )] -pub(super) fn generate_intersection_data( +pub(crate) fn generate_intersection_data( cmap: &CMap2, geometry: &Geometry2, [nx, _ny]: [usize; 2], @@ -213,8 +87,8 @@ pub(super) fn generate_intersection_data( ), ); ( - GridCellId::man_dist(&c1, &c2), - GridCellId::diff(&c1, &c2), + GridCellId::l1_dist(&c1, &c2), + GridCellId::offset(&c1, &c2), v1, v2, v1_id, @@ -516,251 +390,3 @@ pub(super) fn generate_intersection_data( }).collect(); (new_segments, intersection_metadata) } - -pub(super) fn group_intersections_per_edge( - cmap: &mut CMap2, - intersection_metadata: Vec<(DartIdentifier, T)>, -) -> (IntersectionsPerEdge, DartSlices) { - // group intersection data per edge, and associate an ID to each - let mut edge_intersec: HashMap> = - HashMap::new(); - intersection_metadata - .into_iter() - .filter(|(_, t)| !t.is_nan()) - .enumerate() - .for_each(|(idx, (dart_id, mut t))| { - // classify intersections per edge_id & adjust t if needed - let edge_id = cmap.edge_id(dart_id); - // condition works in 2D because edges are 2 darts at most - if edge_id != dart_id { - t = T::one() - t; - } - if let Some(storage) = edge_intersec.get_mut(&edge_id) { - // not the first intersction with this given edge - storage.push((idx, t, dart_id)); - } else { - // first intersction with this given edge - edge_intersec.insert(edge_id, vec![(idx, t, dart_id)]); - } - }); - - // sort per t for later - for vs in edge_intersec.values_mut() { - // panic unreachable because t s.t. t == NaN have been filtered previously - vs.sort_by(|(_, t1, _), (_, t2, _)| t1.partial_cmp(t2).expect("E: unreachable")); - } - - // prealloc darts that will be used for vertex insertion - let n_darts_per_seg: Vec<_> = edge_intersec.values().map(|vs| 2 * vs.len()).collect(); - let n_tot: usize = n_darts_per_seg.iter().sum(); - let tmp = cmap.add_free_darts(n_tot) as usize; - // the prefix sum gives an offset that corresponds to the starting index of each slice, minus - // the location of the allocated dart block (given by `tmp`) - // end of the slice is deduced using these values and the number of darts the current seg needs - let prefix_sum: Vec = n_darts_per_seg - .iter() - .scan(0, |state, &n_d| { - *state += n_d; - Some(*state - n_d) // we want an offset, not the actual sum - }) - .collect(); - - #[allow(clippy::cast_possible_truncation)] - let dart_slices: Vec> = n_darts_per_seg - .iter() - .zip(prefix_sum.iter()) - .map(|(n_d, start)| { - ((tmp + start) as DartIdentifier..(tmp + start + n_d) as DartIdentifier) - .collect::>() - }) - .collect(); - - (edge_intersec, dart_slices) -} - -pub(super) fn insert_intersections( - cmap: &mut CMap2, - edge_intersec: &IntersectionsPerEdge, - dart_slices: &DartSlices, -) { - for ((edge_id, vs), new_darts) in edge_intersec.iter().zip(dart_slices.iter()) { - let _ = splitn_edge_no_alloc( - cmap, - *edge_id, - new_darts, - &vs.iter().map(|(_, t, _)| *t).collect::>(), - ); - } -} - -pub(super) fn compute_intersection_ids( - n_intersec: usize, - edge_intersec: &IntersectionsPerEdge, - dart_slices: &DartSlices, -) -> Vec { - let mut res = vec![NULL_DART_ID; n_intersec]; - for ((edge_id, vs), new_darts) in edge_intersec.iter().zip(dart_slices.iter()) { - // order should be consistent between collection because of the sort_by call - let hl = new_darts.len() / 2; // half-length; also equal to n_intermediate - let fh = &new_darts[..hl]; // first half; used for the side of edge id - let sh = &new_darts[hl..]; // second half; used for the opposite side - for (i, (id, _, old_dart_id)) in vs.iter().enumerate() { - // readjust according to intersection side - res[*id] = if *old_dart_id == *edge_id { - fh[i] - } else { - sh[hl - 1 - i] - }; - } - } - res -} - -pub(super) fn generate_edge_data( - cmap: &CMap2, - geometry: &Geometry2, - new_segments: &Segments, - intersection_darts: &[DartIdentifier], -) -> Vec> { - new_segments - .iter() - .filter(|(k, _)| { - matches!( - k, - GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(..) - ) - }) - .map(|(start, v)| { - let mut end = v; - let mut intermediates = Vec::new(); - // while we land on regular vertices, go to the next - while !matches!( - end, - GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(_) - ) { - match end { - GeometryVertex::PoI(vid) => { - // save the PoI as an intermediate & update end point - intermediates.push(geometry.vertices[*vid]); - end = &new_segments[end]; - } - GeometryVertex::Regular(_) => { - // skip; update end point - end = &new_segments[end]; - } - GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(..) => { - unreachable!() // outer while should prevent this from happening - } - } - } - - let d_start = match start { - GeometryVertex::Intersec(d_start_idx) => { - cmap.beta::<2>(intersection_darts[*d_start_idx]) - } - GeometryVertex::IntersecCorner(d_in) => { - cmap.beta::<2>(cmap.beta::<1>(cmap.beta::<2>(*d_in))) - } - _ => unreachable!(), // unreachable due to filter - }; - let d_end = match end { - GeometryVertex::Intersec(d_end_idx) => intersection_darts[*d_end_idx], - GeometryVertex::IntersecCorner(d_in) => *d_in, - _ => unreachable!(), // unreachable due to filter - }; - - // the data in this structure can be used to entirely deduce the new connections that should be made - // at STEP 3 - - MapEdge { - start: d_start, - intermediates, - end: d_end, - } - }) - .collect() -} - -pub(super) fn insert_edges_in_map(cmap: &mut CMap2, edges: &[MapEdge]) { - // FIXME: minimize allocs & redundant operations - // prealloc all darts needed - let n_darts_per_seg: Vec<_> = edges - .iter() - .map(|e| 2 + 2 * e.intermediates.len()) - .collect(); - let n_tot: usize = n_darts_per_seg.iter().sum(); - let tmp = cmap.add_free_darts(n_tot) as usize; - // the prefix sum gives an offset that corresponds to the starting index of each slice, minus - // the location of the allocated dart block (given by `tmp`) - // end of the slice is deduced using these values and the number of darts the current seg needs - let prefix_sum: Vec = n_darts_per_seg - .iter() - .scan(0, |state, &n_d| { - *state += n_d; - Some(*state - n_d) // we want an offset, not the actual sum - }) - .collect(); - #[allow(clippy::cast_possible_truncation)] - let dart_slices: Vec> = n_darts_per_seg - .iter() - .zip(prefix_sum.iter()) - .map(|(n_d, start)| { - ((tmp + start) as DartIdentifier..(tmp + start + n_d) as DartIdentifier) - .collect::>() - }) - .collect(); - - // insert new edges - for ( - MapEdge { - start, - intermediates, - end, - }, - dslice, - ) in edges.iter().zip(dart_slices.iter()) - { - // remove deprecated connectivities & save what data is necessary - let b1_start_old = cmap.beta::<1>(*start); - let b0_end_old = cmap.beta::<0>(*end); - cmap.one_unlink(*start); - cmap.one_unlink(b0_end_old); - - let &[d_new, b2_d_new] = &dslice[0..2] else { - unreachable!() - }; - cmap.two_link(d_new, b2_d_new); - - // rebuild - this is the final construct if there are no intermediates - cmap.one_link(*start, d_new); - cmap.one_link(b2_d_new, b1_start_old); - cmap.one_link(d_new, *end); - cmap.one_link(b0_end_old, b2_d_new); - - if !intermediates.is_empty() { - // create the topology components - let edge_id = cmap.edge_id(d_new); - let new_darts = &dslice[2..]; - let _ = splitn_edge_no_alloc( - cmap, - edge_id, - new_darts, - &vec![T::from(0.5).unwrap(); intermediates.len()], - ); - // replace placeholder vertices - let mut dart_id = cmap.beta::<1>(edge_id as DartIdentifier); - for v in intermediates { - let vid = cmap.vertex_id(dart_id); - let _ = cmap.replace_vertex(vid, *v); - dart_id = cmap.beta::<1>(dart_id); - } - } - - let mut d_boundary = cmap.beta::<1>(*start); - while d_boundary != *end { - cmap.set_attribute::(d_boundary, Boundary::Left); - cmap.set_attribute::(cmap.beta::<2>(d_boundary), Boundary::Right); - d_boundary = cmap.beta::<1>(d_boundary); - } - } -} diff --git a/honeycomb-kernels/src/grisubal/routines/compute_new_edges.rs b/honeycomb-kernels/src/grisubal/routines/compute_new_edges.rs new file mode 100644 index 000000000..e632a390d --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/compute_new_edges.rs @@ -0,0 +1,78 @@ +//! Step 4 implementation +//! +//! Rebuild information about the edge that will be inserted into the map. This is done by using +//! the list of "atomic" segments to search for connections between intersections, discarding +//! regular points and registering points of interests. + +// ------ IMPORTS + +use super::Segments; +use crate::grisubal::model::{Geometry2, GeometryVertex, MapEdge}; +use honeycomb_core::prelude::{CMap2, CoordsFloat, DartIdentifier}; + +// ------ CONTENT + +pub(crate) fn generate_edge_data( + cmap: &CMap2, + geometry: &Geometry2, + new_segments: &Segments, + intersection_darts: &[DartIdentifier], +) -> Vec> { + new_segments + .iter() + .filter(|(k, _)| { + matches!( + k, + GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(..) + ) + }) + .map(|(start, v)| { + let mut end = v; + let mut intermediates = Vec::new(); + // while we land on regular vertices, go to the next + while !matches!( + end, + GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(_) + ) { + match end { + GeometryVertex::PoI(vid) => { + // save the PoI as an intermediate & update end point + intermediates.push(geometry.vertices[*vid]); + end = &new_segments[end]; + } + GeometryVertex::Regular(_) => { + // skip; update end point + end = &new_segments[end]; + } + GeometryVertex::Intersec(_) | GeometryVertex::IntersecCorner(..) => { + unreachable!() // outer while should prevent this from happening + } + } + } + + let d_start = match start { + GeometryVertex::Intersec(d_start_idx) => { + cmap.beta::<2>(intersection_darts[*d_start_idx]) + } + GeometryVertex::IntersecCorner(d_in) => { + cmap.beta::<2>(cmap.beta::<1>(cmap.beta::<2>(*d_in))) + } + _ => unreachable!(), // unreachable due to filter + }; + let d_end = match end { + GeometryVertex::Intersec(d_end_idx) => intersection_darts[*d_end_idx], + GeometryVertex::IntersecCorner(d_in) => *d_in, + _ => unreachable!(), // unreachable due to filter + }; + + // the data in this structure can be used to entirely deduce the new connections that should be made + // at STEP 3 + + MapEdge { + start: d_start, + intermediates, + end: d_end, + } + }) + .collect() +} diff --git a/honeycomb-kernels/src/grisubal/routines/insert_intersecs.rs b/honeycomb-kernels/src/grisubal/routines/insert_intersecs.rs new file mode 100644 index 000000000..d937241d5 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/insert_intersecs.rs @@ -0,0 +1,26 @@ +//! Step 3 implementation +//! +//! Insert the intersections into the map. + +// ------ IMPORTS + +use super::{DartSlices, IntersectionsPerEdge}; +use crate::splits::splitn_edge_no_alloc; +use honeycomb_core::prelude::{CMap2, CoordsFloat}; + +// ------ CONTENT + +pub(crate) fn insert_intersections( + cmap: &mut CMap2, + edge_intersec: &IntersectionsPerEdge, + dart_slices: &DartSlices, +) { + for ((edge_id, vs), new_darts) in edge_intersec.iter().zip(dart_slices.iter()) { + let _ = splitn_edge_no_alloc( + cmap, + *edge_id, + new_darts, + &vs.iter().map(|(_, t, _)| *t).collect::>(), + ); + } +} diff --git a/honeycomb-kernels/src/grisubal/routines/insert_new_edges.rs b/honeycomb-kernels/src/grisubal/routines/insert_new_edges.rs new file mode 100644 index 000000000..6ffd3446a --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/insert_new_edges.rs @@ -0,0 +1,95 @@ +//! Step 5 implementation +// +//! Use the information computed at step 4 and insert all new edges into the map. + +// ------ IMPORTS + +use crate::grisubal::model::{Boundary, MapEdge}; +use crate::splits::splitn_edge_no_alloc; +use honeycomb_core::prelude::{CMap2, CoordsFloat, DartIdentifier}; + +// ------ CONTENT + +pub(crate) fn insert_edges_in_map(cmap: &mut CMap2, edges: &[MapEdge]) { + // FIXME: minimize allocs & redundant operations + // prealloc all darts needed + let n_darts_per_seg: Vec<_> = edges + .iter() + .map(|e| 2 + 2 * e.intermediates.len()) + .collect(); + let n_tot: usize = n_darts_per_seg.iter().sum(); + let tmp = cmap.add_free_darts(n_tot) as usize; + // the prefix sum gives an offset that corresponds to the starting index of each slice, minus + // the location of the allocated dart block (given by `tmp`) + // end of the slice is deduced using these values and the number of darts the current seg needs + let prefix_sum: Vec = n_darts_per_seg + .iter() + .scan(0, |state, &n_d| { + *state += n_d; + Some(*state - n_d) // we want an offset, not the actual sum + }) + .collect(); + #[allow(clippy::cast_possible_truncation)] + let dart_slices: Vec> = n_darts_per_seg + .iter() + .zip(prefix_sum.iter()) + .map(|(n_d, start)| { + ((tmp + start) as DartIdentifier..(tmp + start + n_d) as DartIdentifier) + .collect::>() + }) + .collect(); + + // insert new edges + for ( + MapEdge { + start, + intermediates, + end, + }, + dslice, + ) in edges.iter().zip(dart_slices.iter()) + { + // remove deprecated connectivities & save what data is necessary + let b1_start_old = cmap.beta::<1>(*start); + let b0_end_old = cmap.beta::<0>(*end); + cmap.one_unlink(*start); + cmap.one_unlink(b0_end_old); + + let &[d_new, b2_d_new] = &dslice[0..2] else { + unreachable!() + }; + cmap.two_link(d_new, b2_d_new); + + // rebuild - this is the final construct if there are no intermediates + cmap.one_link(*start, d_new); + cmap.one_link(b2_d_new, b1_start_old); + cmap.one_link(d_new, *end); + cmap.one_link(b0_end_old, b2_d_new); + + if !intermediates.is_empty() { + // create the topology components + let edge_id = cmap.edge_id(d_new); + let new_darts = &dslice[2..]; + let _ = splitn_edge_no_alloc( + cmap, + edge_id, + new_darts, + &vec![T::from(0.5).unwrap(); intermediates.len()], + ); + // replace placeholder vertices + let mut dart_id = cmap.beta::<1>(edge_id as DartIdentifier); + for v in intermediates { + let vid = cmap.vertex_id(dart_id); + let _ = cmap.replace_vertex(vid, *v); + dart_id = cmap.beta::<1>(dart_id); + } + } + + let mut d_boundary = cmap.beta::<1>(*start); + while d_boundary != *end { + cmap.set_attribute::(d_boundary, Boundary::Left); + cmap.set_attribute::(cmap.beta::<2>(d_boundary), Boundary::Right); + d_boundary = cmap.beta::<1>(d_boundary); + } + } +} diff --git a/honeycomb-kernels/src/grisubal/routines/mod.rs b/honeycomb-kernels/src/grisubal/routines/mod.rs new file mode 100644 index 000000000..5db8318d0 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/mod.rs @@ -0,0 +1,52 @@ +//! Module short description +//! +//! Should you interact with this module directly? +//! +//! Content description if needed + +// ------ MODULES + +mod clip; +mod compute_intersecs; +mod compute_new_edges; +mod insert_intersecs; +mod insert_new_edges; +mod pre_processing; +mod process_intersecs_data; + +// ------ RE-EXPORTS + +// step 0 +pub(crate) use pre_processing::*; + +// step 1 +pub(crate) use compute_intersecs::*; + +// step 2 +pub(crate) use process_intersecs_data::*; + +// step 3 +pub(crate) use insert_intersecs::*; + +// step 4 +pub(crate) use compute_new_edges::*; + +// step 5 +pub(crate) use insert_new_edges::*; + +// optional clipping routines +pub(crate) use clip::{clip_left, clip_right}; + +// ------ IMPORTS + +use crate::grisubal::model::GeometryVertex; +use honeycomb_core::prelude::{DartIdentifier, EdgeIdentifier}; +use std::collections::HashMap; + +// ------ CONTENT + +pub type Segments = HashMap; + +pub type IntersectionsPerEdge = HashMap>; + +pub type DartSlices = Vec>; diff --git a/honeycomb-kernels/src/grisubal/routines/pre_processing.rs b/honeycomb-kernels/src/grisubal/routines/pre_processing.rs new file mode 100644 index 000000000..55997dc80 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/pre_processing.rs @@ -0,0 +1,201 @@ +//! Step 0 implementation + +// ------ IMPORTS + +// ------ CONTENT + +use std::collections::HashSet; + +use honeycomb_core::prelude::{CoordsFloat, Vertex2}; + +use crate::grisubal::{ + model::{Geometry2, GridCellId}, + GrisubalError, +}; + +/// Check for orientation issue **per boundary**. +/// +/// This function check for the most obvious orientation issue; given a boundary, are all segments making it up +/// oriented consistently. If it is not the case, then there is at least one of: +/// +/// - a vertex being the origin of two segment +/// - a vertex being the end-point of two segment +/// +/// This does not cover consistent orientation across distinct boundaries (e.g. a geometry with a hole in it). +pub fn detect_orientation_issue( + geometry: &Geometry2, +) -> Result<(), GrisubalError> { + let mut origins = HashSet::new(); + let mut endpoints = HashSet::new(); + + for (orig, endp) in &geometry.segments { + if !origins.insert(orig) || !endpoints.insert(endp) { + return Err(GrisubalError::InconsistentOrientation( + "in-boundary inconsistency", + )); + } + } + + Ok(()) +} + +#[allow(clippy::cast_precision_loss)] +pub fn compute_overlapping_grid( + geometry: &Geometry2, + [len_cell_x, len_cell_y]: [T; 2], +) -> Result<([usize; 2], Vertex2), GrisubalError> { + // compute the minimum bounding box + let (mut min_x, mut max_x, mut min_y, mut max_y): (T, T, T, T) = { + let Some(tmp) = geometry.vertices.first() else { + return Err(GrisubalError::InvalidShape("no vertex in shape")); + }; + + (tmp.x(), tmp.x(), tmp.y(), tmp.y()) + }; + + geometry.vertices.iter().for_each(|v| { + min_x = min_x.min(v.x()); + max_x = max_x.max(v.x()); // may not be optimal + min_y = min_y.min(v.y()); // don't care + max_y = max_y.max(v.y()); + }); + + if max_x <= min_x { + return Err(GrisubalError::InvalidShape( + "bounding values along X axis are equal", + )); + } + if max_y <= min_y { + return Err(GrisubalError::InvalidShape( + "bounding values along Y axis are equal", + )); + } + + // compute characteristics of the overlapping Cartesian grid + + // create a ~one-and-a-half cell buffer to contain the geometry + // this, along with the `+1` below, guarantees that + // dart at the boundary of the grid are not intersected by the geometry + let mut og_x = min_x - len_cell_x * T::from(1.5).unwrap(); + let mut og_y = min_y - len_cell_y * T::from(1.5).unwrap(); + // we check for some extremely annoying cases here + // if some are detected, the origin is incrementally shifted + let (mut on_corner, mut reflect) = + detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y)); + let mut i = 1; + + while on_corner | reflect { + eprintln!( + "W: land on corner: {on_corner} - reflect on an axis: {reflect}, shifting origin" + ); + og_x += len_cell_x * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap(); + og_y += len_cell_y * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap(); + (on_corner, reflect) = + detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y)); + i += 1; + } + + let n_cells_x = ((max_x - og_x) / len_cell_x).ceil().to_usize().unwrap() + 1; + let n_cells_y = ((max_y - og_y) / len_cell_y).ceil().to_usize().unwrap() + 1; + + Ok(([n_cells_x, n_cells_y], Vertex2(og_x, og_y))) +} + +/// Remove from their geometry points of interest that intersect with a grid of specified dimension. +/// +/// This function works under the assumption that the grid is Cartesian & has its origin on `(0.0, 0.0)`. +pub fn remove_redundant_poi( + geometry: &mut Geometry2, + [cx, cy]: [T; 2], + origin: Vertex2, +) { + // PoI that land on the grid create a number of issues; removing them is ok since we're intersecting the grid + // at their coordinates, so the shape will be captured via intersection anyway + geometry.poi.retain(|idx| { + let v = geometry.vertices[*idx]; + // origin is assumed to be (0.0, 0.0) + let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); + let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); + !(on_x_axis | on_y_axis) + }); +} + +pub fn detect_overlaps( + geometry: &Geometry2, + [cx, cy]: [T; 2], + origin: Vertex2, +) -> (bool, bool) { + let on_corner = geometry + .vertices + .iter() + .map(|v| { + let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); + let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); + on_x_axis && on_y_axis + }) + .any(|a| a); + + let bad_reflection = geometry + .vertices + .iter() + .enumerate() + .filter_map(|(id, v)| { + let on_x_axis = ((v.x() - origin.x()) % cx).is_zero(); + let on_y_axis = ((v.y() - origin.y()) % cy).is_zero(); + if on_x_axis | on_y_axis { + return Some(id); + } + None + }) + // skip vertices that do not belong to the boundary + .filter(|id| { + geometry + .segments + .iter() + .any(|(v1, v2)| (id == v1) || (id == v2)) + }) + .map(|id| { + // if a vertex appear in the boundary, there should be both a segment landing and a + // segment starting on the vertex; hence `.expect()` + let vid_in = geometry + .segments + .iter() + .find_map(|(vin, ref_id)| { + if id == *ref_id { + return Some(*vin); + } + None + }) + .expect("E: found a vertex with no incident segment - is the geometry open?"); + // same + let vid_out = geometry + .segments + .iter() + .find_map(|(ref_id, vout)| { + if id == *ref_id { + return Some(*vout); + } + None + }) + .expect("E: found a vertex with no incident segment - is the geometry open?"); + let v_in = geometry.vertices[vid_in]; + let v_out = geometry.vertices[vid_out]; + let Vertex2(ox, oy) = origin; + let (c_in, c_out) = ( + GridCellId( + ((v_in.x() - ox) / cx).floor().to_usize().unwrap(), + ((v_in.y() - oy) / cy).floor().to_usize().unwrap(), + ), + GridCellId( + ((v_out.x() - ox) / cx).floor().to_usize().unwrap(), + ((v_out.y() - oy) / cy).floor().to_usize().unwrap(), + ), + ); + // if v_in and v_out belong to the same grid cell, there was a "reflection" on one + // of the grid's axis + c_in == c_out + }) + .any(|a| a); + + (on_corner, bad_reflection) +} diff --git a/honeycomb-kernels/src/grisubal/routines/process_intersecs_data.rs b/honeycomb-kernels/src/grisubal/routines/process_intersecs_data.rs new file mode 100644 index 000000000..5002766c1 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/routines/process_intersecs_data.rs @@ -0,0 +1,97 @@ +//! Step 2 implementation +//! +//! The main goal of this step is tp precompute information to: +//! - parallelize step 3 +//! - make step 3 and step 4 independent from each other + +// ------ IMPORTS + +use super::{DartSlices, IntersectionsPerEdge}; +use honeycomb_core::prelude::{CMap2, CoordsFloat, DartIdentifier, EdgeIdentifier, NULL_DART_ID}; +use std::collections::HashMap; + +// ------ CONTENT + +pub(crate) fn group_intersections_per_edge( + cmap: &mut CMap2, + intersection_metadata: Vec<(DartIdentifier, T)>, +) -> (IntersectionsPerEdge, DartSlices) { + // group intersection data per edge, and associate an ID to each + let mut edge_intersec: HashMap> = + HashMap::new(); + intersection_metadata + .into_iter() + .filter(|(_, t)| !t.is_nan()) + .enumerate() + .for_each(|(idx, (dart_id, mut t))| { + // classify intersections per edge_id & adjust t if needed + let edge_id = cmap.edge_id(dart_id); + // condition works in 2D because edges are 2 darts at most + if edge_id != dart_id { + t = T::one() - t; + } + if let Some(storage) = edge_intersec.get_mut(&edge_id) { + // not the first intersction with this given edge + storage.push((idx, t, dart_id)); + } else { + // first intersction with this given edge + edge_intersec.insert(edge_id, vec![(idx, t, dart_id)]); + } + }); + + // sort per t for later + for vs in edge_intersec.values_mut() { + // panic unreachable because t s.t. t == NaN have been filtered previously + vs.sort_by(|(_, t1, _), (_, t2, _)| t1.partial_cmp(t2).expect("E: unreachable")); + } + + // prealloc darts that will be used for vertex insertion + let n_darts_per_seg: Vec<_> = edge_intersec.values().map(|vs| 2 * vs.len()).collect(); + let n_tot: usize = n_darts_per_seg.iter().sum(); + let tmp = cmap.add_free_darts(n_tot) as usize; + // the prefix sum gives an offset that corresponds to the starting index of each slice, minus + // the location of the allocated dart block (given by `tmp`) + // end of the slice is deduced using these values and the number of darts the current seg needs + let prefix_sum: Vec = n_darts_per_seg + .iter() + .scan(0, |state, &n_d| { + *state += n_d; + Some(*state - n_d) // we want an offset, not the actual sum + }) + .collect(); + + #[allow(clippy::cast_possible_truncation)] + let dart_slices: Vec> = n_darts_per_seg + .iter() + .zip(prefix_sum.iter()) + .map(|(n_d, start)| { + ((tmp + start) as DartIdentifier..(tmp + start + n_d) as DartIdentifier) + .collect::>() + }) + .collect(); + + (edge_intersec, dart_slices) +} + +pub(crate) fn compute_intersection_ids( + n_intersec: usize, + edge_intersec: &IntersectionsPerEdge, + dart_slices: &DartSlices, +) -> Vec { + let mut res = vec![NULL_DART_ID; n_intersec]; + for ((edge_id, vs), new_darts) in edge_intersec.iter().zip(dart_slices.iter()) { + // order should be consistent between collection because of the sort_by call + let hl = new_darts.len() / 2; // half-length; also equal to n_intermediate + let fh = &new_darts[..hl]; // first half; used for the side of edge id + let sh = &new_darts[hl..]; // second half; used for the opposite side + for (i, (id, _, old_dart_id)) in vs.iter().enumerate() { + // readjust according to intersection side + res[*id] = if *old_dart_id == *edge_id { + fh[i] + } else { + sh[hl - 1 - i] + }; + } + } + res +} diff --git a/honeycomb-kernels/src/grisubal/tests.rs b/honeycomb-kernels/src/grisubal/tests.rs index 519b44da8..e94dbf198 100644 --- a/honeycomb-kernels/src/grisubal/tests.rs +++ b/honeycomb-kernels/src/grisubal/tests.rs @@ -1,10 +1,10 @@ // ------ IMPORTS -use crate::grisubal::kernel::{ +use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex}; +use crate::grisubal::routines::{ compute_intersection_ids, generate_edge_data, generate_intersection_data, group_intersections_per_edge, insert_edges_in_map, insert_intersections, }; -use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex}; use honeycomb_core::prelude::{CMapBuilder, GridDescriptor, Orbit2, OrbitPolicy, Vertex2}; use vtkio::Vtk; // ------ CONTENT diff --git a/honeycomb-kernels/src/grisubal/timers.rs b/honeycomb-kernels/src/grisubal/timers.rs new file mode 100644 index 000000000..f7068dc6b --- /dev/null +++ b/honeycomb-kernels/src/grisubal/timers.rs @@ -0,0 +1,74 @@ +//! Internal timers +//! +//! **This code is only compiled if the `profiling` feature is enabled.** + +#[cfg(feature = "profiling")] +/// Global timers for execution times per-section. +pub(crate) static mut TIMERS: [Option; 13] = [None; 13]; + +#[cfg(feature = "profiling")] +/// Kernel section. +pub(crate) enum Section { + ImportVTK = 0, + BuildGeometry, + DetectOrientation, + ComputeOverlappingGrid, + RemoveRedundantPoi, + BuildMeshTot, + BuildMeshInit, + BuildMeshIntersecData, + BuildMeshInsertIntersec, + BuildMeshEdgeData, + BuildMeshInsertEdge, + Clip, + Cleanup, +} + +macro_rules! start_timer { + ($inst: ident) => { + #[cfg(feature = "profiling")] + let mut $inst = std::time::Instant::now(); + }; +} + +pub(crate) use start_timer; + +macro_rules! unsafe_time_section { + ($inst: ident, $sec: expr) => { + #[allow(unused_assignments)] + #[cfg(feature = "profiling")] + unsafe { + timers::TIMERS[$sec as usize] = Some($inst.elapsed()); + $inst = std::time::Instant::now(); + } + }; +} + +pub(crate) use unsafe_time_section; + +macro_rules! finish { + ($inst: ident) => { + #[cfg(feature = "profiling")] + unsafe { + timers::TIMERS[timers::Section::Cleanup as usize] = Some($inst.elapsed()); + println!( + "{},{},{},{},{},{},{},{},{},{},{},{},{}", + timers::TIMERS[0].unwrap().as_nanos(), + timers::TIMERS[1].unwrap().as_nanos(), + timers::TIMERS[2].unwrap().as_nanos(), + timers::TIMERS[3].unwrap().as_nanos(), + timers::TIMERS[4].unwrap().as_nanos(), + timers::TIMERS[5].unwrap().as_nanos(), + timers::TIMERS[6].unwrap().as_nanos(), + timers::TIMERS[7].unwrap().as_nanos(), + timers::TIMERS[8].unwrap().as_nanos(), + timers::TIMERS[9].unwrap().as_nanos(), + timers::TIMERS[10].unwrap().as_nanos(), + timers::TIMERS[11].unwrap().as_nanos(), + timers::TIMERS[12].unwrap().as_nanos(), + ); + } + }; +} + +pub(crate) use finish;