diff --git a/cpp/dolfinx/mesh/MeshTags.h b/cpp/dolfinx/mesh/MeshTags.h
index 32cb2de6a3..3026c25815 100644
--- a/cpp/dolfinx/mesh/MeshTags.h
+++ b/cpp/dolfinx/mesh/MeshTags.h
@@ -127,8 +127,8 @@ class MeshTags
 };
 
 /// @brief Create MeshTags from arrays
-/// @param[in] topology Mesh topology that the tags are associated with
-/// @param[in] dim Topological dimension of tagged entities
+/// @param[in] topology Mesh topology that the tags are associated with.
+/// @param[in] dim Topological dimension of tagged entities.
 /// @param[in] entities Local vertex indices for tagged entities.
 /// @param[in] values Tag values for each entity in `entities`. The
 /// length of `values` must be equal to number of rows in `entities`.
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 32f18d8982..430bcc495b 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2010-2023 Garth N. Wells
+// Copyright (C) 2010-2024 Garth N. Wells and Paul T. Kühner
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -6,6 +6,7 @@
 
 #pragma once
 
+#include "dolfinx/graph/AdjacencyList.h"
 #include "dolfinx/mesh/Mesh.h"
 #include "dolfinx/mesh/Topology.h"
 #include "dolfinx/mesh/cell_types.h"
@@ -15,85 +16,72 @@
 #include <algorithm>
 #include <concepts>
 #include <optional>
+#include <spdlog/spdlog.h>
 #include <utility>
 
 namespace dolfinx::refinement
 {
-namespace impl
-{
-/// @brief create the refined mesh by optionally redistributing it
-template <std::floating_point T>
-mesh::Mesh<T>
-create_refined_mesh(const mesh::Mesh<T>& mesh,
-                    const graph::AdjacencyList<std::int64_t>& cell_adj,
-                    std::span<const T> new_vertex_coords,
-                    std::array<std::size_t, 2> xshape, bool redistribute,
-                    mesh::GhostMode ghost_mode)
-{
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    // No parallel construction necessary, i.e. redistribute also has no
-    // effect
-    return mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                             mesh.geometry().cmap(), new_vertex_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    // Let partition handle the parallel construction of the mesh
-    return partition<T>(mesh, cell_adj, new_vertex_coords, xshape, redistribute,
-                        ghost_mode);
-  }
-}
-} // namespace impl
-
-/// @brief Refine with markers, optionally redistributing, and
-/// optionally calculating the parent-child relationships.
+/// @brief Refine a mesh with markers.
+///
+/// The refined mesh can be optionally re-partitioned across processes.
+/// Passing `nullptr` for `partitioner`, refined cells will be on the
+/// same process as the parent cell.
 ///
-/// @param[in] mesh Input mesh to be refined
-/// @param[in] edges Optional indices of the edges that should be split
-/// by this refinement, if optional is not set, a uniform refinement is
-/// performed, same behavior as passing a list of all indices.
-/// @param[in] redistribute Flag to call the Mesh Partitioner to
-/// redistribute after refinement.
-/// @param[in] ghost_mode Ghost mode of the refined mesh.
+/// Parent-child relationships can be optionally computed. Parent-child
+/// relationships can be used to create MeshTags on the refined mesh
+/// from MeshTags on the parent mesh.
+///
+/// @warning Using the default partitioner for a refined mesh, the
+/// refined mesh will **not** include ghosts cells (cells connected by
+/// facet to an owned cell) even if the parent mesh is ghosted.
+///
+/// @warning Passing `nullptr` for `partitioner`, the refined mesh will
+/// **not** have ghosts cells even if the parent mesh is ghosted. The
+/// possibility to not re-partition the refined mesh and include ghost
+/// cells in the refined mesh will be added in a future release.
+///
+/// @param[in] mesh Input mesh to be refined.
+/// @param[in] edges Indices of the edges that should be split in the
+/// refinement. If not provided (`std::nullopt`), uniform refinement is
+/// performed.
+/// @param[in] partitioner Partitioner to be used to distribute the
+/// refined mesh. If not callable, refined cells will be on the same
+/// process as the parent cell.
 /// @param[in] option Control the computation of parent facets, parent
-/// cells. If an option is unselected, an empty list is returned.
-/// @return New Mesh and optional parent cell index, parent facet
+/// cells. If an option is not selected, an empty list is returned.
+/// @return New mesh, and optional parent cell indices and parent facet
 /// indices.
 template <std::floating_point T>
 std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
            std::optional<std::vector<std::int8_t>>>
 refine(const mesh::Mesh<T>& mesh,
-       std::optional<std::span<const std::int32_t>> edges, bool redistribute,
-       mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet,
+       std::optional<std::span<const std::int32_t>> edges,
+       const mesh::CellPartitionFunction& partitioner
+       = mesh::create_cell_partitioner(mesh::GhostMode::none),
        Option option = Option::none)
 {
   auto topology = mesh.topology();
   assert(topology);
-
-  mesh::CellType cell_t = topology->cell_type();
-  if (!mesh::is_simplex(cell_t))
+  if (!mesh::is_simplex(topology->cell_type()))
     throw std::runtime_error("Refinement only defined for simplices");
-  bool oned = topology->cell_type() == mesh::CellType::interval;
+
   auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
-      = oned ? interval::compute_refinement_data(mesh, edges, option)
-             : plaza::compute_refinement_data(mesh, edges, option);
+      = (topology->cell_type() == mesh::CellType::interval)
+            ? interval::compute_refinement_data(mesh, edges, option)
+            : plaza::compute_refinement_data(mesh, edges, option);
 
-  mesh::Mesh<T> refined_mesh = impl::create_refined_mesh(
-      mesh, cell_adj, std::span<const T>(new_vertex_coords), xshape,
-      redistribute, ghost_mode);
+  mesh::Mesh<T> mesh1 = mesh::create_mesh(
+      mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(),
+      mesh.comm(), new_vertex_coords, xshape, partitioner);
 
-  // Report the number of refined cellse
+  // Report the number of refined cells
   const int D = topology->dim();
   const std::int64_t n0 = topology->index_map(D)->size_global();
-  const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global();
+  const std::int64_t n1 = mesh1.topology()->index_map(D)->size_global();
   spdlog::info(
       "Number of cells increased from {} to {} ({}% increase).", n0, n1,
       100.0 * (static_cast<double>(n1) / static_cast<double>(n0) - 1.0));
 
-  return {std::move(refined_mesh), std::move(parent_cell),
-          std::move(parent_facet)};
+  return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)};
 }
-
 } // namespace dolfinx::refinement
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index 275b40d7cc..ec550cb357 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -264,85 +264,48 @@ create_new_vertices(MPI_Comm comm,
           xshape};
 }
 
-/// Use vertex and topology data to partition new mesh across
-/// processes
-/// @param[in] old_mesh
-/// @param[in] cell_topology Topology of cells, (vertex indices)
-/// @param[in] new_coords New coordinates, row-major storage
-/// @param[in] xshape The shape of `new_coords`
-/// @param[in] redistribute Call graph partitioner if true
-/// @param[in] ghost_mode None or shared_facet
-/// @return New mesh
-template <std::floating_point T>
-mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
-                        const graph::AdjacencyList<std::int64_t>& cell_topology,
-                        std::span<const T> new_coords,
-                        std::array<std::size_t, 2> xshape, bool redistribute,
-                        mesh::GhostMode ghost_mode)
-{
-  if (redistribute)
-  {
-    return mesh::create_mesh(old_mesh.comm(), cell_topology.array(),
-                             old_mesh.geometry().cmap(), new_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    auto partitioner
-        = [](MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-             const std::vector<std::span<const std::int64_t>>& cell_topology)
-    {
-      std::int32_t mpi_rank = MPI::rank(comm);
-      std::int32_t num_cell_vertices
-          = mesh::num_cell_vertices(cell_types.front());
-      std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
-      std::vector<std::int32_t> destinations(num_cells, mpi_rank);
-      std::vector<std::int32_t> dest_offsets(num_cells + 1);
-      std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
-      return graph::AdjacencyList(std::move(destinations),
-                                  std::move(dest_offsets));
-    };
-
-    return mesh::create_mesh(old_mesh.comm(), old_mesh.comm(),
-                             cell_topology.array(), old_mesh.geometry().cmap(),
-                             old_mesh.comm(), new_coords, xshape, partitioner);
-  }
-}
-
-/// @brief Given an index map, add "n" extra indices at the end of local range
+/// @todo Improve docstring.
+///
+/// @brief Given an index map, add "n" extra indices at the end of local
+/// range.
 ///
-/// @note The returned global indices (local and ghosts) are adjust for the
-/// addition of new local indices.
-/// @param[in] map Index map
-/// @param[in] n Number of new local indices
-/// @return New global indices for both owned and ghosted indices in input
-/// index map.
+/// @note The returned global indices (local and ghosts) are adjust for
+/// the addition of new local indices.
+///
+/// @param[in] map Index map.
+/// @param[in] n Number of new local indices.
+/// @return New global indices for both owned and ghosted indices in
+/// input index map.
 std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
                                          std::int32_t n);
 
-/// @brief Transfer facet MeshTags from coarse mesh to refined mesh
-/// @note The refined mesh must not have been redistributed during
-/// refinement
-/// @note GhostMode must be GhostMode.none
-/// @param[in] tags0 Tags on the parent mesh
-/// @param[in] topology1 Refined mesh topology
+/// @brief Transfer facet MeshTags from coarse mesh to refined mesh.
+///
+/// @warning The refined mesh must not have been redistributed during
+/// refinement.
+///
+/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none.
+///
+/// @param[in] tags0 Tags on the parent mesh.
+/// @param[in] topology1 Refined mesh topology.
 /// @param[in] cell Parent cell of each cell in refined mesh
-/// @param[in] facet Local facets of parent in each cell in refined mesh
-/// @return (0) entities and (1) values on the refined topology
+/// @param[in] facet Local facets of parent in each cell in refined mesh.
+/// @return (0) entities and (1) values on the refined topology.
 std::array<std::vector<std::int32_t>, 2> transfer_facet_meshtag(
     const mesh::MeshTags<std::int32_t>& tags0, const mesh::Topology& topology1,
     std::span<const std::int32_t> cell, std::span<const std::int8_t> facet);
 
 /// @brief Transfer cell MeshTags from coarse mesh to refined mesh.
 ///
-/// @note The refined mesh must not have been redistributed during
+/// @warning The refined mesh must not have been redistributed during
 /// refinement.
-/// @note GhostMode must be GhostMode.none
 ///
-/// @param[in] tags0 Tags on the parent mesh
-/// @param[in] topology1 Refined mesh topology
-/// @param[in] parent_cell Parent cell of each cell in refined mesh
-/// @return (0) entities and (1) values on the refined topology
+/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none.
+///
+/// @param[in] tags0 Tags on the parent mesh.
+/// @param[in] topology1 Refined mesh topology.
+/// @param[in] parent_cell Parent cell of each cell in refined mesh.
+/// @return (0) entities and (1) values on the refined topology.
 std::array<std::vector<std::int32_t>, 2>
 transfer_cell_meshtag(const mesh::MeshTags<std::int32_t>& tags0,
                       const mesh::Topology& topology1,
diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp
index 763418ac4b..1253fe00b0 100644
--- a/cpp/test/mesh/refinement/interval.cpp
+++ b/cpp/test/mesh/refinement/interval.cpp
@@ -73,8 +73,7 @@ TEMPLATE_TEST_CASE("Interval uniform refinement",
 
   // TODO: parent_facet
   auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
-      mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
-      refinement::Option::parent_cell);
+      mesh, std::nullopt, nullptr, refinement::Option::parent_cell);
 
   std::vector<T> expected_x = {
       /* v_0 */ 0.0, 0.0, 0.0,
@@ -114,7 +113,8 @@ TEMPLATE_TEST_CASE("Interval adaptive refinement",
   std::vector<std::int32_t> edges{1};
   // TODO: parent_facet
   auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
-      mesh, std::span(edges), false, mesh::GhostMode::shared_facet,
+      mesh, std::span(edges),
+      mesh::create_cell_partitioner(mesh::GhostMode::shared_facet),
       refinement::Option::parent_cell);
 
   std::vector<T> expected_x = {
@@ -190,15 +190,13 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)",
   mesh::Mesh<T> mesh = create_mesh();
   mesh.topology()->create_connectivity(1, 0);
 
-  // TODO: parent_facet
   auto [refined_mesh, parent_edges, parent_facet] = refinement::refine(
-      mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
-      refinement::Option::parent_cell);
+      mesh, std::nullopt, nullptr, refinement::Option::parent_cell);
 
   T rank_d = static_cast<T>(rank);
   T comm_size_d = static_cast<T>(comm_size);
 
-  auto x = refined_mesh.geometry().x();
+  std::span x = refined_mesh.geometry().x();
   std::ranges::sort(x);
   std::vector<T> expected_x
       = {rank_d / comm_size_d,
diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp
index 954879abb1..ff2a67979d 100644
--- a/cpp/test/mesh/refinement/rectangle.cpp
+++ b/cpp/test/mesh/refinement/rectangle.cpp
@@ -99,9 +99,9 @@ plotter.show()
   // plaza requires the edges to be pre initialized!
   mesh.topology()->create_entities(1);
 
-  auto [mesh_fine, parent_cell, parent_facet]
-      = refinement::refine(mesh, std::nullopt, false, mesh::GhostMode::none,
-                           refinement::Option::parent_cell_and_facet);
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      mesh, std::nullopt, mesh::create_cell_partitioner(mesh::GhostMode::none),
+      refinement::Option::parent_cell_and_facet);
 
   // vertex layout:
   // 8---7---5
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 28b9eb759e..2203ca8d96 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -276,11 +276,11 @@ class Mesh:
     _geometry: Geometry
     _ufl_domain: typing.Optional[ufl.Mesh]
 
-    def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
+    def __init__(self, msh, domain: typing.Optional[ufl.Mesh]):
         """Initialize mesh from a C++ mesh.
 
         Args:
-            mesh: A C++ mesh object.
+            msh: A C++ mesh object.
             domain: A UFL domain.
 
         Note:
@@ -289,7 +289,7 @@ def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
             This class is combined with different base classes that
             depend on the scalar type used in the Mesh.
         """
-        self._cpp_object = mesh
+        self._cpp_object = msh
         self._topology = Topology(self._cpp_object.topology)
         self._geometry = Geometry(self._cpp_object.geometry)
         self._ufl_domain = domain
@@ -439,25 +439,25 @@ def compute_incident_entities(
     return _cpp.mesh.compute_incident_entities(topology._cpp_object, entities, d0, d1)
 
 
-def compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]):
+def compute_midpoints(msh: Mesh, dim: int, entities: npt.NDArray[np.int32]):
     """Compute the midpoints of a set of mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entities to consider.
         entities: Indices of entities in ``mesh`` to consider.
 
     Returns:
         Midpoints of the entities, shape ``(num_entities, 3)``.
     """
-    return _cpp.mesh.compute_midpoints(mesh._cpp_object, dim, entities)
+    return _cpp.mesh.compute_midpoints(msh._cpp_object, dim, entities)
 
 
-def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
+def locate_entities(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
     """Compute mesh entities satisfying a geometric marking function.
 
     Args:
-        mesh: Mesh to locate entities on.
+        msh: Mesh to locate entities on.
         dim: Topological dimension of the mesh entities to consider.
         marker: A function that takes an array of points ``x`` with
             shape ``(gdim, num_points)`` and returns an array of
@@ -467,10 +467,10 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray
     Returns:
         Indices (local to the process) of marked mesh entities.
     """
-    return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker)
+    return _cpp.mesh.locate_entities(msh._cpp_object, dim, marker)
 
 
-def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
+def locate_entities_boundary(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
     """Compute mesh entities that are connected to an owned boundary
     facet and satisfy a geometric marking function.
 
@@ -484,7 +484,7 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n
     parallel communication.
 
     Args:
-        mesh: Mesh to locate boundary entities on.
+        msh: Mesh to locate boundary entities on.
         dim: Topological dimension of the mesh entities to consider
         marker: Function that takes an array of points ``x`` with shape
             ``(gdim, num_points)`` and returns an array of booleans of
@@ -494,12 +494,12 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n
     Returns:
         Indices (local to the process) of marked mesh entities.
     """
-    return _cpp.mesh.locate_entities_boundary(mesh._cpp_object, dim, marker)
+    return _cpp.mesh.locate_entities_boundary(msh._cpp_object, dim, marker)
 
 
 def transfer_meshtag(
     meshtag: MeshTags,
-    mesh1: Mesh,
+    msh1: Mesh,
     parent_cell: npt.NDArray[np.int32],
     parent_facet: typing.Optional[npt.NDArray[np.int8]] = None,
 ) -> MeshTags:
@@ -507,7 +507,7 @@ def transfer_meshtag(
 
     Args:
         meshtag: Mesh tags on the coarse, parent mesh.
-        mesh1: The refined mesh.
+        msh1: The refined mesh.
         parent_cell: Index of the parent cell for each cell in the
             refined mesh.
         parent_facet: Index of the local parent facet for each cell
@@ -519,13 +519,13 @@ def transfer_meshtag(
     """
     if meshtag.dim == meshtag.topology.dim:
         mt = _cpp.refinement.transfer_cell_meshtag(
-            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell
+            meshtag._cpp_object, msh1.topology._cpp_object, parent_cell
         )
         return MeshTags(mt)
     elif meshtag.dim == meshtag.topology.dim - 1:
         assert parent_facet is not None
         mt = _cpp.refinement.transfer_facet_meshtag(
-            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell, parent_facet
+            meshtag._cpp_object, msh1.topology._cpp_object, parent_cell, parent_facet
         )
         return MeshTags(mt)
     else:
@@ -533,32 +533,47 @@ def transfer_meshtag(
 
 
 def refine(
-    mesh: Mesh,
+    msh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
-    redistribute: bool = True,
-    ghost_mode: GhostMode = GhostMode.shared_facet,
+    partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none),
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
 
+    Passing ``None`` for ``partitioner``, refined cells will be on the
+    same process as the parent cell.
+
+    Note:
+        Using the default partitioner for the refined mesh, the refined
+        mesh will **not** include ghosts cells (cells connected by facet
+        to an owned cells) even if the parent mesh is ghosted.
+
+    Warning:
+        Passing ``None`` for ``partitioner``, the refined mesh will
+        **not** have ghosts cells even if the parent mesh has ghost
+        cells. The possibility to not re-partition the refined mesh and
+        include ghost cells in the refined mesh will be added in a
+        future release.
+
     Args:
-        mesh: Mesh from which to create the refined mesh.
+        msh: Mesh from which to create the refined mesh.
         edges: Indices of edges to split during refinement. If ``None``,
             mesh refinement is uniform.
-        redistribute:
-            Refined mesh is re-partitioned if ``True``.
-        ghost_mode: Ghost mode to use for the refined mesh.
-        option: Controls whether parent cells and/or parent facets are computed.
-
+        partitioner: Partitioner to distribute the refined mesh. If
+            ``None`` no redistribution is performed, i.e. refined cells
+            remain on the same process as the parent cell.
+        option: Controls whether parent cells and/or parent facets are
+            computed.
 
     Returns:
        Refined mesh, (optional) parent cells, (optional) parent facets
     """
     mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
-        mesh._cpp_object, edges, redistribute, ghost_mode, option
+        msh._cpp_object, edges, partitioner, option
     )
-    # Create new ufl domain as it will carry a reference to the C++ mesh in the ufl_cargo
-    ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element())  # type: ignore
+    # Create new ufl domain as it will carry a reference to the C++ mesh
+    # in the ufl_cargo
+    ufl_domain = ufl.Mesh(msh._ufl_domain.ufl_coordinate_element())  # type: ignore
     return Mesh(mesh1, ufl_domain), parent_cell, parent_facet
 
 
@@ -582,8 +597,8 @@ def create_mesh(
         x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, gdim)``.
         e: UFL mesh. The mesh scalar type is determined by the scalar
             type of ``e``.
-        partitioner: Function that computes the parallel distribution of
-            cells across MPI ranks.
+        partitioner: Function that determines the parallel distribution
+            of cells across MPI ranks.
 
     Note:
         If required, the coordinates ``x`` will be cast to the same
@@ -593,7 +608,7 @@ def create_mesh(
         A mesh.
     """
     if partitioner is None and comm.size > 1:
-        partitioner = _cpp.mesh.create_cell_partitioner(GhostMode.none)
+        partitioner = create_cell_partitioner(GhostMode.none)
 
     x = np.asarray(x, order="C")
     if x.ndim == 1:
@@ -635,9 +650,9 @@ def create_mesh(
 
     x = np.asarray(x, dtype=dtype, order="C")
     cells = np.asarray(cells, dtype=np.int64, order="C")
-    mesh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner)
+    msh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner)
 
-    return Mesh(mesh, domain)
+    return Mesh(msh, domain)
 
 
 def create_submesh(
@@ -646,7 +661,7 @@ def create_submesh(
     """Create a mesh with specified entities from another mesh.
 
     Args:
-        mesh: Mesh to create the sub-mesh from.
+        msh: Mesh to create the sub-mesh from.
         dim: Topological dimension of the entities in ``msh`` to include
             in the sub-mesh.
         entities: Indices of entities in ``msh`` to include in the
@@ -675,7 +690,7 @@ def create_submesh(
 
 
 def meshtags(
-    mesh: Mesh,
+    msh: Mesh,
     dim: int,
     entities: npt.NDArray[np.int32],
     values: typing.Union[np.ndarray, int, float],
@@ -683,7 +698,7 @@ def meshtags(
     """Create a MeshTags object that associates data with a subset of mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entity.
         entities: Indices(local to process) of entities to associate
             values with . The array must be sorted and must not contain
@@ -716,19 +731,19 @@ def meshtags(
         raise NotImplementedError(f"Type {values.dtype} not supported.")
 
     return MeshTags(
-        ftype(mesh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values)
+        ftype(msh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values)
     )
 
 
 def meshtags_from_entities(
-    mesh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any]
+    msh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any]
 ):
     """Create a :class:dolfinx.mesh.MeshTags` object that associates
     data with a subset of mesh entities, where the entities are defined
     by their vertices.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entity.
         entities: Entities to associated values with, with entities
             defined by their vertices.
@@ -748,7 +763,7 @@ def meshtags_from_entities(
     elif isinstance(values, float):
         values = np.full(entities.num_nodes, values, dtype=np.double)
     values = np.asarray(values)
-    return MeshTags(_cpp.mesh.create_meshtags(mesh.topology._cpp_object, dim, entities, values))
+    return MeshTags(_cpp.mesh.create_meshtags(msh.topology._cpp_object, dim, entities, values))
 
 
 def create_interval(
@@ -779,12 +794,13 @@ def create_interval(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", "interval", 1, shape=(1,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner)
+        msh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner)
+        msh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_interval(
@@ -847,12 +863,13 @@ def create_rectangle(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(2,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal)
+        msh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal)
+        msh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_square(
@@ -925,12 +942,13 @@ def create_box(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(3,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner)
+        msh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner)
+        msh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_cube(
@@ -973,12 +991,12 @@ def create_unit_cube(
 
 
 def entities_to_geometry(
-    mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False
+    msh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False
 ) -> npt.NDArray[np.int32]:
     """Compute the geometric DOFs associated with the closure of the given mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the entities of interest.
         entities: Entity indices (local to the process).
         permute: Permute the DOFs such that they are consistent with the
@@ -989,7 +1007,7 @@ def entities_to_geometry(
         The geometric DOFs associated with the closure of the entities
         in `entities`.
     """
-    return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, permute)
+    return _cpp.mesh.entities_to_geometry(msh._cpp_object, dim, entities, permute)
 
 
 def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]:
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index f8c3679974..5cf66c9c6a 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -4,11 +4,11 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
+#include "mesh.h"
 #include "MPICommWrapper.h"
 #include "array.h"
 #include "caster_mpi.h"
 #include "numpy_dtype.h"
-#include <cfloat>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/fem/CoordinateElement.h>
 #include <dolfinx/fem/ElementDofLayout.h>
@@ -21,7 +21,6 @@
 #include <dolfinx/mesh/graphbuild.h>
 #include <dolfinx/mesh/topologycomputation.h>
 #include <dolfinx/mesh/utils.h>
-#include <iostream>
 #include <memory>
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
@@ -36,49 +35,10 @@
 
 namespace nb = nanobind;
 
-namespace
-{
-/// Wrap a Python graph partitioning function as a C++ function
-template <typename Functor>
-auto create_partitioner_cpp(Functor p)
+namespace dolfinx_wrappers
 {
-  return [p](MPI_Comm comm, int nparts,
-             const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
-             bool ghosting)
-  {
-    return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph,
-             ghosting);
-  };
-}
-
-/// Wrap a C++ cell partitioning function as a Python function
-template <typename Functor>
-auto create_cell_partitioner_py(Functor p)
+namespace part::impl
 {
-  return [p](dolfinx_wrappers::MPICommWrapper comm, int n,
-             const std::vector<dolfinx::mesh::CellType>& cell_types,
-             std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb)
-  {
-    std::vector<std::span<const std::int64_t>> cells;
-    std::ranges::transform(
-        cells_nb, std::back_inserter(cells), [](auto c)
-        { return std::span<const std::int64_t>(c.data(), c.size()); });
-    return p(comm.get(), n, cell_types, cells);
-  };
-}
-
-using PythonCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        dolfinx_wrappers::MPICommWrapper, int,
-        const std::vector<dolfinx::mesh::CellType>&,
-        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
-
-using CppCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
-        const std::vector<std::span<const std::int64_t>>&)>;
-
-/// Wrap a Python cell graph partitioning function as a C++ function
 CppCellPartitionFunction
 create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
 {
@@ -103,10 +63,8 @@ create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
   else
     return nullptr;
 }
-} // namespace
+} // namespace part::impl
 
-namespace dolfinx_wrappers
-{
 template <typename T>
 void declare_meshtags(nb::module_& m, std::string type)
 {
@@ -260,7 +218,8 @@ void declare_mesh(nb::module_& m, std::string type)
           "__init__",
           [](dolfinx::mesh::Mesh<T>* mesh, MPICommWrapper comm,
              std::shared_ptr<dolfinx::mesh::Topology> topology,
-             dolfinx::mesh::Geometry<T>& geometry) {
+             dolfinx::mesh::Geometry<T>& geometry)
+          {
             new (mesh) dolfinx::mesh::Mesh<T>(comm.get(), topology, geometry);
           },
           nb::arg("comm"), nb::arg("topology"), nb::arg("geometry"))
@@ -279,11 +238,12 @@ void declare_mesh(nb::module_& m, std::string type)
   m.def(
       create_interval.c_str(),
       [](MPICommWrapper comm, std::int64_t n, std::array<T, 2> p,
-         dolfinx::mesh::GhostMode ghost_mode,
-         const PythonCellPartitionFunction& part)
+         dolfinx::mesh::GhostMode mode,
+         const part::impl::PythonCellPartitionFunction& part)
       {
         return dolfinx::mesh::create_interval<T>(
-            comm.get(), n, p, ghost_mode, create_cell_partitioner_cpp(part));
+            comm.get(), n, p, mode,
+            part::impl::create_cell_partitioner_cpp(part));
       },
       nb::arg("comm"), nb::arg("n"), nb::arg("p"), nb::arg("ghost_mode"),
       nb::arg("partitioner").none());
@@ -293,12 +253,12 @@ void declare_mesh(nb::module_& m, std::string type)
       create_rectangle.c_str(),
       [](MPICommWrapper comm, std::array<std::array<T, 2>, 2> p,
          std::array<std::int64_t, 2> n, dolfinx::mesh::CellType celltype,
-         const PythonCellPartitionFunction& part,
+         const part::impl::PythonCellPartitionFunction& part,
          dolfinx::mesh::DiagonalType diagonal)
       {
         return dolfinx::mesh::create_rectangle<T>(
-            comm.get(), p, n, celltype, create_cell_partitioner_cpp(part),
-            diagonal);
+            comm.get(), p, n, celltype,
+            part::impl::create_cell_partitioner_cpp(part), diagonal);
       },
       nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"),
       nb::arg("partitioner").none(), nb::arg("diagonal"));
@@ -308,11 +268,12 @@ void declare_mesh(nb::module_& m, std::string type)
       create_box.c_str(),
       [](MPICommWrapper comm, std::array<std::array<T, 3>, 2> p,
          std::array<std::int64_t, 3> n, dolfinx::mesh::CellType celltype,
-         const PythonCellPartitionFunction& part)
+         const part::impl::PythonCellPartitionFunction& part)
       {
         MPI_Comm _comm = comm.get();
-        return dolfinx::mesh::create_box<T>(_comm, _comm, p, n, celltype,
-                                            create_cell_partitioner_cpp(part));
+        return dolfinx::mesh::create_box<T>(
+            _comm, _comm, p, n, celltype,
+            part::impl::create_cell_partitioner_cpp(part));
       },
       nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"),
       nb::arg("partitioner").none());
@@ -323,7 +284,7 @@ void declare_mesh(nb::module_& m, std::string type)
                                          nb::c_contig>>& cells_nb,
            const std::vector<dolfinx::fem::CoordinateElement<T>>& elements,
            nb::ndarray<const T, nb::c_contig> x,
-           const PythonCellPartitionFunction& p)
+           const part::impl::PythonCellPartitionFunction& p)
         {
           std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1);
 
@@ -365,7 +326,7 @@ void declare_mesh(nb::module_& m, std::string type)
          nb::ndarray<const std::int64_t, nb::ndim<2>, nb::c_contig> cells,
          const dolfinx::fem::CoordinateElement<T>& element,
          nb::ndarray<const T, nb::c_contig> x,
-         const PythonCellPartitionFunction& p)
+         const part::impl::PythonCellPartitionFunction& p)
       {
         std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1);
         if (p)
@@ -599,7 +560,8 @@ void mesh(nb::module_& m)
   m.def(
       "compute_entities",
       [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology, int dim,
-         int index) {
+         int index)
+      {
         return dolfinx::mesh::compute_entities(comm.get(), topology, dim,
                                                index);
       },
@@ -721,8 +683,6 @@ void mesh(nb::module_& m)
               ghost_owners_span, boundary_vertices_span);
         });
 
-  // dolfinx::mesh::MeshTags
-
   declare_meshtags<std::int8_t>(m, "int8");
   declare_meshtags<std::int32_t>(m, "int32");
   declare_meshtags<std::int64_t>(m, "int64");
@@ -733,10 +693,11 @@ void mesh(nb::module_& m)
 
   m.def(
       "create_cell_partitioner",
-      [](dolfinx::mesh::GhostMode gm) -> PythonCellPartitionFunction
+      [](dolfinx::mesh::GhostMode mode)
+          -> part::impl::PythonCellPartitionFunction
       {
-        return create_cell_partitioner_py(
-            dolfinx::mesh::create_cell_partitioner(gm));
+        return part::impl::create_cell_partitioner_py(
+            dolfinx::mesh::create_cell_partitioner(mode));
       },
       "Create default cell partitioner.");
   m.def(
@@ -746,11 +707,12 @@ void mesh(nb::module_& m)
              const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
              bool ghosting)>
              part,
-         dolfinx::mesh::GhostMode ghost_mode) -> PythonCellPartitionFunction
+         dolfinx::mesh::GhostMode mode)
+          -> part::impl::PythonCellPartitionFunction
       {
-        return create_cell_partitioner_py(
+        return part::impl::create_cell_partitioner_py(
             dolfinx::mesh::create_cell_partitioner(
-                ghost_mode, create_partitioner_cpp(part)));
+                mode, part::impl::create_partitioner_cpp(part)));
       },
       nb::arg("part"), nb::arg("ghost_mode") = dolfinx::mesh::GhostMode::none,
       "Create a cell partitioner from a graph partitioning function.");
diff --git a/python/dolfinx/wrappers/mesh.h b/python/dolfinx/wrappers/mesh.h
new file mode 100644
index 0000000000..5471dc1ab2
--- /dev/null
+++ b/python/dolfinx/wrappers/mesh.h
@@ -0,0 +1,59 @@
+// Copyright (C) 2017-2024 Chris N. Richardson and Garth N. Wells
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include "MPICommWrapper.h"
+#include <dolfinx/mesh/cell_types.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/vector.h>
+
+namespace nb = nanobind;
+
+namespace dolfinx_wrappers::part::impl
+{
+/// Wrap a Python graph partitioning function as a C++ function
+template <typename Functor>
+auto create_partitioner_cpp(Functor p)
+{
+  return [p](MPI_Comm comm, int nparts,
+             const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
+             bool ghosting)
+  {
+    return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph,
+             ghosting);
+  };
+}
+
+/// Wrap a C++ cell partitioning function as a Python function
+template <typename Functor>
+auto create_cell_partitioner_py(Functor p)
+{
+  return [p](dolfinx_wrappers::MPICommWrapper comm, int n,
+             const std::vector<dolfinx::mesh::CellType>& cell_types,
+             std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb)
+  {
+    std::vector<std::span<const std::int64_t>> cells;
+    std::ranges::transform(
+        cells_nb, std::back_inserter(cells), [](auto c)
+        { return std::span<const std::int64_t>(c.data(), c.size()); });
+    return p(comm.get(), n, cell_types, cells);
+  };
+}
+
+using PythonCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        dolfinx_wrappers::MPICommWrapper, int,
+        const std::vector<dolfinx::mesh::CellType>&,
+        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
+
+using CppCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
+        const std::vector<std::span<const std::int64_t>>&)>;
+
+/// Wrap a Python cell graph partitioning function as a C++ function
+CppCellPartitionFunction
+create_cell_partitioner_cpp(const PythonCellPartitionFunction& p);
+} // namespace dolfinx_wrappers::part::impl
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 37ca393799..d9b86fcf78 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -1,21 +1,24 @@
-// Copyright (C) 2018 Chris N. Richardson and Garth N. Wells
+// Copyright (C) 2018-2024 Chris N. Richardson and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
+#include "MPICommWrapper.h"
 #include "array.h"
+#include "mesh.h"
 #include <concepts>
 #include <cstdint>
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/MeshTags.h>
 #include <dolfinx/refinement/interval.h>
 #include <dolfinx/refinement/option.h>
-#include <dolfinx/refinement/plaza.h>
 #include <dolfinx/refinement/refine.h>
 #include <dolfinx/refinement/utils.h>
+#include <functional>
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
+#include <nanobind/stl/function.h>
 #include <nanobind/stl/optional.h>
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/tuple.h>
@@ -25,11 +28,10 @@
 
 namespace nb = nanobind;
 
-namespace dolfinx_wrappers
+namespace
 {
-
 template <std::floating_point T>
-void export_refinement_with_variable_mesh_type(nb::module_& m)
+void export_refinement(nb::module_& m)
 {
   m.def(
       "refine",
@@ -37,7 +39,9 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
          std::optional<
              nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
              edges,
-         bool redistribute, dolfinx::mesh::GhostMode ghost_mode,
+         std::optional<
+             dolfinx_wrappers::part::impl::PythonCellPartitionFunction>
+             partitioner,
          dolfinx::refinement::Option option)
       {
         std::optional<std::span<const std::int32_t>> cpp_edges(std::nullopt);
@@ -47,30 +51,45 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
               std::span(edges.value().data(), edges.value().size()));
         }
 
+        dolfinx_wrappers::part::impl::CppCellPartitionFunction cpp_partitioner
+            = partitioner.has_value()
+                  ? dolfinx_wrappers::part::impl::create_cell_partitioner_cpp(
+                        partitioner.value())
+                  : nullptr;
         auto [mesh1, cell, facet] = dolfinx::refinement::refine(
-            mesh, cpp_edges, redistribute, ghost_mode, option);
+            mesh, cpp_edges, cpp_partitioner, option);
 
         std::optional<nb::ndarray<std::int32_t, nb::numpy>> python_cell(
             std::nullopt);
         if (cell.has_value())
-          python_cell.emplace(as_nbarray(std::move(cell.value())));
+        {
+          python_cell.emplace(
+              dolfinx_wrappers::as_nbarray(std::move(cell.value())));
+        }
 
         std::optional<nb::ndarray<std::int8_t, nb::numpy>> python_facet(
             std::nullopt);
         if (facet.has_value())
-          python_facet.emplace(as_nbarray(std::move(facet.value())));
+        {
+          python_facet.emplace(
+              dolfinx_wrappers::as_nbarray(std::move(facet.value())));
+        }
 
         return std::tuple{std::move(mesh1), std::move(python_cell),
                           std::move(python_facet)};
       },
-      nb::arg("mesh"), nb::arg("edges") = nb::none(), nb::arg("redistribute"),
-      nb::arg("ghost_mode"), nb::arg("option"));
+      nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(),
+      nb::arg("option"));
 }
+} // namespace
+
+namespace dolfinx_wrappers
+{
 
 void refinement(nb::module_& m)
 {
-  export_refinement_with_variable_mesh_type<float>(m);
-  export_refinement_with_variable_mesh_type<double>(m);
+  export_refinement<float>(m);
+  export_refinement<double>(m);
 
   nb::enum_<dolfinx::refinement::Option>(m, "RefinementOption")
       .value("none", dolfinx::refinement::Option::none)
diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py
index 2877964484..27040fb0b9 100644
--- a/python/test/unit/refinement/test_interval.py
+++ b/python/test/unit/refinement/test_interval.py
@@ -20,16 +20,14 @@
     "ghost_mode_refined",
     [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet],
 )
-@pytest.mark.parametrize("redistribute", [True, False])
 @pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell])
-def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option):
+def test_refine_interval(n, ghost_mode, ghost_mode_refined, option):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
     msh_refined, edges, vertices = mesh.refine(
         msh,
-        redistribute=redistribute,
-        ghost_mode=ghost_mode_refined,
         option=option,
     )
+    # TODO: add create_cell_partitioner(ghost_mode) when works
 
     # vertex count
     assert msh_refined.topology.index_map(0).size_global == 2 * n + 1
@@ -52,16 +50,14 @@ def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option
     "ghost_mode_refined",
     [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet],
 )
-@pytest.mark.parametrize("redistribute", [True, False])
-def test_refine_interval_adaptive(n, ghost_mode, redistribute, ghost_mode_refined):
+def test_refine_interval_adaptive(n, ghost_mode, ghost_mode_refined):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
     msh_refined, edges, vertices = mesh.refine(
         msh,
         np.arange(10, dtype=np.int32),
-        redistribute=redistribute,
-        ghost_mode=ghost_mode_refined,
         option=mesh.RefinementOption.parent_cell,
     )
+    # TODO: add create_cell_partitioner(ghost_mode) when works
 
     # vertex count
     assert msh_refined.topology.index_map(0).size_global == n + 1 + 10 * MPI.COMM_WORLD.size
diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 6aaa4090df..0a8ee67b79 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -11,6 +11,7 @@
 from numpy import isclose
 
 import ufl
+from dolfinx.cpp.mesh import create_cell_partitioner
 from dolfinx.fem import assemble_matrix, form, functionspace
 from dolfinx.mesh import (
     CellType,
@@ -32,7 +33,7 @@ def test_refine_create_unit_square():
     """Refine mesh of unit square."""
     mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none)
     mesh.topology.create_entities(1)
-    mesh_refined, _, _ = refine(mesh, redistribute=False)
+    mesh_refined, _, _ = refine(mesh)
     assert mesh_refined.topology.index_map(0).size_global == 165
     assert mesh_refined.topology.index_map(2).size_global == 280
 
@@ -46,7 +47,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute):
     """Refine mesh of unit cube."""
     mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode)
     mesh.topology.create_entities(1)
-    mesh, _, _ = refine(mesh, redistribute=redistribute)
+    mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode))
     assert mesh.topology.index_map(0).size_global == 3135
     assert mesh.topology.index_map(3).size_global == 15120
 
@@ -58,8 +59,7 @@ def test_refine_create_form():
     """Check that forms can be assembled on refined mesh"""
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)
     mesh.topology.create_entities(1)
-    mesh, _, _ = refine(mesh, redistribute=True)
-
+    mesh, _, _ = refine(mesh)
     V = functionspace(mesh, ("Lagrange", 1))
 
     # Define variational problem
@@ -71,20 +71,20 @@ def test_refine_create_form():
 
 def test_sub_refine():
     """Test that refinement of a subset of edges works"""
-    mesh = create_unit_square(
+    msh = create_unit_square(
         MPI.COMM_WORLD, 3, 4, diagonal=DiagonalType.left, ghost_mode=GhostMode.none
     )
-    mesh.topology.create_entities(1)
+    msh.topology.create_entities(1)
 
     def left_corner_edge(x, tol=1e-7):
         return isclose(x[0], 0) & (x[1] < 1 / 4 + tol)
 
-    edges = locate_entities_boundary(mesh, 1, left_corner_edge)
-    if MPI.COMM_WORLD.size == 0:
+    edges = locate_entities_boundary(msh, 1, left_corner_edge)
+    if MPI.COMM_WORLD.size == 1:
         assert edges == 1
 
-    mesh2, _, _ = refine(mesh, edges, redistribute=False)
-    assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global
+    msh1, _, _ = refine(msh, edges)
+    assert msh.topology.index_map(2).size_global + 3 == msh1.topology.index_map(2).size_global
 
 
 def test_refine_from_cells():
@@ -100,12 +100,12 @@ def left_side(x, tol=1e-16):
         return x[0] <= 0.5 + tol
 
     cells = locate_entities(msh, msh.topology.dim, left_side)
-    if MPI.COMM_WORLD.size == 0:
-        assert cells.__len__() == Nx * Ny
+    if MPI.COMM_WORLD.size == 1:
+        assert len(cells) == Nx * Ny
     edges = compute_incident_entities(msh.topology, cells, 2, 1)
-    if MPI.COMM_WORLD.size == 0:
-        assert edges.__len__() == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny
-    mesh2, _, _ = refine(msh, edges, redistribute=True)
+    if MPI.COMM_WORLD.size == 1:
+        assert len(edges) == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny
+    mesh2, _, _ = refine(msh, edges)
 
     num_cells_global = mesh2.topology.index_map(2).size_global
     actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny
@@ -116,61 +116,58 @@ def left_side(x, tol=1e-16):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(
-            mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet
-        ),
-        lambda mesh: refine(
-            mesh,
-            np.arange(mesh.topology.index_map(1).size_local),
-            redistribute=False,
+        lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
+        lambda msh: refine(
+            msh,
+            edges=np.arange(msh.topology.index_map(1).size_local),
+            partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],
 )
 def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
     if tdim == 3:
-        mesh = create_unit_cube(
+        msh = create_unit_cube(
             MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none
         )
     else:
-        mesh = create_unit_square(
-            MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none
-        )
-    mesh.topology.create_entities(tdim - 1)
-    mesh.topology.create_connectivity(tdim - 1, tdim)
-    mesh.topology.create_entities(1)
-    f_to_c = mesh.topology.connectivity(tdim - 1, tdim)
+        msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none)
+    msh.topology.create_entities(tdim - 1)
+    msh.topology.create_connectivity(tdim - 1, tdim)
+    msh.topology.create_entities(1)
+    f_to_c = msh.topology.connectivity(tdim - 1, tdim)
     facet_indices = []
-    for f in range(mesh.topology.index_map(tdim - 1).size_local):
+    for f in range(msh.topology.index_map(tdim - 1).size_local):
         if len(f_to_c.links(f)) == 1:
             facet_indices += [f]
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim - 1,
         np.array(facet_indices, dtype=np.int32),
         np.arange(len(facet_indices), dtype=np.int32),
     )
 
-    fine_mesh, parent_cell, parent_facet = refine_plaza_wrapper(mesh)
-    fine_mesh.topology.create_entities(tdim - 1)
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet)
+    msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh)
+
+    msh1.topology.create_entities(tdim - 1)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet)
     assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices)
 
     # New tags should be on facets with one cell (i.e. exterior)
-    fine_mesh.topology.create_connectivity(tdim - 1, tdim)
-    new_f_to_c = fine_mesh.topology.connectivity(tdim - 1, tdim)
+    msh1.topology.create_connectivity(tdim - 1, tdim)
+    new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim)
     for f in new_meshtag.indices:
         assert len(new_f_to_c.links(f)) == 1
 
     # Now mark all facets (including internal)
-    facet_indices = np.arange(mesh.topology.index_map(tdim - 1).size_local)
+    facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local)
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim - 1,
         np.array(facet_indices, dtype=np.int32),
         np.arange(len(facet_indices), dtype=np.int32),
     )
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet)
     assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices)
 
 
@@ -178,44 +175,40 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(
-            mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet
-        ),
-        lambda mesh: refine(
-            mesh,
-            np.arange(mesh.topology.index_map(1).size_local),
-            redistribute=False,
+        lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
+        lambda msh: refine(
+            msh,
+            np.arange(msh.topology.index_map(1).size_local),
+            partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],
 )
 def test_refine_cell_meshtag(tdim, refine_plaza_wrapper):
     if tdim == 3:
-        mesh = create_unit_cube(
+        msh = create_unit_cube(
             MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none
         )
     else:
-        mesh = create_unit_square(
-            MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none
-        )
+        msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none)
 
-    mesh.topology.create_entities(1)
-    cell_indices = np.arange(mesh.topology.index_map(tdim).size_local)
+    msh.topology.create_entities(1)
+    cell_indices = np.arange(msh.topology.index_map(tdim).size_local)
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim,
         np.array(cell_indices, dtype=np.int32),
         np.arange(len(cell_indices), dtype=np.int32),
     )
 
-    fine_mesh, parent_cell, _ = refine_plaza_wrapper(mesh)
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell)
+    msh1, parent_cell, _ = refine_plaza_wrapper(msh)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell)
     assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values)
     assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices)
 
 
 def test_refine_ufl_cargo():
-    mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3)
-    mesh.topology.create_entities(1)
-    refined_mesh, _, _ = refine(mesh, redistribute=False)
-    assert refined_mesh.ufl_domain().ufl_cargo() != mesh.ufl_domain().ufl_cargo()
+    msh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3)
+    msh.topology.create_entities(1)
+    msh1, _, _ = refine(msh)
+    assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo()