diff --git a/Cargo.lock b/Cargo.lock
index 04975f27..28a1273b 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -161,6 +161,7 @@ dependencies = [
  "num",
  "rand",
  "rand_pcg",
+ "rayon",
 ]
 
 [[package]]
diff --git a/tools/Cargo.toml b/tools/Cargo.toml
index 4c1133d0..77a41c52 100644
--- a/tools/Cargo.toml
+++ b/tools/Cargo.toml
@@ -30,3 +30,4 @@ anyhow = { version = "1", default-features = false, features = ["std"] }
 # Other utilities
 itertools = { version = "0.10", default-features = false }
 num = { version = "0.4", default-features = false }
+rayon = "1.3"
diff --git a/tools/src/bin/mesh-part/main.rs b/tools/src/bin/mesh-part/main.rs
index 9a44c088..faa4e539 100644
--- a/tools/src/bin/mesh-part/main.rs
+++ b/tools/src/bin/mesh-part/main.rs
@@ -1,75 +1,29 @@
 use anyhow::Context as _;
 use anyhow::Result;
-use coupe::RunInfo;
+use coupe::Partition as _;
+use coupe::PointND;
+use mesh_io::medit::Mesh;
 use mesh_io::weight;
+use mesh_io::weight::Array;
+use rayon::iter::IntoParallelRefIterator as _;
+use rayon::iter::ParallelIterator as _;
 use std::env;
 use std::fs;
 use std::io;
+use std::iter::Sum;
 
-macro_rules! coupe_part {
-    ( $algo:expr, $dimension:expr, $points:expr, $weights:expr ) => {{
-        match $dimension {
-            2 => {
-                coupe::Partitioner::<2>::partition(&$algo, $points.as_slice(), $weights).into_ids()
-            }
-            3 => {
-                coupe::Partitioner::<3>::partition(&$algo, $points.as_slice(), $weights).into_ids()
-            }
-            _ => anyhow::bail!("only 2D and 3D meshes are supported"),
-        }
-    }};
-    ( $algo:expr, $problem:expr ) => {{
-        let weights: Vec<f64> = match &$problem.weights {
-            weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(),
-            weight::Array::Integers(_) => anyhow::bail!("integers are not supported by coupe yet"),
-        };
-        coupe_part!($algo, $problem.dimension, $problem.points, &weights)
-    }};
+struct Problem<const D: usize, W> {
+    points: Vec<PointND<D>>,
+    weights: Vec<Vec<W>>,
 }
 
-macro_rules! coupe_part_improve {
-    ( $algo:expr, $partition:expr, $dimension:expr, $points:expr, $weights:expr ) => {{
-        let ids = match $dimension {
-            2 => {
-                let points = coupe::PointsView::to_points_nd($points.as_slice());
-                let partition =
-                    coupe::partition::Partition::from_ids(points, $weights, $partition.to_vec());
-                coupe::PartitionImprover::<2>::improve_partition(&$algo, partition).into_ids()
-            }
-            3 => {
-                let points = coupe::PointsView::to_points_nd($points.as_slice());
-                let partition =
-                    coupe::partition::Partition::from_ids(points, $weights, $partition.to_vec());
-                coupe::PartitionImprover::<3>::improve_partition(&$algo, partition).into_ids()
-            }
-            _ => anyhow::bail!("only 2D and 3D meshes are supported"),
-        };
-        $partition.copy_from_slice(&ids);
-    }};
-    ( $algo:expr, $partition:expr, $problem:expr ) => {{
-        let weights: Vec<f64> = match &$problem.weights {
-            weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(),
-            weight::Array::Integers(_) => anyhow::bail!("integers are not supported by coupe yet"),
-        };
-        coupe_part_improve!(
-            $algo,
-            $partition,
-            $problem.dimension,
-            $problem.points,
-            &weights
-        )
-    }};
-}
+type Algorithm<const D: usize, W> = Box<dyn FnMut(&mut [usize], &Problem<D, W>) -> Result<()>>;
 
-struct Problem {
-    dimension: usize,
-    points: Vec<f64>,
-    weights: weight::Array,
-}
-
-type Algorithm = Box<dyn FnMut(&mut [usize], &Problem) -> Result<RunInfo>>;
-
-fn parse_algorithm(spec: String) -> Result<Algorithm> {
+fn parse_algorithm<const D: usize, W>(spec: String) -> Result<Algorithm<D, W>>
+where
+    W: Copy + Send + Sync + std::fmt::Debug + Default,
+    W: Ord + Sum + num::traits::NumAssign + num::ToPrimitive,
+{
     let mut args = spec.split(',');
     let name = args.next().context("empty algorithm spec")?;
 
@@ -100,115 +54,94 @@ fn parse_algorithm(spec: String) -> Result<Algorithm> {
                 bytes.resize(32_usize, 0_u8);
                 bytes.try_into().unwrap()
             };
-            let mut rng = rand_pcg::Pcg64::from_seed(seed);
+            let rng = rand_pcg::Pcg64::from_seed(seed);
+            let algo = coupe::Random { rng, part_count };
             Box::new(move |partition, problem| {
-                let algo = coupe::Random::new(&mut rng, part_count);
-                let weights = vec![0.0; problem.points.len()];
-                let ids = coupe_part!(algo, problem.dimension, problem.points, &weights);
-                partition.copy_from_slice(&ids);
-                Ok(RunInfo::default())
+                algo.partition(partition, ())?;
+                Ok(())
             })
         }
         "greedy" => {
-            let part_count = required(usize_arg(args.next()))?;
+            let algo = coupe::Greedy {
+                part_count: required(usize_arg(args.next()))?,
+            };
             Box::new(move |partition, problem| {
-                let ids = coupe_part!(coupe::Greedy::new(part_count), problem);
-                partition.copy_from_slice(&ids);
-                Ok(RunInfo::default())
+                algo.partition(partition, problem.weights.iter().map(|w| w[0]))?;
+                Ok(())
             })
         }
         "kk" => {
-            let part_count = required(usize_arg(args.next()))?;
+            let algo = coupe::KarmarkarKarp {
+                part_count: required(usize_arg(args.next()))?,
+            };
             Box::new(move |partition, problem| {
-                let ids = coupe_part!(coupe::KarmarkarKarp::new(part_count), problem);
-                partition.copy_from_slice(&ids);
-                Ok(RunInfo::default())
+                algo.partition(partition, problem.weights.iter().map(|w| w[0]))?;
+                Ok(())
             })
         }
         "vn-best" => {
-            let part_count = required(usize_arg(args.next()))?;
+            let algo = coupe::VnBest {
+                part_count: required(usize_arg(args.next()))?,
+            };
             Box::new(move |partition, problem| {
-                coupe_part_improve!(coupe::VnBest::new(part_count), partition, problem);
-                Ok(RunInfo::default())
+                algo.partition(partition, problem.weights.iter().map(|w| w[0]))?;
+                Ok(())
             })
         }
         "rcb" => {
-            let iter_count = required(usize_arg(args.next()))?;
+            let algo = coupe::Rcb {
+                iter_count: required(usize_arg(args.next()))?,
+            };
             Box::new(move |partition, problem| {
-                let ids = coupe_part!(coupe::Rcb::new(iter_count), problem);
-                partition.copy_from_slice(&ids);
-                Ok(RunInfo::default())
+                algo.partition(
+                    partition,
+                    (
+                        problem.points.par_iter().cloned(),
+                        problem.weights.par_iter().map(|w| w[0]),
+                    ),
+                )?;
+                Ok(())
             })
         }
         "hilbert" => {
-            let num_partitions = required(usize_arg(args.next()))?;
-            let order = optional(usize_arg(args.next()), 12)? as u32;
-            Box::new(move |partition, problem| {
-                let algo = coupe::HilbertCurve {
-                    num_partitions,
-                    order,
-                };
-                let weights: Vec<f64> = match &problem.weights {
-                    weight::Array::Floats(ws) => ws.iter().map(|weight| weight[0]).collect(),
-                    weight::Array::Integers(_) => {
-                        anyhow::bail!("hilbert is not implemented for integers")
-                    }
-                };
-                let res = match problem.dimension {
-                    2 => coupe::Partitioner::<2>::partition(
-                        &algo,
-                        problem.points.as_slice(),
-                        &weights,
-                    )
-                    .into_ids(),
-                    _ => anyhow::bail!("hilbert is only implemented in 2D"),
-                };
-                partition.copy_from_slice(&res);
-                Ok(RunInfo::default())
-            })
+            let algo = coupe::HilbertCurve {
+                part_count: required(usize_arg(args.next()))?,
+                order: optional(usize_arg(args.next()), 12)? as u32,
+            };
+            if D == 2 {
+                Box::new(move |partition, problem| {
+                    let weights: Vec<_> = problem
+                        .weights
+                        .iter()
+                        .map(|w| w[0].to_f64().unwrap())
+                        .collect();
+                    //algo.partition(partition, (&problem.points, &weights))?;
+                    Ok(())
+                })
+            } else {
+                Box::new(move |partition, problem| anyhow::bail!("hilbert is restricted to 2D"))
+            }
         }
         _ => anyhow::bail!("unknown algorithm {:?}", name),
     })
 }
 
-fn main() -> Result<()> {
-    let mut options = getopts::Options::new();
-    options.optflag("h", "help", "print this help menu");
-    options.optmulti(
-        "a",
-        "algorithm",
-        "name of the algorithm to run, see ALGORITHMS",
-        "NAME",
-    );
-    options.optopt("m", "mesh", "mesh file", "FILE");
-    options.optopt("w", "weights", "weight file", "FILE");
-
-    let matches = options.parse(env::args().skip(1))?;
-
-    if matches.opt_present("h") {
-        eprintln!("{}", options.usage("Usage: mesh-part [options]"));
-        eprint!(include_str!("help_after.txt"));
-        return Ok(());
-    }
-
-    let algorithms: Vec<_> = matches
-        .opt_strs("a")
-        .into_iter()
-        .map(parse_algorithm)
-        .collect::<Result<_>>()?;
-
-    let mesh_file = matches
-        .opt_str("m")
-        .context("missing required option 'mesh'")?;
-    let mesh = mesh_io::medit::Mesh::from_file(mesh_file).context("failed to read mesh file")?;
-
+fn main_d_w<const D: usize, W>(
+    matches: getopts::Matches,
+    mesh: Mesh,
+    weights: Vec<Vec<W>>,
+) -> Result<()>
+where
+    W: Copy + Send + Sync + std::fmt::Debug + Default,
+    W: Ord + Sum + num::traits::NumAssign + num::ToPrimitive,
+{
     let points: Vec<_> = mesh
         .elements()
         .filter_map(|(element_type, nodes, _element_ref)| {
             if element_type.dimension() != mesh.dimension() {
                 return None;
             }
-            let mut barycentre = vec![0.0; mesh.dimension()];
+            let mut barycentre = [0.0; D];
             for node_idx in nodes {
                 let node_coordinates = mesh.node(*node_idx);
                 for (bc_coord, node_coord) in barycentre.iter_mut().zip(node_coordinates) {
@@ -218,11 +151,27 @@ fn main() -> Result<()> {
             for bc_coord in &mut barycentre {
                 *bc_coord = *bc_coord / nodes.len() as f64;
             }
-            Some(barycentre)
+            Some(PointND::from(barycentre))
         })
-        .flatten()
         .collect();
 
+    let mut partition = vec![0; points.len()];
+    let problem = Problem { points, weights };
+
+    for algorithm in matches.opt_strs("a") {
+        let mut algorithm = parse_algorithm(algorithm)?;
+        algorithm(&mut partition, &problem).context("failed to apply algorithm")?;
+    }
+
+    let stdout = io::stdout();
+    let stdout = stdout.lock();
+    let stdout = io::BufWriter::new(stdout);
+    mesh_io::partition::write(stdout, partition).context("failed to print partition")?;
+
+    Ok(())
+}
+
+fn main_d<const D: usize>(matches: getopts::Matches, mesh: Mesh) -> Result<()> {
     let weight_file = matches
         .opt_str("w")
         .context("missing required option 'weights'")?;
@@ -230,21 +179,42 @@ fn main() -> Result<()> {
     let weight_file = io::BufReader::new(weight_file);
     let weights = weight::read(weight_file).context("failed to read weight file")?;
 
-    let problem = Problem {
-        dimension: mesh.dimension(),
-        points,
-        weights,
-    };
-    let mut partition = vec![0; problem.points.len() / problem.dimension];
+    match weights {
+        Array::Integers(is) => main_d_w(matches, mesh, is),
+        Array::Floats(fs) => main_d_w(matches, mesh, fs),
+    }
+}
 
-    for mut algorithm in algorithms {
-        algorithm(&mut partition, &problem).context("failed to apply algorithm")?;
+fn main() -> Result<()> {
+    let mut options = getopts::Options::new();
+    options.optflag("h", "help", "print this help menu");
+    options.optmulti(
+        "a",
+        "algorithm",
+        "name of the algorithm to run, see ALGORITHMS",
+        "NAME",
+    );
+    options.optopt("m", "mesh", "mesh file", "FILE");
+    options.optopt("w", "weights", "weight file", "FILE");
+
+    let matches = options.parse(env::args().skip(1))?;
+
+    if matches.opt_present("h") {
+        eprintln!("{}", options.usage("Usage: mesh-part [options]"));
+        eprint!(include_str!("help_after.txt"));
+        return Ok(());
     }
 
-    let stdout = io::stdout();
-    let stdout = stdout.lock();
-    let stdout = io::BufWriter::new(stdout);
-    mesh_io::partition::write(stdout, partition).context("failed to print partition")?;
+    let mesh_file = matches
+        .opt_str("m")
+        .context("missing required option 'mesh'")?;
+    let mesh = Mesh::from_file(mesh_file).context("failed to read mesh file")?;
+
+    match mesh.dimension() {
+        2 => main_d::<2>(matches, mesh)?,
+        3 => main_d::<3>(matches, mesh)?,
+        n => anyhow::bail!("expected 2D or 3D mesh, got a {n}D mesh"),
+    }
 
     Ok(())
 }