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/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/mod.rs b/honeycomb-kernels/src/dvr/mod.rs
new file mode 100644
index 000000000..124bc158e
--- /dev/null
+++ b/honeycomb-kernels/src/dvr/mod.rs
@@ -0,0 +1,88 @@
+//! Directional Vertex Relaxation implementation
+
+// ------ MODULE DECLARATIONS
+
+mod lambda;
+mod node;
+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<T: CoordsFloat>(
+    map: &mut CMap2<T>,
+    n_relaxations: usize,
+    dirs_relaxations: Option<Vec<Vector2<T>>>,
+) {
+    // 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<VertexIdentifier> = 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<T: CoordsFloat>(
+    map: &mut CMap2<T>,
+    n_relaxations: usize,
+    dirs_relaxations: Option<Vec<Vector2<T>>>,
+) {
+    // 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<VertexIdentifier> = 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/node.rs b/honeycomb-kernels/src/dvr/node.rs
new file mode 100644
index 000000000..032e5f4f8
--- /dev/null
+++ b/honeycomb-kernels/src/dvr/node.rs
@@ -0,0 +1,130 @@
+use enum_dispatch::enum_dispatch;
+use honeycomb_core::{
+    cmap::{CMap2, DartIdentifier, Orbit2, OrbitPolicy, VertexIdentifier, NULL_DART_ID},
+    prelude::{CoordsFloat, Vector2, Vertex2},
+};
+
+use super::quality;
+
+/// Node defined by a vertex and its incident 2-cells.
+pub struct Node {
+    pub id: VertexIdentifier,
+    cuts: Vec<Cut>,
+}
+
+impl Node {
+    pub fn new<T: CoordsFloat>(id: VertexIdentifier, map: &CMap2<T>) -> 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<T: CoordsFloat>(&self, map: &CMap2<T>) -> 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<T: CoordsFloat>(&self, map: &CMap2<T>, shift: &Vector2<T>, 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<T: CoordsFloat>(&self, map: &CMap2<T>, v: &Vertex2<T>) -> T;
+}
+
+#[enum_dispatch(CutQuality)]
+pub enum Cut {
+    CutTri,
+    CutQuad,
+    CutAny,
+}
+
+/// Triangular node cell.
+pub struct CutTri([VertexIdentifier; 2]);
+
+impl CutQuality for CutTri {
+    fn quality<T: CoordsFloat>(&self, map: &CMap2<T>, v: &Vertex2<T>) -> 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<T: CoordsFloat>(&self, map: &CMap2<T>, v: &Vertex2<T>) -> 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<VertexIdentifier>);
+
+impl CutQuality for CutAny {
+    fn quality<T: CoordsFloat>(&self, map: &CMap2<T>, v: &Vertex2<T>) -> 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
new file mode 100644
index 000000000..e19aa3b18
--- /dev/null
+++ b/honeycomb-kernels/src/dvr/quality.rs
@@ -0,0 +1,44 @@
+use honeycomb_core::prelude::{CoordsFloat, Vertex2};
+
+/// Compute the quality factor of a given triangle.
+pub fn quality_tri<T: CoordsFloat>(v1: &Vertex2<T>, v2: &Vertex2<T>, v3: &Vertex2<T>) -> 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_any<T: CoordsFloat>(vertices: &[Vertex2<T>]) -> 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
+}
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;