diff --git a/applications/shift_vertices/cli.rs b/applications/shift_vertices/cli.rs
index 947933e6a..bc22fca69 100644
--- a/applications/shift_vertices/cli.rs
+++ b/applications/shift_vertices/cli.rs
@@ -1,6 +1,6 @@
 use std::{num::NonZero, path::PathBuf};
 
-use clap::Parser;
+use clap::{Parser, Subcommand};
 
 use applications::FileFormat;
 
@@ -22,4 +22,25 @@ pub struct Cli {
     /// Execute benchmarks using `f32` instead of the default `f64`
     #[arg(long("simple-precision"))]
     pub simple_precision: bool,
+    #[command(subcommand)]
+    pub command: Option<SmoothingType>,
+}
+
+#[derive(Subcommand)]
+pub enum SmoothingType {
+    /// Laplace smoothing
+    Laplace {
+        /// Scale coefficient
+        #[arg(long, default_value_t = 0.5)]
+        lambda: f64,
+    },
+    /// Taubin smoothing
+    Taubin {
+        /// Scale coefficient
+        #[arg(long, default_value_t = 0.6307)]
+        lambda: f64,
+        /// Pass-band parameter
+        #[arg(short, long("pass-band"), default_value_t = 0.1)]
+        k: f64,
+    },
 }
diff --git a/applications/shift_vertices/internals.rs b/applications/shift_vertices/internals.rs
index 00104a757..064b1d6e8 100644
--- a/applications/shift_vertices/internals.rs
+++ b/applications/shift_vertices/internals.rs
@@ -1,7 +1,7 @@
 use honeycomb::{
     prelude::{
         CMap2, CoordsFloat, DartIdType, NULL_DART_ID, OrbitPolicy, VertexIdType,
-        remeshing::move_vertex_to_average,
+        remeshing::{move_vertex_to_average, neighbor_based_smooth},
     },
     stm::atomically,
 };
@@ -68,3 +68,70 @@ pub fn shift<T: CoordsFloat>(
         }
     }
 }
+
+pub fn laplace<T: CoordsFloat>(
+    map: &CMap2<T>,
+    graph: &[(VertexIdType, Vec<VertexIdType>)],
+    n_rounds: usize,
+    lambda: T,
+) {
+    println!(" Round | process_time | throughput(vertex/s)");
+    // main loop
+    let mut round = 0;
+    let mut process_time;
+    let n_v = graph.len();
+    loop {
+        let instant = std::time::Instant::now();
+        graph.par_iter().for_each(|(vid, neigh)| {
+            atomically(|t| neighbor_based_smooth(t, map, *vid, neigh, lambda));
+        });
+        process_time = instant.elapsed().as_secs_f64();
+        println!(
+            " {:>5} | {:>12.6e} | {:>20.6e}",
+            round,
+            process_time,
+            n_v as f64 / process_time,
+        );
+
+        round += 1;
+        if round >= n_rounds {
+            break;
+        }
+    }
+}
+
+pub fn taubin<T: CoordsFloat>(
+    map: &CMap2<T>,
+    graph: &[(VertexIdType, Vec<VertexIdType>)],
+    n_rounds: usize,
+    lambda: T,
+    k: T,
+) {
+    println!(" Round | process_time | throughput(vertex/s)");
+    // main loop
+    let mut round = 0;
+    let mut process_time;
+    let n_v = graph.len();
+
+    let mu = T::one() / (k - T::one() / lambda);
+
+    loop {
+        let instant = std::time::Instant::now();
+        let scale = if round % 2 == 0 { lambda } else { mu };
+        graph.par_iter().for_each(|(vid, neigh)| {
+            atomically(|t| neighbor_based_smooth(t, map, *vid, neigh, scale));
+        });
+        process_time = instant.elapsed().as_secs_f64();
+        println!(
+            " {:>5} | {:>12.6e} | {:>20.6e}",
+            round,
+            process_time,
+            n_v as f64 / process_time,
+        );
+
+        round += 1;
+        if round >= n_rounds {
+            break;
+        }
+    }
+}
diff --git a/applications/shift_vertices/main.rs b/applications/shift_vertices/main.rs
index 9db905509..3ed733cff 100644
--- a/applications/shift_vertices/main.rs
+++ b/applications/shift_vertices/main.rs
@@ -17,9 +17,21 @@ fn main() {
     let cli = cli::Cli::parse();
 
     if cli.simple_precision {
-        run_bench::<f32>(cli.input, cli.n_rounds.get(), cli.sort, cli.save_as);
+        run_bench::<f32>(
+            cli.input,
+            cli.n_rounds.get(),
+            cli.sort,
+            cli.save_as,
+            cli.command,
+        );
     } else {
-        run_bench::<f64>(cli.input, cli.n_rounds.get(), cli.sort, cli.save_as);
+        run_bench::<f64>(
+            cli.input,
+            cli.n_rounds.get(),
+            cli.sort,
+            cli.save_as,
+            cli.command,
+        );
     }
 }
 
@@ -28,6 +40,7 @@ fn run_bench<T: CoordsFloat>(
     n_rounds: usize,
     sort: bool,
     save: Option<FileFormat>,
+    smoothing: Option<cli::SmoothingType>,
 ) {
     let (map, input_hash, init_time) = init_2d_map_from_file::<T>(input.clone());
 
@@ -41,12 +54,51 @@ fn run_bench<T: CoordsFloat>(
         get_num_threads().unwrap_or(1)
     );
     println!("|-> # of rounds: {}", n_rounds);
+
+    match smoothing {
+        None => {
+            println!("|-> smoothing  : neighbor average",);
+        }
+        Some(cli::SmoothingType::Laplace { lambda }) => {
+            println!("|-> smoothing  : Laplace");
+            println!("|-> scale (λ)  : {}", lambda);
+            assert!(
+                lambda > 0.0 && lambda < 1.0,
+                "lambda must verify `0 < lambda < 1`"
+            );
+        }
+        Some(cli::SmoothingType::Taubin { lambda, k }) => {
+            println!("|-> smoothing  : Taubin");
+            println!("|-> scale (λ)  : {lambda}");
+            println!("|-> pass (k)   : {k}");
+            assert!(
+                lambda > 0.0 && lambda < 1.0,
+                "lambda must verify `0 < lambda < 1`"
+            );
+        }
+    }
     println!("|-+ init time  :");
     println!("| |->   map built in {}ms", init_time.as_millis());
 
     prof_start!("HCBENCH_SHIFT");
     let graph = internals::build_vertex_graph(&map, sort);
-    internals::shift(&map, &graph, n_rounds);
+    match smoothing {
+        None => {
+            internals::shift(&map, &graph, n_rounds);
+        }
+        Some(cli::SmoothingType::Laplace { lambda }) => {
+            internals::laplace(&map, &graph, n_rounds, T::from(lambda).unwrap());
+        }
+        Some(cli::SmoothingType::Taubin { lambda, k }) => {
+            internals::taubin(
+                &map,
+                &graph,
+                n_rounds,
+                T::from(lambda).unwrap(),
+                T::from(k).unwrap(),
+            );
+        }
+    }
     prof_stop!("HCBENCH_SHIFT");
 
     finalize_2d(map, save);
diff --git a/honeycomb-kernels/src/remeshing/mod.rs b/honeycomb-kernels/src/remeshing/mod.rs
index 22721777e..b8b81a791 100644
--- a/honeycomb-kernels/src/remeshing/mod.rs
+++ b/honeycomb-kernels/src/remeshing/mod.rs
@@ -17,7 +17,7 @@ mod swap;
 pub use capture::{ClassificationError, capture_geometry, classify_capture};
 pub use collapse::{EdgeCollapseError, collapse_edge};
 pub use cut::{cut_inner_edge, cut_outer_edge};
-pub use relaxation::move_vertex_to_average;
+pub use relaxation::{move_vertex_to_average, neighbor_based_smooth};
 pub use swap::{EdgeSwapError, swap_edge};
 
 #[cfg(test)]
diff --git a/honeycomb-kernels/src/remeshing/relaxation.rs b/honeycomb-kernels/src/remeshing/relaxation.rs
index 546de9c78..46f3dd1b0 100644
--- a/honeycomb-kernels/src/remeshing/relaxation.rs
+++ b/honeycomb-kernels/src/remeshing/relaxation.rs
@@ -1,6 +1,6 @@
 use honeycomb_core::{
     cmap::{CMap2, VertexIdType},
-    geometry::{CoordsFloat, Vertex2},
+    geometry::{CoordsFloat, Vector2},
     stm::{StmClosureResult, Transaction},
 };
 
@@ -36,17 +36,64 @@ pub fn move_vertex_to_average<T: CoordsFloat>(
     vid: VertexIdType,
     others: &[VertexIdType],
 ) -> StmClosureResult<()> {
-    if others.is_empty() {
-        return Ok(());
-    }
-    let mut new_val = Vertex2::default();
-    for v in others {
-        let vertex = map.read_vertex(t, *v)?.unwrap();
-        new_val.0 += vertex.0;
-        new_val.1 += vertex.1;
+    neighbor_based_smooth(t, map, vid, others, T::one())
+}
+
+/// Generic neighbor-based vertex smoothing function.
+///
+/// This function smooths the vertex position by moving it toward the average of its neighbors'
+/// positions weighted by lambda.
+///
+/// Note that it is up to the user to provide a correct list of neighbor IDs, and "acceptable"
+/// lambda parameter. For example:
+///
+/// - `lambda == 1` nullifies the influence of the original vertex position,
+/// - `0 < lambda < 1` results in a Laplacian smoothing.
+///
+/// # Arguments
+///
+/// - `t: &mut Transaction` -- Associated transaction.
+/// - `map: &mut CMap2` -- Edited map.
+/// - `vid: VertexIdType` -- Vertex to move.
+/// - `neighbors_id: &[VertexIdType]` -- List of vertex to compute the average from.
+/// - `lambda: T` -- Coefficient weighting the applied offset.
+///
+/// # Errors
+///
+/// This function will abort and raise an error if the transaction cannot be completed.
+///
+/// # Panics
+///
+/// This function may panic if one vertex in the `neighbors_id` list has no associated coordinates.
+#[inline]
+pub fn neighbor_based_smooth<T: CoordsFloat>(
+    t: &mut Transaction,
+    map: &CMap2<T>,
+    vid: VertexIdType,
+    neighbors_id: &[VertexIdType],
+    lambda: T,
+) -> StmClosureResult<()> {
+    let p = map
+        .read_vertex(t, vid)?
+        .expect("E: no coordinates associated to vertex ID");
+
+    let n = neighbors_id.len();
+    let mut neighbors: smallvec::SmallVec<_, 16> = smallvec::SmallVec::with_capacity(n);
+    for &nid in neighbors_id {
+        neighbors.push(
+            map.read_vertex(t, nid)?
+                .expect("E: no coordinates associated to vertex ID"),
+        );
     }
-    new_val.0 /= T::from(others.len()).unwrap();
-    new_val.1 /= T::from(others.len()).unwrap();
-    map.write_vertex(t, vid, new_val)?;
+
+    let delta = neighbors
+        .into_iter()
+        .map(|v| v - p)
+        .fold(Vector2::default(), |a, b| a + b)
+        * lambda
+        / T::from(n).unwrap();
+
+    map.write_vertex(t, vid, p + delta)?;
+
     Ok(())
 }