From b5d98f803121c3a9ce76165471601d5d1184b55f Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Tue, 29 Oct 2024 14:32:14 +0100 Subject: [PATCH 1/3] add dvr module & main functions skeletons --- honeycomb-kernels/src/dvr/maximizers.rs | 0 honeycomb-kernels/src/dvr/mod.rs | 87 +++++++++++++++++++++++++ honeycomb-kernels/src/dvr/quality.rs | 0 honeycomb-kernels/src/lib.rs | 1 + 4 files changed, 88 insertions(+) create mode 100644 honeycomb-kernels/src/dvr/maximizers.rs create mode 100644 honeycomb-kernels/src/dvr/mod.rs create mode 100644 honeycomb-kernels/src/dvr/quality.rs diff --git a/honeycomb-kernels/src/dvr/maximizers.rs b/honeycomb-kernels/src/dvr/maximizers.rs new file mode 100644 index 000000000..e69de29bb diff --git a/honeycomb-kernels/src/dvr/mod.rs b/honeycomb-kernels/src/dvr/mod.rs new file mode 100644 index 000000000..c28f142e0 --- /dev/null +++ b/honeycomb-kernels/src/dvr/mod.rs @@ -0,0 +1,87 @@ +//! Directional Vertex Relaxation implementation + +// ------ MODULE DECLARATIONS + +mod maximizers; +mod quality; + +// ------ IMPORTS + +use honeycomb_core::{ + cmap::{CMap2, DartIdentifier, Orbit2, OrbitPolicy, VertexIdentifier}, + prelude::{CoordsFloat, Vector2}, +}; + +// ------ CONTENT + +/// Regular DVR implementation. +#[allow(clippy::needless_for_each)] +pub fn dvr( + map: &mut CMap2, + n_relaxations: usize, + dirs_relaxations: Option>>, +) { + // DVR does not affect topology, so these IDs will stay valid during the entire run + // we can filter-out vertices on the boundary from the get-go + let vertices: Vec = map + .fetch_vertices() + .identifiers + .iter() + .filter(|vid| { + !Orbit2::new(map, OrbitPolicy::Vertex, **vid as DartIdentifier) + .any(|dart_id| map.is_i_free::<2>(dart_id)) + }) + .copied() + .collect(); + + // use arg relaxation dirs if passed, else default to cartesian dirs + let dirs = dirs_relaxations.unwrap_or(vec![Vector2::unit_x(), Vector2::unit_y()]); + + // using a for outside since this should not be parallelized + // iterator inside since it can be with changes to vertex selection + for k in 0..n_relaxations { + let _dir = dirs[k % dirs.len()]; + vertices.iter().for_each(|_vid| { + // compute the set of quality maximizers + // if card == 1 => set lambda_opt + // else find lambda_opt s.t. ... + // v <= v + dir * lambda_opt + todo!() + }); + } +} + +/// Approximate DVR implementation. +#[allow(clippy::needless_for_each)] +pub fn dvr_approximate( + map: &mut CMap2, + n_relaxations: usize, + dirs_relaxations: Option>>, +) { + // DVR does not affect topology, so these IDs will stay valid during the entire run + // we can filter-out vertices on the boundary from the get-go + let vertices: Vec = map + .fetch_vertices() + .identifiers + .iter() + .filter(|vid| { + !Orbit2::new(map, OrbitPolicy::Vertex, **vid as DartIdentifier) + .any(|dart_id| map.is_i_free::<2>(dart_id)) + }) + .copied() + .collect(); + + // use arg relaxation dirs if passed, else default to cartesian dirs + let dirs = dirs_relaxations.unwrap_or(vec![Vector2::unit_x(), Vector2::unit_y()]); + + // using a for outside since this should not be parallelized + // iterator inside since it can be with changes to vertex selection + for k in 0..n_relaxations { + let _dir = dirs[k % dirs.len()]; + vertices.iter().for_each(|_vid| { + // find lambda s.t. quality(lambda) > quality(0) + // v <= v + dir * lambda + todo!() + }); + } +} diff --git a/honeycomb-kernels/src/dvr/quality.rs b/honeycomb-kernels/src/dvr/quality.rs new file mode 100644 index 000000000..e69de29bb diff --git a/honeycomb-kernels/src/lib.rs b/honeycomb-kernels/src/lib.rs index ecb0f1b3b..ba4939dc0 100644 --- a/honeycomb-kernels/src/lib.rs +++ b/honeycomb-kernels/src/lib.rs @@ -21,6 +21,7 @@ // ------ MODULE DECLARATIONS +pub mod dvr; pub mod grisubal; pub mod splits; pub mod triangulation; From f51d5ee6f693b4905b5443950d9dc1f939bb7847 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 7 Nov 2024 08:03:12 +0100 Subject: [PATCH 2/3] WIP --- honeycomb-kernels/src/dvr/lambda.rs | 3 ++ honeycomb-kernels/src/dvr/maximizers.rs | 0 honeycomb-kernels/src/dvr/mod.rs | 3 +- honeycomb-kernels/src/dvr/node.rs | 11 +++++++ honeycomb-kernels/src/dvr/quality.rs | 44 +++++++++++++++++++++++++ 5 files changed, 60 insertions(+), 1 deletion(-) create mode 100644 honeycomb-kernels/src/dvr/lambda.rs delete mode 100644 honeycomb-kernels/src/dvr/maximizers.rs create mode 100644 honeycomb-kernels/src/dvr/node.rs diff --git a/honeycomb-kernels/src/dvr/lambda.rs b/honeycomb-kernels/src/dvr/lambda.rs new file mode 100644 index 000000000..ebda89465 --- /dev/null +++ b/honeycomb-kernels/src/dvr/lambda.rs @@ -0,0 +1,3 @@ +pub fn find_best_lambda() {} + +pub fn find_common_lambda() {} diff --git a/honeycomb-kernels/src/dvr/maximizers.rs b/honeycomb-kernels/src/dvr/maximizers.rs deleted file mode 100644 index e69de29bb..000000000 diff --git a/honeycomb-kernels/src/dvr/mod.rs b/honeycomb-kernels/src/dvr/mod.rs index c28f142e0..124bc158e 100644 --- a/honeycomb-kernels/src/dvr/mod.rs +++ b/honeycomb-kernels/src/dvr/mod.rs @@ -2,7 +2,8 @@ // ------ MODULE DECLARATIONS -mod maximizers; +mod lambda; +mod node; mod quality; // ------ IMPORTS diff --git a/honeycomb-kernels/src/dvr/node.rs b/honeycomb-kernels/src/dvr/node.rs new file mode 100644 index 000000000..767958ced --- /dev/null +++ b/honeycomb-kernels/src/dvr/node.rs @@ -0,0 +1,11 @@ +use honeycomb_core::cmap::VertexIdentifier; + +pub struct NodeTri { + v: VertexIdentifier, + neighbors: [VertexIdentifier; 2], +} + +pub struct Node { + v: VertexIdentifier, + neighbors: Vec, +} diff --git a/honeycomb-kernels/src/dvr/quality.rs b/honeycomb-kernels/src/dvr/quality.rs index e69de29bb..3e40a579d 100644 --- a/honeycomb-kernels/src/dvr/quality.rs +++ b/honeycomb-kernels/src/dvr/quality.rs @@ -0,0 +1,44 @@ +use honeycomb_core::prelude::{CoordsFloat, Vertex2}; + +/// Compute the quality factor of a given simple polygon. +pub fn quality(vertices: &[Vertex2]) -> T { + let last = &vertices[vertices.len() - 1]; + let blast = &vertices[vertices.len() - 2]; + let first = &vertices[0]; + // shoelace formula + let mut signed_area = vertices + .windows(3) + .map(|sl| { + let [v1, v2, v3] = sl else { unreachable!() }; + v2.x() * (v3.y() - v1.y()) + }) + .fold(T::zero(), |t1, t2| t1 + t2); + signed_area += last.x() * (first.y() - blast.y()); + signed_area += first.x() * (vertices[1].y() - last.y()); + signed_area /= T::from(2.0).unwrap(); + + let mut squared_edge_sum = vertices + .windows(2) + .map(|sl| { + let [v1, v2] = sl else { unreachable!() }; + (v2.x() - v1.x()).powi(2) + (v2.y() - v1.y()).powi(2) + }) + .fold(T::zero(), |t1, t2| t1 + t2); + squared_edge_sum += (first.x() - last.x()).powi(2) + (first.y() - last.y()).powi(2); + T::from(4.0 * 3.0_f64.sqrt()).unwrap() * signed_area / squared_edge_sum +} + +/// Compute the quality factor of a given triangle. +pub fn quality_tri(v1: &Vertex2, v2: &Vertex2, v3: &Vertex2) -> T { + let mut signed_area = + v1.x() * (v2.y() - v3.y()) + v2.x() * (v3.y() - v1.y()) + v3.x() * (v1.y() - v2.y()); + signed_area /= T::from(2.0).unwrap(); + + let squared_edge_sum = (v2.x() - v1.x()).powi(2) + + (v2.y() - v1.y()).powi(2) + + (v3.x() - v2.x()).powi(2) + + (v3.y() - v2.y()).powi(2) + + (v1.x() - v3.x()).powi(2) + + (v1.y() - v3.y()).powi(2); + T::from(4.0 * 3.0_f64.sqrt()).unwrap() * signed_area / squared_edge_sum +} From 1415534619b3712ffd1a6b5d3929f6efc6b3b01f Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Fri, 15 Nov 2024 12:04:25 +0100 Subject: [PATCH 3/3] implement node & cuts for quality computation --- Cargo.toml | 3 + honeycomb-kernels/Cargo.toml | 1 + honeycomb-kernels/src/dvr/node.rs | 133 +++++++++++++++++++++++++-- honeycomb-kernels/src/dvr/quality.rs | 32 +++---- 4 files changed, 146 insertions(+), 23 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 3cf830883..967b9ef09 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -45,6 +45,9 @@ dyn-clone = "1.0.17" num-traits = "0.2.19" vtkio = { version = "0.7.0-rc1", default-features = false } +# kernel +enum_dispatch = "0.3.13" + # benchmarks criterion = "0.5.1" iai-callgrind = "0.14.0" diff --git a/honeycomb-kernels/Cargo.toml b/honeycomb-kernels/Cargo.toml index ec6d71f3e..3b6c6402e 100644 --- a/honeycomb-kernels/Cargo.toml +++ b/honeycomb-kernels/Cargo.toml @@ -13,6 +13,7 @@ authors.workspace = true publish = true [dependencies] +enum_dispatch.workspace = true honeycomb-core = { workspace = true, features = ["io", "utils"] } num-traits.workspace = true thiserror.workspace = true diff --git a/honeycomb-kernels/src/dvr/node.rs b/honeycomb-kernels/src/dvr/node.rs index 767958ced..032e5f4f8 100644 --- a/honeycomb-kernels/src/dvr/node.rs +++ b/honeycomb-kernels/src/dvr/node.rs @@ -1,11 +1,130 @@ -use honeycomb_core::cmap::VertexIdentifier; +use enum_dispatch::enum_dispatch; +use honeycomb_core::{ + cmap::{CMap2, DartIdentifier, Orbit2, OrbitPolicy, VertexIdentifier, NULL_DART_ID}, + prelude::{CoordsFloat, Vector2, Vertex2}, +}; -pub struct NodeTri { - v: VertexIdentifier, - neighbors: [VertexIdentifier; 2], -} +use super::quality; +/// Node defined by a vertex and its incident 2-cells. pub struct Node { - v: VertexIdentifier, - neighbors: Vec, + pub id: VertexIdentifier, + cuts: Vec, +} + +impl Node { + pub fn new(id: VertexIdentifier, map: &CMap2) -> Self { + let cuts = Orbit2::new(map, OrbitPolicy::Vertex, id as DartIdentifier) + .map(|dart| { + let mut tmp = vec![]; + let mut crt = map.beta::<1>(dart); + while (crt != dart) || (crt != NULL_DART_ID) { + tmp.push(map.vertex_id(crt)); + crt = map.beta::<1>(dart) + } + match tmp.len() { + // TODO: raise a proper error + 0 | 1 => panic!(), + // 2 is a tri because we don't count the center vertex + 2 => CutTri(tmp.try_into().unwrap()).into(), + // ditto + 3 => CutQuad(tmp.try_into().unwrap()).into(), + _ => CutAny(tmp).into(), + } + }) + .collect(); + + Self { id, cuts } + } + + /// Compute an upper bound for displacement length, independently from the shift's direction. + pub fn lambda_max(&self, map: &CMap2) -> T { + let v_cent = map.vertex(self.id).unwrap(); + Orbit2::new(map, OrbitPolicy::Vertex, self.id as DartIdentifier) + .map(|dart| { + let v_neigh = map.vertex(map.vertex_id(map.beta::<2>(dart))).unwrap(); + (v_cent - v_neigh).norm() + }) + .max_by(|x, y| x.partial_cmp(y).unwrap()) + .unwrap() + } + + /// Compute the quality of the node. + /// + /// The original paper provides two definitions for quality: + /// - + /// - + /// + /// Control over definition usage is done through the boolean argument, `l2_mean`. + pub fn quality(&self, map: &CMap2, shift: &Vector2, l2_mean: bool) -> T { + let v = map.vertex(self.id).unwrap() + *shift; + if l2_mean { + self.cuts + .iter() + .map(|cut| cut.quality(map, &v)) + .map(|q| q.powi(2)) + .fold(T::zero(), |q1, q2| q1 + q2) + .sqrt() + } else { + self.cuts + .iter() + .map(|cut| cut.quality(map, &v)) + .min_by(|x, y| x.partial_cmp(y).unwrap()) + .unwrap() + } + } +} + +// crate shenanigans to go full functional programming while avoiding dyn objects collection + +#[enum_dispatch] +pub trait CutQuality { + fn quality(&self, map: &CMap2, v: &Vertex2) -> T; +} + +#[enum_dispatch(CutQuality)] +pub enum Cut { + CutTri, + CutQuad, + CutAny, +} + +/// Triangular node cell. +pub struct CutTri([VertexIdentifier; 2]); + +impl CutQuality for CutTri { + fn quality(&self, map: &CMap2, v: &Vertex2) -> T { + let [v1, v2] = [ + map.vertex(self.0[0]).unwrap(), + map.vertex(self.0[1]).unwrap(), + ]; + quality::quality_tri(v, &v1, &v2) + } +} + +/// Quadangular node cell. +pub struct CutQuad([VertexIdentifier; 3]); + +impl CutQuality for CutQuad { + fn quality(&self, map: &CMap2, v: &Vertex2) -> T { + // TODO: implement a proper hardcoded version for quads + let sl = [ + *v, + map.vertex(self.0[0]).unwrap(), + map.vertex(self.0[1]).unwrap(), + map.vertex(self.0[2]).unwrap(), + ]; + quality::quality_any(&sl) + } +} + +/// N-polygonal node cell, `n >= 5`. +pub struct CutAny(Vec); + +impl CutQuality for CutAny { + fn quality(&self, map: &CMap2, v: &Vertex2) -> T { + let mut vertices = vec![*v]; + vertices.extend(self.0.iter().map(|vid| map.vertex(*vid).unwrap())); + quality::quality_any(&vertices) + } } diff --git a/honeycomb-kernels/src/dvr/quality.rs b/honeycomb-kernels/src/dvr/quality.rs index 3e40a579d..e19aa3b18 100644 --- a/honeycomb-kernels/src/dvr/quality.rs +++ b/honeycomb-kernels/src/dvr/quality.rs @@ -1,7 +1,22 @@ use honeycomb_core::prelude::{CoordsFloat, Vertex2}; +/// Compute the quality factor of a given triangle. +pub fn quality_tri(v1: &Vertex2, v2: &Vertex2, v3: &Vertex2) -> T { + let mut signed_area = + v1.x() * (v2.y() - v3.y()) + v2.x() * (v3.y() - v1.y()) + v3.x() * (v1.y() - v2.y()); + signed_area /= T::from(2.0).unwrap(); + + let squared_edge_sum = (v2.x() - v1.x()).powi(2) + + (v2.y() - v1.y()).powi(2) + + (v3.x() - v2.x()).powi(2) + + (v3.y() - v2.y()).powi(2) + + (v1.x() - v3.x()).powi(2) + + (v1.y() - v3.y()).powi(2); + T::from(4.0 * 3.0_f64.sqrt()).unwrap() * signed_area / squared_edge_sum +} + /// Compute the quality factor of a given simple polygon. -pub fn quality(vertices: &[Vertex2]) -> T { +pub fn quality_any(vertices: &[Vertex2]) -> T { let last = &vertices[vertices.len() - 1]; let blast = &vertices[vertices.len() - 2]; let first = &vertices[0]; @@ -27,18 +42,3 @@ pub fn quality(vertices: &[Vertex2]) -> T { squared_edge_sum += (first.x() - last.x()).powi(2) + (first.y() - last.y()).powi(2); T::from(4.0 * 3.0_f64.sqrt()).unwrap() * signed_area / squared_edge_sum } - -/// Compute the quality factor of a given triangle. -pub fn quality_tri(v1: &Vertex2, v2: &Vertex2, v3: &Vertex2) -> T { - let mut signed_area = - v1.x() * (v2.y() - v3.y()) + v2.x() * (v3.y() - v1.y()) + v3.x() * (v1.y() - v2.y()); - signed_area /= T::from(2.0).unwrap(); - - let squared_edge_sum = (v2.x() - v1.x()).powi(2) - + (v2.y() - v1.y()).powi(2) - + (v3.x() - v2.x()).powi(2) - + (v3.y() - v2.y()).powi(2) - + (v1.x() - v3.x()).powi(2) - + (v1.y() - v3.y()).powi(2); - T::from(4.0 * 3.0_f64.sqrt()).unwrap() * signed_area / squared_edge_sum -}