From 0d94a93808bcd89c3936902e7bdad5df601bd4fd Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 08:37:50 +0100 Subject: [PATCH 01/12] replace file with module folder --- honeycomb-kernels/src/grisubal/{kernel.rs => kernel/mod.rs} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename honeycomb-kernels/src/grisubal/{kernel.rs => kernel/mod.rs} (100%) diff --git a/honeycomb-kernels/src/grisubal/kernel.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel.rs rename to honeycomb-kernels/src/grisubal/kernel/mod.rs From f97f43adc51fbd514a3490af365288e562892a25 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 08:55:34 +0100 Subject: [PATCH 02/12] move internal timers def to dedicated module --- honeycomb-kernels/src/grisubal/kernel/mod.rs | 24 ++++--------- honeycomb-kernels/src/grisubal/mod.rs | 37 +++----------------- honeycomb-kernels/src/grisubal/timers.rs | 37 ++++++++++++++++++++ 3 files changed, 49 insertions(+), 49 deletions(-) create mode 100644 honeycomb-kernels/src/grisubal/timers.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/mod.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs index fbcbbe480..3d4c6c270 100644 --- a/honeycomb-kernels/src/grisubal/kernel/mod.rs +++ b/honeycomb-kernels/src/grisubal/kernel/mod.rs @@ -6,19 +6,19 @@ // ------ IMPORTS -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, }; +use std::{ + cmp::{max, min}, + collections::{HashMap, VecDeque}, +}; + +#[cfg(feature = "profiling")] +use super::timers::{unsafe_time_section, Section, TIMERS}; // ------ CONTENT @@ -60,16 +60,6 @@ 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 diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 26f2f6434..de70f3680 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -56,6 +56,8 @@ pub(crate) mod clip; pub(crate) mod grid; pub(crate) mod kernel; pub(crate) mod model; +#[cfg(feature = "profiling")] +pub(crate) mod timers; // ------ RE-EXPORTS @@ -71,6 +73,9 @@ use honeycomb_core::prelude::{CMap2, CoordsFloat}; use thiserror::Error; use vtkio::Vtk; +#[cfg(feature = "profiling")] +use timers::{unsafe_time_section, Section, TIMERS}; + // ------ CONTENT #[derive(Error, Debug)] @@ -93,38 +98,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. /// diff --git a/honeycomb-kernels/src/grisubal/timers.rs b/honeycomb-kernels/src/grisubal/timers.rs new file mode 100644 index 000000000..c25786bad --- /dev/null +++ b/honeycomb-kernels/src/grisubal/timers.rs @@ -0,0 +1,37 @@ +//! Internal timers +//! +//! **This code is only compiled if the `profiling` feature is enabled.** + +/// Global timers for execution times per-section. +#[cfg(feature = "profiling")] +pub(crate) static mut TIMERS: [Option; 13] = [None; 13]; + +/// Kernel section. +#[cfg(feature = "profiling")] +pub(crate) 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(); + } + }; +} + +pub(crate) use unsafe_time_section; From 13e82b10824842e7fe7cb3de47764f69566201ad Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 09:10:12 +0100 Subject: [PATCH 03/12] remove the grid module --- honeycomb-kernels/src/grisubal/grid.rs | 30 --------------- honeycomb-kernels/src/grisubal/kernel/mod.rs | 4 +- honeycomb-kernels/src/grisubal/mod.rs | 18 +++++++-- honeycomb-kernels/src/grisubal/model.rs | 39 +++++++++++--------- 4 files changed, 39 insertions(+), 52 deletions(-) delete mode 100644 honeycomb-kernels/src/grisubal/grid.rs 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/kernel/mod.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs index 3d4c6c270..723b88aff 100644 --- a/honeycomb-kernels/src/grisubal/kernel/mod.rs +++ b/honeycomb-kernels/src/grisubal/kernel/mod.rs @@ -6,8 +6,8 @@ // ------ IMPORTS -use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex, MapEdge}; -use crate::{grisubal::grid::GridCellId, splits::splitn_edge_no_alloc}; +use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex, GridCellId, MapEdge}; +use crate::splits::splitn_edge_no_alloc; use honeycomb_core::prelude::{ CMap2, CMapBuilder, CoordsFloat, DartIdentifier, EdgeIdentifier, GridDescriptor, Vertex2, NULL_DART_ID, diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index de70f3680..9989c5391 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -53,7 +53,6 @@ // ------ MODULE DECLARATIONS pub(crate) mod clip; -pub(crate) mod grid; pub(crate) mod kernel; pub(crate) mod model; #[cfg(feature = "profiling")] @@ -61,8 +60,6 @@ pub(crate) mod timers; // ------ RE-EXPORTS -pub use model::Clip; - // ------ IMPORTS use crate::grisubal::clip::{clip_left, clip_right}; @@ -78,6 +75,21 @@ use timers::{unsafe_time_section, Section, TIMERS}; // ------ 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. /// diff --git a/honeycomb-kernels/src/grisubal/model.rs b/honeycomb-kernels/src/grisubal/model.rs index c9fbb89e4..267642abd 100644 --- a/honeycomb-kernels/src/grisubal/model.rs +++ b/honeycomb-kernels/src/grisubal/model.rs @@ -6,38 +6,43 @@ // ------ 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 std::collections::HashSet; 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 [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, + ) + } } /// Geometry representation structure. From d388abaecd1215eda343e0af92397cfd526c88f3 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 09:20:33 +0100 Subject: [PATCH 04/12] inline kernel steps --- honeycomb-kernels/src/grisubal/kernel/mod.rs | 100 ------------------- honeycomb-kernels/src/grisubal/mod.rs | 88 +++++++++++++++- 2 files changed, 84 insertions(+), 104 deletions(-) diff --git a/honeycomb-kernels/src/grisubal/kernel/mod.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs index 723b88aff..85cd6248b 100644 --- a/honeycomb-kernels/src/grisubal/kernel/mod.rs +++ b/honeycomb-kernels/src/grisubal/kernel/mod.rs @@ -60,106 +60,6 @@ macro_rules! up_intersec { }}; } -/// 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; diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 9989c5391..dacbf597e 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -66,7 +66,12 @@ use crate::grisubal::clip::{clip_left, clip_right}; use crate::grisubal::model::{ compute_overlapping_grid, detect_orientation_issue, remove_redundant_poi, Boundary, Geometry2, }; +use honeycomb_core::cmap::{CMapBuilder, GridDescriptor}; use honeycomb_core::prelude::{CMap2, CoordsFloat}; +use kernel::{ + compute_intersection_ids, generate_edge_data, generate_intersection_data, + group_intersections_per_edge, insert_edges_in_map, insert_intersections, +}; use thiserror::Error; use vtkio::Vtk; @@ -184,7 +189,8 @@ pub fn grisubal( unsafe_time_section!(instant, Section::DetectOrientation); // compute an overlapping grid & remove redundant PoIs - let (grid_n_cells, origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; + let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; + let [cx, cy] = grid_cell_sizes; #[cfg(feature = "profiling")] unsafe_time_section!(instant, Section::ComputeOverlappingGrid); @@ -194,11 +200,85 @@ pub fn grisubal( #[cfg(feature = "profiling")] unsafe_time_section!(instant, Section::RemoveRedundantPoi); - // build the map - let mut cmap = kernel::build_mesh(&geometry, grid_cell_sizes, grid_n_cells, origin); + // 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); + + #[cfg(feature = "profiling")] + let mut kernel = std::time::Instant::now(); + + // 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()); + } #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshTot); + unsafe_time_section!(kernel, Section::BuildMeshTot); // optional post-processing match clip { From 4b80d733c084b9da9831225ab4034477c4bc508f Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 11:17:20 +0100 Subject: [PATCH 05/12] move each kernel step into dedicated files --- honeycomb-kernels/src/grisubal/kernel/mod.rs | 657 +----------------- .../src/grisubal/kernel/step0.rs | 0 .../src/grisubal/kernel/step1.rs | 390 +++++++++++ .../src/grisubal/kernel/step2.rs | 93 +++ .../src/grisubal/kernel/step3.rs | 25 + .../src/grisubal/kernel/step4.rs | 74 ++ .../src/grisubal/kernel/step5.rs | 93 +++ honeycomb-kernels/src/grisubal/mod.rs | 8 +- 8 files changed, 697 insertions(+), 643 deletions(-) create mode 100644 honeycomb-kernels/src/grisubal/kernel/step0.rs create mode 100644 honeycomb-kernels/src/grisubal/kernel/step1.rs create mode 100644 honeycomb-kernels/src/grisubal/kernel/step2.rs create mode 100644 honeycomb-kernels/src/grisubal/kernel/step3.rs create mode 100644 honeycomb-kernels/src/grisubal/kernel/step4.rs create mode 100644 honeycomb-kernels/src/grisubal/kernel/step5.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/mod.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs index 85cd6248b..8eedd5fe7 100644 --- a/honeycomb-kernels/src/grisubal/kernel/mod.rs +++ b/honeycomb-kernels/src/grisubal/kernel/mod.rs @@ -4,653 +4,34 @@ //! //! Content description if needed -// ------ IMPORTS - -use crate::grisubal::model::{Boundary, Geometry2, GeometryVertex, GridCellId, MapEdge}; -use crate::splits::splitn_edge_no_alloc; -use honeycomb_core::prelude::{ - CMap2, CMapBuilder, CoordsFloat, DartIdentifier, EdgeIdentifier, GridDescriptor, Vertex2, - NULL_DART_ID, -}; -use std::{ - cmp::{max, min}, - collections::{HashMap, VecDeque}, -}; +// ------ MODULES -#[cfg(feature = "profiling")] -use super::timers::{unsafe_time_section, Section, TIMERS}; - -// ------ CONTENT +mod step0; +mod step1; +mod step2; +mod step3; +mod step4; +mod step5; -macro_rules! make_geometry_vertex { - ($g: ident, $vid: ident) => { - if $g.poi.contains(&$vid) { - GeometryVertex::PoI($vid) - } else { - GeometryVertex::Regular($vid) - } - }; -} +// ------ RE-EXPORTS -macro_rules! left_intersec { - ($va: ident, $vb: ident, $vdart:ident, $cy: ident) => {{ - let s = ($vdart.x() - $va.x()) / ($vb.x() - $va.x()); - (s, ($vdart.y() - $va.y() - ($vb.y() - $va.y()) * s) / $cy) - }}; -} +pub(crate) use step0::*; +pub(crate) use step1::*; +pub(crate) use step2::*; +pub(crate) use step3::*; +pub(crate) use step4::*; +pub(crate) use step5::*; -macro_rules! right_intersec { - ($va: ident, $vb: ident, $vdart:ident, $cy: ident) => {{ - let s = ($vdart.x() - $va.x()) / ($vb.x() - $va.x()); - (s, (($vb.y() - $va.y()) * s - ($vdart.y() - $va.y())) / $cy) - }}; -} - -macro_rules! down_intersec { - ($va: ident, $vb: ident, $vdart:ident, $cx: ident) => {{ - let s = ($vdart.y() - $va.y()) / ($vb.y() - $va.y()); - (s, (($vb.x() - $va.x()) * s - ($vdart.x() - $va.x())) / $cx) - }}; -} +// ------ IMPORTS -macro_rules! up_intersec { - ($va: ident, $vb: ident, $vdart:ident, $cx: ident) => {{ - let s = ($vdart.y() - $va.y()) / ($vb.y() - $va.y()); - (s, (($vdart.x() - $va.x()) - ($vb.x() - $va.x()) * s) / $cx) - }}; -} +use crate::grisubal::model::GeometryVertex; +use honeycomb_core::prelude::{DartIdentifier, EdgeIdentifier}; +use std::collections::HashMap; -// --- type aliases for clarity +// ------ CONTENT pub type Segments = HashMap; pub type IntersectionsPerEdge = HashMap>; pub type DartSlices = Vec>; - -// --- main kernels steps - -#[allow( - clippy::too_many_lines, - clippy::cast_possible_truncation, - clippy::cast_possible_wrap, - clippy::cast_sign_loss -)] -pub(super) fn generate_intersection_data( - cmap: &CMap2, - geometry: &Geometry2, - [nx, _ny]: [usize; 2], - [cx, cy]: [T; 2], - origin: Vertex2, -) -> (Segments, Vec<(DartIdentifier, T)>) { - let tmp: Vec<_> = geometry - .segments - .iter() - .map(|&(v1_id, v2_id)| { - // fetch vertices of the segment - let Vertex2(ox, oy) = origin; - let (v1, v2) = (&geometry.vertices[v1_id], &geometry.vertices[v2_id]); - // compute their position in the grid - // we assume that the origin of the grid is at (0., 0.) - let (c1, c2) = ( - GridCellId( - ((v1.x() - ox) / cx).floor().to_usize().unwrap(), - ((v1.y() - oy) / cy).floor().to_usize().unwrap(), - ), - GridCellId( - ((v2.x() - ox) / cx).floor().to_usize().unwrap(), - ((v2.y() - oy) / cy).floor().to_usize().unwrap(), - ), - ); - ( - GridCellId::man_dist(&c1, &c2), - GridCellId::diff(&c1, &c2), - v1, - v2, - v1_id, - v2_id, - c1, - ) - }) - .collect(); - // total number of intersection - let n_intersec: usize = tmp.iter().map(|(dist, _, _, _, _, _, _)| dist).sum(); - // we're using the prefix sum to compute an offset from the start. that's why we need a 0 at the front - // we'll cut off the last element later - let prefix_sum = tmp - .iter() - .map(|(dist, _, _, _, _, _, _)| dist) - .scan(0, |state, &dist| { - *state += dist; - Some(*state - dist) // we want an offset, not the actual sum - }); - // preallocate the intersection vector - let mut intersection_metadata = vec![(NULL_DART_ID, T::nan()); n_intersec]; - - let new_segments: Segments = tmp.iter().zip(prefix_sum).flat_map(|(&(dist, diff, v1, v2, v1_id, v2_id, c1), start)| { - let transform = Box::new(|seg: &[GeometryVertex]| { - assert_eq!(seg.len(), 2); - (seg[0].clone(), seg[1].clone()) - }); - // check neighbor status - match dist { - // trivial case: - // v1 & v2 belong to the same cell - 0 => { - vec![(make_geometry_vertex!(geometry, v1_id), make_geometry_vertex!(geometry, v2_id))] - } - // ok case: - // v1 & v2 belong to neighboring cells - 1 => { - // fetch base dart of the cell of v1 - #[allow(clippy::cast_possible_truncation)] - let d_base = (1 + 4 * c1.0 + nx * 4 * c1.1) as DartIdentifier; - // which dart does this correspond to? - #[rustfmt::skip] - let dart_id = match diff { - (-1, 0) => d_base + 3, - ( 1, 0) => d_base + 1, - ( 0, -1) => d_base , - ( 0, 1) => d_base + 2, - _ => unreachable!(), - }; - // what's the vertex associated to the dart? - let v_dart = cmap - .vertex(cmap.vertex_id(dart_id)) - .expect("E: found a topological vertex with no associated coordinates"); - // compute relative position of the intersection on the interecting edges - // `s` is relative to the segment `v1v2`, `t` to the grid's segment (the origin being `v_dart`) - #[rustfmt::skip] - let (_s, t) = match diff { - (-1, 0) => left_intersec!(v1, v2, v_dart, cy), - ( 1, 0) => right_intersec!(v1, v2, v_dart, cy), - ( 0, -1) => down_intersec!(v1, v2, v_dart, cx), - ( 0, 1) => up_intersec!(v1, v2, v_dart, cx), - _ => unreachable!(), - }; - - let id = start; - intersection_metadata[id] = (dart_id, t); - - vec![ - (make_geometry_vertex!(geometry, v1_id), GeometryVertex::Intersec(id)), - (GeometryVertex::Intersec(id), make_geometry_vertex!(geometry, v2_id)), - ] - } - // highly annoying case: - // v1 & v2 do not belong to neighboring cell - _ => { - // pure vertical / horizontal traversal are treated separately because it ensures we're not trying - // to compute intersections of parallel segments (which results at best in a division by 0) - let i_ids = start..start+dist; - match diff { - (i, 0) => { - // we can solve the intersection equation - // for each vertical edge of the grid we cross (i times) - let i_base = c1.0 as isize; - let tmp = - // the range is either - // i > 0: i_base..i_base + i - // or - // i < 0: i_base + 1 + i..i_base + 1 - (min(i_base, i_base + 1 + i)..max(i_base + i, i_base + 1)).zip(i_ids).map(|(x, id)| { - // cell base dart - let d_base = - (1 + 4 * x + (nx * 4 * c1.1) as isize) as DartIdentifier; - // intersected dart - let dart_id = if i.is_positive() { - d_base + 1 - } else { - d_base + 3 - }; - // vertex associated to the intersected dart - let v_dart = cmap.vertex(cmap.vertex_id(dart_id)) - .expect("E: found a topological vertex with no associated coordinates"); - // compute intersection - let (_s, t) = if i.is_positive() { - right_intersec!(v1, v2, v_dart, cy) - } else { - left_intersec!(v1, v2, v_dart, cy) - }; - - intersection_metadata[id] = (dart_id, t); - - GeometryVertex::Intersec(id) - }); - - // because of how the range is written, we need to reverse the iterator in one case - // to keep intersection ordered from v1 to v2 (i.e. ensure the segments we build are correct) - let mut vs: VecDeque = if i > 0 { - tmp.collect() - } else { - tmp.rev().collect() - }; - - // complete the vertex list - vs.push_front(make_geometry_vertex!(geometry, v1_id)); - vs.push_back(make_geometry_vertex!(geometry, v2_id)); - - vs.make_contiguous() - .windows(2) - .map(transform) - .collect::>() - } - (0, j) => { - // we can solve the intersection equation - // for each horizontal edge of the grid we cross (j times) - let j_base = c1.1 as isize; - let tmp = - // the range is either - // j > 0: j_base..j_base + j - // or - // j < 0: j_base + 1 + j..j_base + 1 - (min(j_base, j_base + 1 + j)..max(j_base + j, j_base + 1)).zip(i_ids).map(|(y, id)| { - // cell base dart - let d_base = (1 + 4 * c1.0 + nx * 4 * y as usize) as DartIdentifier; - // intersected dart - let dart_id = if j.is_positive() { d_base + 2 } else { d_base }; - // vertex associated to the intersected dart - let v_dart = cmap.vertex(cmap.vertex_id(dart_id)) - .expect("E: found a topological vertex with no associated coordinates"); - // compute intersection - let (_s, t) = if j.is_positive() { - up_intersec!(v1, v2, v_dart, cx) - } else { - down_intersec!(v1, v2, v_dart, cx) - }; - - intersection_metadata[id] = (dart_id, t); - - GeometryVertex::Intersec(id) - }); - - // because of how the range is written, we need to reverse the iterator in one case - // to keep intersection ordered from v1 to v2 (i.e. ensure the segments we build are correct) - let mut vs: VecDeque = if j > 0 { - tmp.collect() - } else { - tmp.rev().collect() - }; - - // complete the vertex list - vs.push_front(make_geometry_vertex!(geometry, v1_id)); - vs.push_back(make_geometry_vertex!(geometry, v2_id)); - - vs.make_contiguous() - .windows(2) - .map(transform) - .collect::>() - } - (i, j) => { - // in order to process this, we'll consider a "sub-grid" & use the direction of the segment to - // deduce which pair of dart we are supposed to intersect - // we also have to consider corner traversal; this corresponds to intersecting both darts of - // the pair at respective relative positions 1 and 0 (or 0 and 1) - let i_base = c1.0 as isize; - let j_base = c1.1 as isize; - let i_cell_range = min(i_base, i_base + i)..=max(i_base + i, i_base); - let j_cell_range = min(j_base, j_base + j)..=max(j_base + j, j_base); - let subgrid_cells = - i_cell_range.flat_map(|x| j_cell_range.clone().map(move |y| (x, y))); - - let mut intersec_data: Vec<(T, T, DartIdentifier)> = subgrid_cells - .map(|(x, y)| { - // cell base dart - let d_base = (1 + 4 * x + nx as isize * 4 * y) as DartIdentifier; - // (potentially) intersected darts - let vdart_id = if i.is_positive() { - d_base + 1 - } else { - d_base + 3 - }; - let hdart_id = if j.is_positive() { d_base + 2 } else { d_base }; - // associated vertices - let v_vdart = cmap.vertex(cmap.vertex_id(vdart_id)) - .expect("E: found a topological vertex with no associated coordinates"); - let v_hdart = cmap.vertex(cmap.vertex_id(hdart_id)) - .expect("E: found a topological vertex with no associated coordinates"); - // compute (potential) intersections - let v_coeffs = if i.is_positive() { - right_intersec!(v1, v2, v_vdart, cy) - } else { - left_intersec!(v1, v2, v_vdart, cy) - }; - let h_coeffs = if j.is_positive() { - up_intersec!(v1, v2, v_hdart, cx) - } else { - down_intersec!(v1, v2, v_hdart, cx) - }; - - (hdart_id, vdart_id, v_coeffs, h_coeffs) - }) - .filter_map(|(hdart_id, vdart_id, (vs, vt), (hs, ht))| { - let zero = T::zero(); - let one = T::one(); - // there is one corner intersection to check per (i, j) quadrant - match (i.is_positive(), j.is_positive()) { - // check - (true, true) | (false, false) => { - if ((vt - one).abs() < T::epsilon()) - && (ht.abs() < T::epsilon()) - { - return Some((hs, zero, hdart_id)); - } - } - (false, true) | (true, false) => { - if (vt.abs() < T::epsilon()) - && ((ht - one).abs() < T::epsilon()) - { - return Some((vs, zero, vdart_id)); - } - } - } - - // we can deduce if and which side is intersected using s and t values - // these should be comprised strictly between 0 and 1 for regular intersections - if (T::epsilon() <= vs) - & (vs <= one - T::epsilon()) - & (T::epsilon() <= vt) - & (vt <= one - T::epsilon()) - { - return Some((vs, vt, vdart_id)); // intersect vertical side - } - if (T::epsilon() < hs) - & (hs <= one - T::epsilon()) - & (T::epsilon() <= ht) - & (ht <= one - T::epsilon()) - { - return Some((hs, ht, hdart_id)); // intersect horizontal side - } - - // intersect none; this is possible since we're looking at cells of a subgrid, - // not following through the segment's intersections - None - }) - .collect(); - - // sort intersections from v1 to v2 - intersec_data.retain(|(s, _, _)| (T::zero() <= *s) && (*s <= T::one())); - // panic unreachable because of the retain above; there's no s s.t. s == NaN - intersec_data.sort_by(|(s1, _, _), (s2, _, _)| s1.partial_cmp(s2) - .expect("E: unreachable")); - - // collect geometry vertices - let mut vs = vec![make_geometry_vertex!(geometry, v1_id)]; - vs.extend(intersec_data.iter_mut().zip(i_ids).map(|((_, t, dart_id), id)| { - if t.is_zero() { - // we assume that the segment fully goes through the corner and does not land exactly - // on it, this allows us to compute directly the dart from which the next segment - // should start: the one incident to the vertex in the opposite quadrant - - // in that case, the preallocated intersection metadata slot will stay as (0, Nan) - // this is ok, we can simply ignore the entry when processing the data later - - let dart_in = *dart_id; - GeometryVertex::IntersecCorner(dart_in) - } else { - intersection_metadata[id] = (*dart_id, *t); - - GeometryVertex::Intersec(id) - } - })); - - vs.push(make_geometry_vertex!(geometry, v2_id)); - - vs.windows(2) - .map(transform) - .collect::>() - } - } - } - } - }).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/kernel/step0.rs b/honeycomb-kernels/src/grisubal/kernel/step0.rs new file mode 100644 index 000000000..e69de29bb diff --git a/honeycomb-kernels/src/grisubal/kernel/step1.rs b/honeycomb-kernels/src/grisubal/kernel/step1.rs new file mode 100644 index 000000000..6d6eacca1 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/kernel/step1.rs @@ -0,0 +1,390 @@ +//! Step 1 implementation +//! +//! + +// ------ 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::VecDeque, +}; + +macro_rules! make_geometry_vertex { + ($g: ident, $vid: ident) => { + if $g.poi.contains(&$vid) { + GeometryVertex::PoI($vid) + } else { + GeometryVertex::Regular($vid) + } + }; +} + +macro_rules! left_intersec { + ($va: ident, $vb: ident, $vdart:ident, $cy: ident) => {{ + let s = ($vdart.x() - $va.x()) / ($vb.x() - $va.x()); + (s, ($vdart.y() - $va.y() - ($vb.y() - $va.y()) * s) / $cy) + }}; +} + +macro_rules! right_intersec { + ($va: ident, $vb: ident, $vdart:ident, $cy: ident) => {{ + let s = ($vdart.x() - $va.x()) / ($vb.x() - $va.x()); + (s, (($vb.y() - $va.y()) * s - ($vdart.y() - $va.y())) / $cy) + }}; +} + +macro_rules! down_intersec { + ($va: ident, $vb: ident, $vdart:ident, $cx: ident) => {{ + let s = ($vdart.y() - $va.y()) / ($vb.y() - $va.y()); + (s, (($vb.x() - $va.x()) * s - ($vdart.x() - $va.x())) / $cx) + }}; +} + +macro_rules! up_intersec { + ($va: ident, $vb: ident, $vdart:ident, $cx: ident) => {{ + let s = ($vdart.y() - $va.y()) / ($vb.y() - $va.y()); + (s, (($vdart.x() - $va.x()) - ($vb.x() - $va.x()) * s) / $cx) + }}; +} + +// ------ CONTENT + +#[allow( + clippy::too_many_lines, + clippy::cast_possible_truncation, + clippy::cast_possible_wrap, + clippy::cast_sign_loss +)] +pub(crate) fn generate_intersection_data( + cmap: &CMap2, + geometry: &Geometry2, + [nx, _ny]: [usize; 2], + [cx, cy]: [T; 2], + origin: Vertex2, +) -> (Segments, Vec<(DartIdentifier, T)>) { + let tmp: Vec<_> = geometry + .segments + .iter() + .map(|&(v1_id, v2_id)| { + // fetch vertices of the segment + let Vertex2(ox, oy) = origin; + let (v1, v2) = (&geometry.vertices[v1_id], &geometry.vertices[v2_id]); + // compute their position in the grid + // we assume that the origin of the grid is at (0., 0.) + let (c1, c2) = ( + GridCellId( + ((v1.x() - ox) / cx).floor().to_usize().unwrap(), + ((v1.y() - oy) / cy).floor().to_usize().unwrap(), + ), + GridCellId( + ((v2.x() - ox) / cx).floor().to_usize().unwrap(), + ((v2.y() - oy) / cy).floor().to_usize().unwrap(), + ), + ); + ( + GridCellId::man_dist(&c1, &c2), + GridCellId::diff(&c1, &c2), + v1, + v2, + v1_id, + v2_id, + c1, + ) + }) + .collect(); + // total number of intersection + let n_intersec: usize = tmp.iter().map(|(dist, _, _, _, _, _, _)| dist).sum(); + // we're using the prefix sum to compute an offset from the start. that's why we need a 0 at the front + // we'll cut off the last element later + let prefix_sum = tmp + .iter() + .map(|(dist, _, _, _, _, _, _)| dist) + .scan(0, |state, &dist| { + *state += dist; + Some(*state - dist) // we want an offset, not the actual sum + }); + // preallocate the intersection vector + let mut intersection_metadata = vec![(NULL_DART_ID, T::nan()); n_intersec]; + + let new_segments: Segments = tmp.iter().zip(prefix_sum).flat_map(|(&(dist, diff, v1, v2, v1_id, v2_id, c1), start)| { + let transform = Box::new(|seg: &[GeometryVertex]| { + assert_eq!(seg.len(), 2); + (seg[0].clone(), seg[1].clone()) + }); + // check neighbor status + match dist { + // trivial case: + // v1 & v2 belong to the same cell + 0 => { + vec![(make_geometry_vertex!(geometry, v1_id), make_geometry_vertex!(geometry, v2_id))] + } + // ok case: + // v1 & v2 belong to neighboring cells + 1 => { + // fetch base dart of the cell of v1 + #[allow(clippy::cast_possible_truncation)] + let d_base = (1 + 4 * c1.0 + nx * 4 * c1.1) as DartIdentifier; + // which dart does this correspond to? + #[rustfmt::skip] + let dart_id = match diff { + (-1, 0) => d_base + 3, + ( 1, 0) => d_base + 1, + ( 0, -1) => d_base , + ( 0, 1) => d_base + 2, + _ => unreachable!(), + }; + // what's the vertex associated to the dart? + let v_dart = cmap + .vertex(cmap.vertex_id(dart_id)) + .expect("E: found a topological vertex with no associated coordinates"); + // compute relative position of the intersection on the interecting edges + // `s` is relative to the segment `v1v2`, `t` to the grid's segment (the origin being `v_dart`) + #[rustfmt::skip] + let (_s, t) = match diff { + (-1, 0) => left_intersec!(v1, v2, v_dart, cy), + ( 1, 0) => right_intersec!(v1, v2, v_dart, cy), + ( 0, -1) => down_intersec!(v1, v2, v_dart, cx), + ( 0, 1) => up_intersec!(v1, v2, v_dart, cx), + _ => unreachable!(), + }; + + let id = start; + intersection_metadata[id] = (dart_id, t); + + vec![ + (make_geometry_vertex!(geometry, v1_id), GeometryVertex::Intersec(id)), + (GeometryVertex::Intersec(id), make_geometry_vertex!(geometry, v2_id)), + ] + } + // highly annoying case: + // v1 & v2 do not belong to neighboring cell + _ => { + // pure vertical / horizontal traversal are treated separately because it ensures we're not trying + // to compute intersections of parallel segments (which results at best in a division by 0) + let i_ids = start..start+dist; + match diff { + (i, 0) => { + // we can solve the intersection equation + // for each vertical edge of the grid we cross (i times) + let i_base = c1.0 as isize; + let tmp = + // the range is either + // i > 0: i_base..i_base + i + // or + // i < 0: i_base + 1 + i..i_base + 1 + (min(i_base, i_base + 1 + i)..max(i_base + i, i_base + 1)).zip(i_ids).map(|(x, id)| { + // cell base dart + let d_base = + (1 + 4 * x + (nx * 4 * c1.1) as isize) as DartIdentifier; + // intersected dart + let dart_id = if i.is_positive() { + d_base + 1 + } else { + d_base + 3 + }; + // vertex associated to the intersected dart + let v_dart = cmap.vertex(cmap.vertex_id(dart_id)) + .expect("E: found a topological vertex with no associated coordinates"); + // compute intersection + let (_s, t) = if i.is_positive() { + right_intersec!(v1, v2, v_dart, cy) + } else { + left_intersec!(v1, v2, v_dart, cy) + }; + + intersection_metadata[id] = (dart_id, t); + + GeometryVertex::Intersec(id) + }); + + // because of how the range is written, we need to reverse the iterator in one case + // to keep intersection ordered from v1 to v2 (i.e. ensure the segments we build are correct) + let mut vs: VecDeque = if i > 0 { + tmp.collect() + } else { + tmp.rev().collect() + }; + + // complete the vertex list + vs.push_front(make_geometry_vertex!(geometry, v1_id)); + vs.push_back(make_geometry_vertex!(geometry, v2_id)); + + vs.make_contiguous() + .windows(2) + .map(transform) + .collect::>() + } + (0, j) => { + // we can solve the intersection equation + // for each horizontal edge of the grid we cross (j times) + let j_base = c1.1 as isize; + let tmp = + // the range is either + // j > 0: j_base..j_base + j + // or + // j < 0: j_base + 1 + j..j_base + 1 + (min(j_base, j_base + 1 + j)..max(j_base + j, j_base + 1)).zip(i_ids).map(|(y, id)| { + // cell base dart + let d_base = (1 + 4 * c1.0 + nx * 4 * y as usize) as DartIdentifier; + // intersected dart + let dart_id = if j.is_positive() { d_base + 2 } else { d_base }; + // vertex associated to the intersected dart + let v_dart = cmap.vertex(cmap.vertex_id(dart_id)) + .expect("E: found a topological vertex with no associated coordinates"); + // compute intersection + let (_s, t) = if j.is_positive() { + up_intersec!(v1, v2, v_dart, cx) + } else { + down_intersec!(v1, v2, v_dart, cx) + }; + + intersection_metadata[id] = (dart_id, t); + + GeometryVertex::Intersec(id) + }); + + // because of how the range is written, we need to reverse the iterator in one case + // to keep intersection ordered from v1 to v2 (i.e. ensure the segments we build are correct) + let mut vs: VecDeque = if j > 0 { + tmp.collect() + } else { + tmp.rev().collect() + }; + + // complete the vertex list + vs.push_front(make_geometry_vertex!(geometry, v1_id)); + vs.push_back(make_geometry_vertex!(geometry, v2_id)); + + vs.make_contiguous() + .windows(2) + .map(transform) + .collect::>() + } + (i, j) => { + // in order to process this, we'll consider a "sub-grid" & use the direction of the segment to + // deduce which pair of dart we are supposed to intersect + // we also have to consider corner traversal; this corresponds to intersecting both darts of + // the pair at respective relative positions 1 and 0 (or 0 and 1) + let i_base = c1.0 as isize; + let j_base = c1.1 as isize; + let i_cell_range = min(i_base, i_base + i)..=max(i_base + i, i_base); + let j_cell_range = min(j_base, j_base + j)..=max(j_base + j, j_base); + let subgrid_cells = + i_cell_range.flat_map(|x| j_cell_range.clone().map(move |y| (x, y))); + + let mut intersec_data: Vec<(T, T, DartIdentifier)> = subgrid_cells + .map(|(x, y)| { + // cell base dart + let d_base = (1 + 4 * x + nx as isize * 4 * y) as DartIdentifier; + // (potentially) intersected darts + let vdart_id = if i.is_positive() { + d_base + 1 + } else { + d_base + 3 + }; + let hdart_id = if j.is_positive() { d_base + 2 } else { d_base }; + // associated vertices + let v_vdart = cmap.vertex(cmap.vertex_id(vdart_id)) + .expect("E: found a topological vertex with no associated coordinates"); + let v_hdart = cmap.vertex(cmap.vertex_id(hdart_id)) + .expect("E: found a topological vertex with no associated coordinates"); + // compute (potential) intersections + let v_coeffs = if i.is_positive() { + right_intersec!(v1, v2, v_vdart, cy) + } else { + left_intersec!(v1, v2, v_vdart, cy) + }; + let h_coeffs = if j.is_positive() { + up_intersec!(v1, v2, v_hdart, cx) + } else { + down_intersec!(v1, v2, v_hdart, cx) + }; + + (hdart_id, vdart_id, v_coeffs, h_coeffs) + }) + .filter_map(|(hdart_id, vdart_id, (vs, vt), (hs, ht))| { + let zero = T::zero(); + let one = T::one(); + // there is one corner intersection to check per (i, j) quadrant + match (i.is_positive(), j.is_positive()) { + // check + (true, true) | (false, false) => { + if ((vt - one).abs() < T::epsilon()) + && (ht.abs() < T::epsilon()) + { + return Some((hs, zero, hdart_id)); + } + } + (false, true) | (true, false) => { + if (vt.abs() < T::epsilon()) + && ((ht - one).abs() < T::epsilon()) + { + return Some((vs, zero, vdart_id)); + } + } + } + + // we can deduce if and which side is intersected using s and t values + // these should be comprised strictly between 0 and 1 for regular intersections + if (T::epsilon() <= vs) + & (vs <= one - T::epsilon()) + & (T::epsilon() <= vt) + & (vt <= one - T::epsilon()) + { + return Some((vs, vt, vdart_id)); // intersect vertical side + } + if (T::epsilon() < hs) + & (hs <= one - T::epsilon()) + & (T::epsilon() <= ht) + & (ht <= one - T::epsilon()) + { + return Some((hs, ht, hdart_id)); // intersect horizontal side + } + + // intersect none; this is possible since we're looking at cells of a subgrid, + // not following through the segment's intersections + None + }) + .collect(); + + // sort intersections from v1 to v2 + intersec_data.retain(|(s, _, _)| (T::zero() <= *s) && (*s <= T::one())); + // panic unreachable because of the retain above; there's no s s.t. s == NaN + intersec_data.sort_by(|(s1, _, _), (s2, _, _)| s1.partial_cmp(s2) + .expect("E: unreachable")); + + // collect geometry vertices + let mut vs = vec![make_geometry_vertex!(geometry, v1_id)]; + vs.extend(intersec_data.iter_mut().zip(i_ids).map(|((_, t, dart_id), id)| { + if t.is_zero() { + // we assume that the segment fully goes through the corner and does not land exactly + // on it, this allows us to compute directly the dart from which the next segment + // should start: the one incident to the vertex in the opposite quadrant + + // in that case, the preallocated intersection metadata slot will stay as (0, Nan) + // this is ok, we can simply ignore the entry when processing the data later + + let dart_in = *dart_id; + GeometryVertex::IntersecCorner(dart_in) + } else { + intersection_metadata[id] = (*dart_id, *t); + + GeometryVertex::Intersec(id) + } + })); + + vs.push(make_geometry_vertex!(geometry, v2_id)); + + vs.windows(2) + .map(transform) + .collect::>() + } + } + } + } + }).collect(); + (new_segments, intersection_metadata) +} diff --git a/honeycomb-kernels/src/grisubal/kernel/step2.rs b/honeycomb-kernels/src/grisubal/kernel/step2.rs new file mode 100644 index 000000000..24f2ddb72 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/kernel/step2.rs @@ -0,0 +1,93 @@ +//! Step 2 implementation + +// ------ 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/kernel/step3.rs b/honeycomb-kernels/src/grisubal/kernel/step3.rs new file mode 100644 index 000000000..737c6ea3c --- /dev/null +++ b/honeycomb-kernels/src/grisubal/kernel/step3.rs @@ -0,0 +1,25 @@ +//! Step 3 implementation +//! + +// ------ 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/kernel/step4.rs b/honeycomb-kernels/src/grisubal/kernel/step4.rs new file mode 100644 index 000000000..6499ddb26 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/kernel/step4.rs @@ -0,0 +1,74 @@ +//! Step 4 implementation + +// ------ 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/kernel/step5.rs b/honeycomb-kernels/src/grisubal/kernel/step5.rs new file mode 100644 index 000000000..9eaea1131 --- /dev/null +++ b/honeycomb-kernels/src/grisubal/kernel/step5.rs @@ -0,0 +1,93 @@ +//! Step 5 implementation + +// ------ 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/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index dacbf597e..53a103559 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -273,13 +273,11 @@ pub fn grisubal( insert_edges_in_map(&mut cmap, &edges); #[cfg(feature = "profiling")] - unsafe { - TIMERS[Section::BuildMeshInsertEdge as usize] = Some(instant.elapsed()); + { + unsafe_time_section!(instant, Section::BuildMeshInsertEdge); + unsafe_time_section!(kernel, Section::BuildMeshTot); } - #[cfg(feature = "profiling")] - unsafe_time_section!(kernel, Section::BuildMeshTot); - // optional post-processing match clip { Clip::Left => clip_left(&mut cmap)?, From 4a1449dc0be65103168a38f837ed7ee0fa90ff29 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 11:27:52 +0100 Subject: [PATCH 06/12] move preproc function to step0.rs --- .../src/grisubal/kernel/step0.rs | 201 ++++++++++++++++++ honeycomb-kernels/src/grisubal/mod.rs | 9 +- honeycomb-kernels/src/grisubal/model.rs | 191 ----------------- 3 files changed, 205 insertions(+), 196 deletions(-) diff --git a/honeycomb-kernels/src/grisubal/kernel/step0.rs b/honeycomb-kernels/src/grisubal/kernel/step0.rs index e69de29bb..55997dc80 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step0.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step0.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/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 53a103559..c0eb4a815 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -63,14 +63,13 @@ 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}; use honeycomb_core::cmap::{CMapBuilder, GridDescriptor}; use honeycomb_core::prelude::{CMap2, CoordsFloat}; use kernel::{ - compute_intersection_ids, generate_edge_data, generate_intersection_data, - group_intersections_per_edge, insert_edges_in_map, insert_intersections, + 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, }; use thiserror::Error; use vtkio::Vtk; diff --git a/honeycomb-kernels/src/grisubal/model.rs b/honeycomb-kernels/src/grisubal/model.rs index 267642abd..9f3eacb04 100644 --- a/honeycomb-kernels/src/grisubal/model.rs +++ b/honeycomb-kernels/src/grisubal/model.rs @@ -11,15 +11,11 @@ use honeycomb_core::attributes::AttrSparseVec; use honeycomb_core::prelude::{ AttributeBind, AttributeUpdate, CoordsFloat, DartIdentifier, OrbitPolicy, Vertex2, }; -use std::collections::HashSet; use vtkio::{ model::{CellType, DataSet, VertexNumbers}, IOBuffer, Vtk, }; -#[cfg(doc)] -use honeycomb_core::prelude::CMap2; - // ------ CONTENT /// Structure used to index the overlapping grid's cases. @@ -207,193 +203,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. From 6bf637c6d2f38d5774a7ddbe38342f4ccb650b8d Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 11:48:35 +0100 Subject: [PATCH 07/12] cleaner timers --- honeycomb-kernels/src/grisubal/mod.rs | 71 ++++++------------------ honeycomb-kernels/src/grisubal/timers.rs | 45 +++++++++++++-- 2 files changed, 57 insertions(+), 59 deletions(-) diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index c0eb4a815..da16ca2af 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -55,15 +55,13 @@ pub(crate) mod clip; pub(crate) mod kernel; pub(crate) mod model; -#[cfg(feature = "profiling")] pub(crate) mod timers; -// ------ RE-EXPORTS - // ------ IMPORTS use crate::grisubal::clip::{clip_left, clip_right}; use crate::grisubal::model::{Boundary, Geometry2}; +use crate::grisubal::timers::{finish, start_timer, unsafe_time_section}; use honeycomb_core::cmap::{CMapBuilder, GridDescriptor}; use honeycomb_core::prelude::{CMap2, CoordsFloat}; use kernel::{ @@ -74,9 +72,6 @@ use kernel::{ use thiserror::Error; use vtkio::Vtk; -#[cfg(feature = "profiling")] -use timers::{unsafe_time_section, Section, TIMERS}; - // ------ CONTENT /// Post-processing clip operation. @@ -164,8 +159,7 @@ pub fn grisubal( grid_cell_sizes: [T; 2], clip: Clip, ) -> Result, GrisubalError> { - #[cfg(feature = "profiling")] - let mut instant = std::time::Instant::now(); + start_timer!(instant); // load geometry from file let geometry_vtk = match Vtk::import(file_path) { @@ -173,31 +167,26 @@ pub fn grisubal( Err(e) => panic!("E: could not open specified vtk file - {e}"), }; - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::ImportVTK); + unsafe_time_section!(instant, timers::Section::ImportVTK); // pre-processing let mut geometry = Geometry2::try_from(geometry_vtk)?; - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildGeometry); + unsafe_time_section!(instant, timers::Section::BuildGeometry); detect_orientation_issue(&geometry)?; - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::DetectOrientation); + unsafe_time_section!(instant, timers::Section::DetectOrientation); // compute an overlapping grid & remove redundant PoIs let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; let [cx, cy] = grid_cell_sizes; - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::ComputeOverlappingGrid); + unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid); remove_redundant_poi(&mut geometry, grid_cell_sizes, origin); - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::RemoveRedundantPoi); + unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi); // compute grid characteristics // build grid descriptor @@ -208,8 +197,7 @@ pub fn grisubal( .len_per_cell_y(cy) .origin(origin); - #[cfg(feature = "profiling")] - let mut kernel = std::time::Instant::now(); + start_timer!(kernel); // build initial map let mut cmap = CMapBuilder::default() @@ -218,8 +206,7 @@ pub fn grisubal( .build() .expect("E: unreachable"); // urneachable because grid dims are valid - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshInit); + unsafe_time_section!(instant, timers::Section::BuildMeshInit); // process the geometry @@ -231,8 +218,7 @@ pub fn grisubal( let (new_segments, intersection_metadata) = generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin); - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshIntersecData); + unsafe_time_section!(instant, timers::Section::BuildMeshIntersecData); // STEP 1.5 // precompute stuff to @@ -250,8 +236,7 @@ pub fn grisubal( insert_intersections(&mut cmap, &edge_intersec, &dart_slices); - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshInsertIntersec); + unsafe_time_section!(instant, timers::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 @@ -260,8 +245,7 @@ pub fn grisubal( let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts); - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::BuildMeshEdgeData); + unsafe_time_section!(instant, timers::Section::BuildMeshEdgeData); // STEP 4 // now that we have some segments that are directly defined between intersections, we can use some N-maps' @@ -271,11 +255,8 @@ pub fn grisubal( insert_edges_in_map(&mut cmap, &edges); - #[cfg(feature = "profiling")] - { - unsafe_time_section!(instant, Section::BuildMeshInsertEdge); - unsafe_time_section!(kernel, Section::BuildMeshTot); - } + unsafe_time_section!(instant, timers::Section::BuildMeshInsertEdge); + unsafe_time_section!(kernel, timers::Section::BuildMeshTot); // optional post-processing match clip { @@ -284,32 +265,12 @@ pub fn grisubal( Clip::None => {} } - #[cfg(feature = "profiling")] - unsafe_time_section!(instant, Section::Clip); + unsafe_time_section!(instant, timers::Section::Clip); // remove attribute used for clipping 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/timers.rs b/honeycomb-kernels/src/grisubal/timers.rs index c25786bad..f7068dc6b 100644 --- a/honeycomb-kernels/src/grisubal/timers.rs +++ b/honeycomb-kernels/src/grisubal/timers.rs @@ -2,12 +2,12 @@ //! //! **This code is only compiled if the `profiling` feature is enabled.** -/// Global timers for execution times per-section. #[cfg(feature = "profiling")] +/// Global timers for execution times per-section. pub(crate) static mut TIMERS: [Option; 13] = [None; 13]; -/// Kernel section. #[cfg(feature = "profiling")] +/// Kernel section. pub(crate) enum Section { ImportVTK = 0, BuildGeometry, @@ -24,14 +24,51 @@ pub(crate) enum Section { Cleanup, } -#[cfg(feature = "profiling")] +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[$sec as usize] = Some($inst.elapsed()); + 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; From f31907a7e3833e53bf82e1eb920783cad7eef4b3 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 13:55:13 +0100 Subject: [PATCH 08/12] move clip to kernel module --- .../src/grisubal/{ => kernel}/clip.rs | 1 + honeycomb-kernels/src/grisubal/kernel/mod.rs | 2 ++ honeycomb-kernels/src/grisubal/mod.rs | 23 +++++++++++-------- 3 files changed, 16 insertions(+), 10 deletions(-) rename honeycomb-kernels/src/grisubal/{ => kernel}/clip.rs (99%) diff --git a/honeycomb-kernels/src/grisubal/clip.rs b/honeycomb-kernels/src/grisubal/kernel/clip.rs similarity index 99% rename from honeycomb-kernels/src/grisubal/clip.rs rename to honeycomb-kernels/src/grisubal/kernel/clip.rs index 251751455..6b5c80ee9 100644 --- a/honeycomb-kernels/src/grisubal/clip.rs +++ b/honeycomb-kernels/src/grisubal/kernel/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/mod.rs b/honeycomb-kernels/src/grisubal/kernel/mod.rs index 8eedd5fe7..13bda5875 100644 --- a/honeycomb-kernels/src/grisubal/kernel/mod.rs +++ b/honeycomb-kernels/src/grisubal/kernel/mod.rs @@ -6,6 +6,7 @@ // ------ MODULES +mod clip; mod step0; mod step1; mod step2; @@ -15,6 +16,7 @@ mod step5; // ------ RE-EXPORTS +pub(crate) use clip::{clip_left, clip_right}; pub(crate) use step0::*; pub(crate) use step1::*; pub(crate) use step2::*; diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index da16ca2af..2f1c3f859 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -52,22 +52,25 @@ // ------ MODULE DECLARATIONS -pub(crate) mod clip; pub(crate) mod kernel; pub(crate) mod model; pub(crate) mod timers; // ------ IMPORTS -use crate::grisubal::clip::{clip_left, clip_right}; -use crate::grisubal::model::{Boundary, Geometry2}; -use crate::grisubal::timers::{finish, start_timer, unsafe_time_section}; -use honeycomb_core::cmap::{CMapBuilder, GridDescriptor}; -use honeycomb_core::prelude::{CMap2, CoordsFloat}; -use kernel::{ - 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, +use crate::grisubal::{ + kernel::{ + 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, + }, + model::{Boundary, Geometry2}, + timers::{finish, start_timer, unsafe_time_section}, +}; +use honeycomb_core::{ + cmap::{CMapBuilder, GridDescriptor}, + prelude::{CMap2, CoordsFloat}, }; use thiserror::Error; use vtkio::Vtk; From 060f14f98b0717c458ada6444950373ea6ac0a1f Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 14:22:08 +0100 Subject: [PATCH 09/12] cleanup main function code --- .../src/grisubal/kernel/step1.rs | 4 +- .../src/grisubal/kernel/step2.rs | 4 + .../src/grisubal/kernel/step3.rs | 1 + .../src/grisubal/kernel/step4.rs | 4 + .../src/grisubal/kernel/step5.rs | 2 + honeycomb-kernels/src/grisubal/mod.rs | 89 ++++++++----------- 6 files changed, 50 insertions(+), 54 deletions(-) diff --git a/honeycomb-kernels/src/grisubal/kernel/step1.rs b/honeycomb-kernels/src/grisubal/kernel/step1.rs index 6d6eacca1..1878505e7 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step1.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step1.rs @@ -1,6 +1,8 @@ //! Step 1 implementation //! -//! +//! 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 diff --git a/honeycomb-kernels/src/grisubal/kernel/step2.rs b/honeycomb-kernels/src/grisubal/kernel/step2.rs index 24f2ddb72..5002766c1 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step2.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step2.rs @@ -1,4 +1,8 @@ //! 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 diff --git a/honeycomb-kernels/src/grisubal/kernel/step3.rs b/honeycomb-kernels/src/grisubal/kernel/step3.rs index 737c6ea3c..d937241d5 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step3.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step3.rs @@ -1,5 +1,6 @@ //! Step 3 implementation //! +//! Insert the intersections into the map. // ------ IMPORTS diff --git a/honeycomb-kernels/src/grisubal/kernel/step4.rs b/honeycomb-kernels/src/grisubal/kernel/step4.rs index 6499ddb26..e632a390d 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step4.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step4.rs @@ -1,4 +1,8 @@ //! 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 diff --git a/honeycomb-kernels/src/grisubal/kernel/step5.rs b/honeycomb-kernels/src/grisubal/kernel/step5.rs index 9eaea1131..6ffd3446a 100644 --- a/honeycomb-kernels/src/grisubal/kernel/step5.rs +++ b/honeycomb-kernels/src/grisubal/kernel/step5.rs @@ -1,4 +1,6 @@ //! Step 5 implementation +// +//! Use the information computed at step 4 and insert all new edges into the map. // ------ IMPORTS diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 2f1c3f859..71edd0866 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -162,118 +162,101 @@ pub fn grisubal( grid_cell_sizes: [T; 2], clip: Clip, ) -> Result, GrisubalError> { + // 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); + //----/ - // pre-processing + // --- BUILD OUR MODEL FROM THE VTK IMPORT let mut geometry = Geometry2::try_from(geometry_vtk)?; - unsafe_time_section!(instant, timers::Section::BuildGeometry); + //----/ + // --- FIRST DETECTION OF ORIENTATION ISSUES detect_orientation_issue(&geometry)?; - unsafe_time_section!(instant, timers::Section::DetectOrientation); + //----/ - // compute an overlapping grid & remove redundant PoIs + // --- FIND AN OVERLAPPING GRID let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes)?; let [cx, cy] = grid_cell_sizes; - - unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid); - - remove_redundant_poi(&mut geometry, grid_cell_sizes, origin); - - unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi); - - // 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); + unsafe_time_section!(instant, timers::Section::ComputeOverlappingGrid); + //----/ + // --- REMOVE REDUNDANT PoIs + remove_redundant_poi(&mut geometry, grid_cell_sizes, origin); + unsafe_time_section!(instant, timers::Section::RemoveRedundantPoi); + //----/ + + // ------ START MAIN KERNEL TIMER start_timer!(kernel); - // build initial map + // --- BUILD THE GRID 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 - + .expect("E: unreachable"); // unreachable because grid dims are valid unsafe_time_section!(instant, timers::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. - + // --- STEP 1 & 2 + // (1) let (new_segments, intersection_metadata) = generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin); - - unsafe_time_section!(instant, timers::Section::BuildMeshIntersecData); - - // STEP 1.5 - // precompute stuff to - // - parallelize step 2 - // - make step 2 and step 3 independent from each other - + // (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); + //----/ - // 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. - + // --- STEP 3 insert_intersections(&mut cmap, &edge_intersec, &dart_slices); - unsafe_time_section!(instant, timers::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. - + // --- STEP 4 let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts); - unsafe_time_section!(instant, timers::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. - + // --- STEP 5 insert_edges_in_map(&mut cmap, &edges); - unsafe_time_section!(instant, timers::Section::BuildMeshInsertEdge); + //----/ + unsafe_time_section!(kernel, timers::Section::BuildMeshTot); + //-------/ - // optional post-processing + // --- CLIP match clip { Clip::Left => clip_left(&mut cmap)?, Clip::Right => clip_right(&mut cmap)?, Clip::None => {} } - unsafe_time_section!(instant, timers::Section::Clip); + //----/ - // remove attribute used for clipping + // CLEANUP cmap.remove_attribute_storage::(); - finish!(instant); + //-/ Ok(cmap) } From 84632bb270a0219eae16877d2a1663177cb6a474 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 31 Oct 2024 14:29:52 +0100 Subject: [PATCH 10/12] rename `kernel` module to `routines` --- honeycomb-kernels/src/grisubal/mod.rs | 6 +++--- honeycomb-kernels/src/grisubal/{kernel => routines}/clip.rs | 0 honeycomb-kernels/src/grisubal/{kernel => routines}/mod.rs | 0 .../src/grisubal/{kernel => routines}/step0.rs | 0 .../src/grisubal/{kernel => routines}/step1.rs | 0 .../src/grisubal/{kernel => routines}/step2.rs | 0 .../src/grisubal/{kernel => routines}/step3.rs | 0 .../src/grisubal/{kernel => routines}/step4.rs | 0 .../src/grisubal/{kernel => routines}/step5.rs | 0 honeycomb-kernels/src/grisubal/tests.rs | 4 ++-- 10 files changed, 5 insertions(+), 5 deletions(-) rename honeycomb-kernels/src/grisubal/{kernel => routines}/clip.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/mod.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step0.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step1.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step2.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step3.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step4.rs (100%) rename honeycomb-kernels/src/grisubal/{kernel => routines}/step5.rs (100%) diff --git a/honeycomb-kernels/src/grisubal/mod.rs b/honeycomb-kernels/src/grisubal/mod.rs index 71edd0866..ea4036188 100644 --- a/honeycomb-kernels/src/grisubal/mod.rs +++ b/honeycomb-kernels/src/grisubal/mod.rs @@ -52,20 +52,20 @@ // ------ MODULE DECLARATIONS -pub(crate) mod kernel; pub(crate) mod model; +pub(crate) mod routines; pub(crate) mod timers; // ------ IMPORTS use crate::grisubal::{ - kernel::{ + 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, }, - model::{Boundary, Geometry2}, timers::{finish, start_timer, unsafe_time_section}, }; use honeycomb_core::{ diff --git a/honeycomb-kernels/src/grisubal/kernel/clip.rs b/honeycomb-kernels/src/grisubal/routines/clip.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/clip.rs rename to honeycomb-kernels/src/grisubal/routines/clip.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/mod.rs b/honeycomb-kernels/src/grisubal/routines/mod.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/mod.rs rename to honeycomb-kernels/src/grisubal/routines/mod.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step0.rs b/honeycomb-kernels/src/grisubal/routines/step0.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step0.rs rename to honeycomb-kernels/src/grisubal/routines/step0.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step1.rs b/honeycomb-kernels/src/grisubal/routines/step1.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step1.rs rename to honeycomb-kernels/src/grisubal/routines/step1.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step2.rs b/honeycomb-kernels/src/grisubal/routines/step2.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step2.rs rename to honeycomb-kernels/src/grisubal/routines/step2.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step3.rs b/honeycomb-kernels/src/grisubal/routines/step3.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step3.rs rename to honeycomb-kernels/src/grisubal/routines/step3.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step4.rs b/honeycomb-kernels/src/grisubal/routines/step4.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step4.rs rename to honeycomb-kernels/src/grisubal/routines/step4.rs diff --git a/honeycomb-kernels/src/grisubal/kernel/step5.rs b/honeycomb-kernels/src/grisubal/routines/step5.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/kernel/step5.rs rename to honeycomb-kernels/src/grisubal/routines/step5.rs 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 From 1e73afb2bd47afd55f989083e0c7d6615b49fee8 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Fri, 8 Nov 2024 14:34:09 +0100 Subject: [PATCH 11/12] rename modules to describe steps --- .../{step1.rs => compute_intersecs.rs} | 0 .../{step4.rs => compute_new_edges.rs} | 0 .../{step3.rs => insert_intersecs.rs} | 0 .../{step5.rs => insert_new_edges.rs} | 0 .../src/grisubal/routines/mod.rs | 37 +++++++++++++------ .../routines/{step0.rs => pre_processing.rs} | 0 .../{step2.rs => process_intersecs_data.rs} | 0 7 files changed, 25 insertions(+), 12 deletions(-) rename honeycomb-kernels/src/grisubal/routines/{step1.rs => compute_intersecs.rs} (100%) rename honeycomb-kernels/src/grisubal/routines/{step4.rs => compute_new_edges.rs} (100%) rename honeycomb-kernels/src/grisubal/routines/{step3.rs => insert_intersecs.rs} (100%) rename honeycomb-kernels/src/grisubal/routines/{step5.rs => insert_new_edges.rs} (100%) rename honeycomb-kernels/src/grisubal/routines/{step0.rs => pre_processing.rs} (100%) rename honeycomb-kernels/src/grisubal/routines/{step2.rs => process_intersecs_data.rs} (100%) diff --git a/honeycomb-kernels/src/grisubal/routines/step1.rs b/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step1.rs rename to honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs diff --git a/honeycomb-kernels/src/grisubal/routines/step4.rs b/honeycomb-kernels/src/grisubal/routines/compute_new_edges.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step4.rs rename to honeycomb-kernels/src/grisubal/routines/compute_new_edges.rs diff --git a/honeycomb-kernels/src/grisubal/routines/step3.rs b/honeycomb-kernels/src/grisubal/routines/insert_intersecs.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step3.rs rename to honeycomb-kernels/src/grisubal/routines/insert_intersecs.rs diff --git a/honeycomb-kernels/src/grisubal/routines/step5.rs b/honeycomb-kernels/src/grisubal/routines/insert_new_edges.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step5.rs rename to honeycomb-kernels/src/grisubal/routines/insert_new_edges.rs diff --git a/honeycomb-kernels/src/grisubal/routines/mod.rs b/honeycomb-kernels/src/grisubal/routines/mod.rs index 13bda5875..5db8318d0 100644 --- a/honeycomb-kernels/src/grisubal/routines/mod.rs +++ b/honeycomb-kernels/src/grisubal/routines/mod.rs @@ -7,22 +7,35 @@ // ------ MODULES mod clip; -mod step0; -mod step1; -mod step2; -mod step3; -mod step4; -mod step5; +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}; -pub(crate) use step0::*; -pub(crate) use step1::*; -pub(crate) use step2::*; -pub(crate) use step3::*; -pub(crate) use step4::*; -pub(crate) use step5::*; // ------ IMPORTS diff --git a/honeycomb-kernels/src/grisubal/routines/step0.rs b/honeycomb-kernels/src/grisubal/routines/pre_processing.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step0.rs rename to honeycomb-kernels/src/grisubal/routines/pre_processing.rs diff --git a/honeycomb-kernels/src/grisubal/routines/step2.rs b/honeycomb-kernels/src/grisubal/routines/process_intersecs_data.rs similarity index 100% rename from honeycomb-kernels/src/grisubal/routines/step2.rs rename to honeycomb-kernels/src/grisubal/routines/process_intersecs_data.rs From 3f7228ce9bfa45469c3c8d23f3cc0ce2019ead03 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Fri, 8 Nov 2024 14:38:18 +0100 Subject: [PATCH 12/12] rename man_dist & diff methods --- honeycomb-kernels/src/grisubal/model.rs | 8 +++++--- .../src/grisubal/routines/compute_intersecs.rs | 4 ++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/honeycomb-kernels/src/grisubal/model.rs b/honeycomb-kernels/src/grisubal/model.rs index 9f3eacb04..2ed10c95b 100644 --- a/honeycomb-kernels/src/grisubal/model.rs +++ b/honeycomb-kernels/src/grisubal/model.rs @@ -26,14 +26,16 @@ use vtkio::{ pub struct GridCellId(pub usize, pub usize); impl GridCellId { - /// Compute the [Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry) between + /// Compute the [L1 / Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry) between /// two cells. - pub fn man_dist(lhs: &Self, rhs: &Self) -> usize { + 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 diff(lhs: &Self, rhs: &Self) -> (isize, isize) { + 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, diff --git a/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs b/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs index 1878505e7..07496be96 100644 --- a/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs +++ b/honeycomb-kernels/src/grisubal/routines/compute_intersecs.rs @@ -87,8 +87,8 @@ pub(crate) 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,