diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt
index 5254b1dd066..9560e45c823 100644
--- a/cpp/dolfinx/refinement/CMakeLists.txt
+++ b/cpp/dolfinx/refinement/CMakeLists.txt
@@ -4,6 +4,7 @@ set(HEADERS_refinement
     ${CMAKE_CURRENT_SOURCE_DIR}/plaza.h
     ${CMAKE_CURRENT_SOURCE_DIR}/refine.h
     ${CMAKE_CURRENT_SOURCE_DIR}/utils.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/option.h
     PARENT_SCOPE
 )
 
diff --git a/cpp/dolfinx/refinement/interval.h b/cpp/dolfinx/refinement/interval.h
index 50bb290d22a..5348260fd3a 100644
--- a/cpp/dolfinx/refinement/interval.h
+++ b/cpp/dolfinx/refinement/interval.h
@@ -19,23 +19,33 @@
 #include <stdexcept>
 #include <vector>
 
-namespace dolfinx::refinement
-{
-namespace impl
+#include "dolfinx/refinement/option.h"
+#include "dolfinx/refinement/utils.h"
+
+namespace dolfinx::refinement::interval
 {
 /// @brief Refine with markers returning new mesh data.
 ///
-/// @param[in] mesh Input mesh to be refined.
-/// @param[in] cells Indices of the cells that are marked for refinement.
-///
+/// @param[in] mesh Input mesh to be refined
+/// @param[in] cells Indices of the cells that are marked for refinement
+/// @param[in] option Refinement option indicating if parent cells
+/// and/or facets are to be computed.
 /// @return New mesh data: cell topology, vertex coordinates and parent
 /// cell indices.
 template <std::floating_point T>
 std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
-           std::array<std::size_t, 2>, std::vector<std::int32_t>>
-compute_interval_refinement(const mesh::Mesh<T>& mesh,
-                            std::optional<std::span<const std::int32_t>> cells)
+           std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
+           std::optional<std::vector<std::int8_t>>>
+compute_refinement_data(const mesh::Mesh<T>& mesh,
+                        std::optional<std::span<const std::int32_t>> cells,
+                        Option option)
 {
+  bool compute_parent_facet = option_parent_facet(option);
+  bool compute_parent_cell = option_parent_cell(option);
+
+  if (compute_parent_facet)
+    throw std::runtime_error("Parent facet computation not yet supported!");
+
   auto topology = mesh.topology();
   assert(topology);
   assert(topology->dim() == 1);
@@ -124,8 +134,14 @@ compute_interval_refinement(const mesh::Mesh<T>& mesh,
 
   std::vector<std::int64_t> cell_topology;
   cell_topology.reserve(refined_cell_count * 2);
-  std::vector<std::int32_t> parent_cell;
-  parent_cell.reserve(refined_cell_count);
+
+  std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
+  if (compute_parent_cell)
+  {
+    parent_cell.emplace();
+    parent_cell->reserve(refined_cell_count);
+  }
+
   for (std::int32_t cell = 0; cell < map_c->size_local(); ++cell)
   {
     const auto& vertices = c_to_v->links(cell);
@@ -145,71 +161,30 @@ compute_interval_refinement(const mesh::Mesh<T>& mesh,
 
       // Add new cells/edges to refined topology
       cell_topology.insert(cell_topology.end(), {a, c, c, b});
-      parent_cell.insert(parent_cell.end(), {cell, cell});
+
+      if (compute_parent_cell)
+        parent_cell->insert(parent_cell->end(), {cell, cell});
     }
     else
     {
       // Copy the previous cell
       cell_topology.insert(cell_topology.end(), {a, b});
-      parent_cell.push_back(cell);
+
+      if (compute_parent_cell)
+        parent_cell->push_back(cell);
     }
   }
 
   assert(cell_topology.size() == 2 * refined_cell_count);
-  assert(parent_cell.size() == refined_cell_count);
+  assert(parent_cell->size() == (compute_parent_cell ? refined_cell_count : 0));
 
   std::vector<std::int32_t> offsets(refined_cell_count + 1);
   std::ranges::generate(offsets, [i = 0]() mutable { return 2 * i++; });
 
-  return {graph::AdjacencyList(std::move(cell_topology), std::move(offsets)),
-          std::move(new_vertex_coords), xshape, std::move(parent_cell)};
-}
-
-} // namespace impl
+  graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
 
-/// @brief Refine a (topologically) one-dimensional mesh by splitting
-/// cells, i.e. edges.
-///
-/// @param[in] mesh Mesh to be refined.
-/// @param[in] cells Optional indices of the cells that should be
-/// refined by this refinement. If not provided, all cells are
-/// considered marked for refinement, i.e. a uniform refinement is
-/// performed.
-/// @param[in] redistribute Option to enable redistribution of the
-/// refined mesh across processes.
-/// @param[in] ghost_mode Ghost mode of the refined mesh, default is
-/// ghost mode none.
-///
-/// @return Refined mesh, and list of parent cells and an array mapping
-/// the child cell index of the refined mesh to its parent cell index in
-/// the unrefined mesh.
-template <std::floating_point T>
-std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>>
-refine_interval(const mesh::Mesh<T>& mesh,
-                std::optional<std::span<const std::int32_t>> cells,
-                bool redistribute,
-                mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet)
-{
-  if (mesh.topology()->cell_type() != mesh::CellType::interval)
-    throw std::runtime_error("Cell type not supported");
-  assert(mesh.topology()->dim() == 1);
-  assert(mesh.topology()->index_map(1));
-
-  auto [cell_adj, new_coords, xshape, parent_cell]
-      = impl::compute_interval_refinement(mesh, cells);
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                              mesh.geometry().cmap(), new_coords, xshape,
-                              mesh::GhostMode::none),
-            std::move(parent_cell)};
-  }
-  else
-  {
-    return {partition<T>(mesh, cell_adj, std::span(new_coords), xshape,
-                         redistribute, ghost_mode),
-            std::move(parent_cell)};
-  }
+  return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
+          std::move(parent_cell), std::nullopt};
 }
 
-} // namespace dolfinx::refinement
+} // namespace dolfinx::refinement::interval
diff --git a/cpp/dolfinx/refinement/option.h b/cpp/dolfinx/refinement/option.h
new file mode 100644
index 00000000000..a1c5acdbb1c
--- /dev/null
+++ b/cpp/dolfinx/refinement/option.h
@@ -0,0 +1,52 @@
+// Copyright (C) 2024 Paul Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#pragma once
+
+#include <cstdint>
+#include <type_traits>
+
+namespace dolfinx::refinement
+{
+/// @brief Options for data to compute during mesh refinement.
+enum class Option : std::uint8_t
+{
+  none = 0b00,         /*!< No extra data */
+  parent_facet = 0b01, /*!< Compute list of the cell-local facet indices in the
+          parent cell of each facet in each new cell (or -1 if no match)  */
+  parent_cell
+  = 0b10, /*!< Compute list with the parent cell index for each new cell */
+  parent_cell_and_facet = 0b11 /*< Both cell and facet parent data */
+};
+
+/// @brief Combine two refinement options into one, both flags will be
+/// set for the resulting option.
+inline constexpr Option operator|(Option a, Option b)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t b_native = static_cast<bitmask_t>(b);
+  return static_cast<Option>(a_native | b_native);
+}
+
+/// @brief Check if parent_facet flag is set
+inline constexpr bool option_parent_facet(Option a)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t facet_native = static_cast<bitmask_t>(Option::parent_facet);
+  return (a_native & facet_native) == facet_native;
+}
+
+/// @brief Check if parent_cell flag is set
+inline constexpr bool option_parent_cell(Option a)
+{
+  using bitmask_t = std::underlying_type_t<Option>;
+  bitmask_t a_native = static_cast<bitmask_t>(a);
+  bitmask_t facet_native = static_cast<bitmask_t>(Option::parent_cell);
+  return (a_native & facet_native) == facet_native;
+}
+} // namespace dolfinx::refinement
diff --git a/cpp/dolfinx/refinement/plaza.h b/cpp/dolfinx/refinement/plaza.h
index 9231144a84a..65cdec06645 100644
--- a/cpp/dolfinx/refinement/plaza.h
+++ b/cpp/dolfinx/refinement/plaza.h
@@ -4,16 +4,16 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
+#include "dolfinx/graph/AdjacencyList.h"
+#include "dolfinx/mesh/Mesh.h"
+#include "dolfinx/mesh/Topology.h"
+#include "dolfinx/mesh/utils.h"
+#include "option.h"
 #include "utils.h"
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
-#include <dolfinx/common/MPI.h>
-#include <dolfinx/graph/AdjacencyList.h>
-#include <dolfinx/mesh/Geometry.h>
-#include <dolfinx/mesh/Mesh.h>
-#include <dolfinx/mesh/Topology.h>
-#include <dolfinx/mesh/utils.h>
+#include <optional>
 #include <span>
 #include <tuple>
 #include <utility>
@@ -28,19 +28,6 @@
 /// Applied Numerical Mathematics 32 (2000), 195-218.
 namespace dolfinx::refinement::plaza
 {
-
-/// @brief Options for data to compute during mesh refinement.
-enum class Option : int
-{
-  none = 0, /*!< No extra data */
-  parent_cell
-  = 1, /*!< Compute list with the parent cell index for each new cell  */
-  parent_facet
-  = 2, /*!< Compute list of the cell-local facet indices in the parent cell of
-          each facet in each new cell (or -1 if no match) */
-  parent_cell_and_facet = 3 /*!< Both cell and facet parent data */
-};
-
 namespace impl
 {
 /// Computes the parent-child facet relationship
@@ -116,21 +103,23 @@ auto compute_parent_facets(std::span<const std::int32_t> simplex_set)
   return parent_facet;
 }
 
-/// Get the subdivision of an original simplex into smaller simplices,
+/// @brief Get the subdivision of an original simplex into smaller simplices,
 /// for a given set of marked edges, and the longest edge of each facet
-/// (cell local indexing). A flag indicates if a uniform subdivision is
-/// preferable in 2D.
+/// (cell local indexing).
 ///
-/// @param[in] indices Vector containing the global indices for the original
-/// vertices and potential new vertices at each edge. Size (num_vertices +
-/// num_edges). If an edge is not refined its corresponding entry is -1
+/// A flag indicates if a uniform subdivision is preferable in 2D.
+///
+/// @param[in] indices Vector containing the global indices for the
+/// original vertices and potential new vertices at each edge. Size
+/// (num_vertices + num_edges). If an edge is not refined its
+/// corresponding entry is -1.
 /// @param[in] longest_edge Vector indicating the longest edge for each
-/// triangle in the cell. For triangular cells (2D) there is only one value,
-/// and for tetrahedra (3D) there are four values, one for each facet. The
-/// values give the local edge indices of the cell.
-/// @param[in] tdim Topological dimension (2 or 3)
-/// @param[in] uniform Make a "uniform" subdivision with all triangles being
-/// similar shape
+/// triangle in the cell. For triangular cells (2D) there is only one
+/// value, and for tetrahedra (3D) there are four values, one for each
+/// facet. The values give the local edge indices of the cell.
+/// @param[in] tdim Topological dimension (2 or 3).
+/// @param[in] uniform Make a "uniform" subdivision with all triangles.
+/// being similar shape.
 /// @return
 std::pair<std::array<std::int32_t, 32>, std::size_t>
 get_simplices(std::span<const std::int64_t> indices,
@@ -144,15 +133,15 @@ void enforce_rules(MPI_Comm comm, const graph::AdjacencyList<int>& shared_edges,
                    const mesh::Topology& topology,
                    std::span<const std::int32_t> long_edge);
 
-/// @brief Get the longest edge of each face (using local mesh index)
+/// @brief Get the longest edge of each face (using local mesh index).
 ///
-/// @note Edge ratio ok returns an array in 2D, where it checks if the ratio
-/// between the shortest and longest edge of a cell is bigger than sqrt(2)/2. In
-/// 3D an empty array is returned
+/// @note Edge ratio ok returns an array in 2D, where it checks if the
+/// ratio between the shortest and longest edge of a cell is bigger than
+/// sqrt(2)/2. In 3D an empty array is returned.
 ///
 /// @param[in] mesh The mesh
-/// @return A tuple (longest edge, edge ratio ok) where longest edge gives the
-/// local index of the longest edge for each face.
+/// @return A tuple (longest edge, edge ratio ok) where longest edge
+/// gives the local index of the longest edge for each face.
 template <std::floating_point T>
 std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
 face_long_edge(const mesh::Mesh<T>& mesh)
@@ -295,30 +284,34 @@ face_long_edge(const mesh::Mesh<T>& mesh)
 /// and (5) map from refined facets to parent facets.
 template <std::floating_point T>
 std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
-           std::array<std::size_t, 2>, std::vector<std::int32_t>,
-           std::vector<std::int8_t>>
+           std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
+           std::optional<std::vector<std::int8_t>>>
 compute_refinement(MPI_Comm neighbor_comm,
                    std::span<const std::int8_t> marked_edges,
                    const graph::AdjacencyList<int>& shared_edges,
                    const mesh::Mesh<T>& mesh,
                    std::span<const std::int32_t> long_edge,
-                   std::span<const std::int8_t> edge_ratio_ok,
-                   plaza::Option option)
+                   std::span<const std::int8_t> edge_ratio_ok, Option option)
 {
   int tdim = mesh.topology()->dim();
   int num_cell_edges = tdim * 3 - 3;
   int num_cell_vertices = tdim + 1;
-  bool compute_facets = option == plaza::Option::parent_facet
-                        or option == plaza::Option::parent_cell_and_facet;
-  bool compute_parent_cell = option == plaza::Option::parent_cell
-                             or option == plaza::Option::parent_cell_and_facet;
+
+  bool compute_facets = option_parent_facet(option);
+  bool compute_parent_cell = option_parent_cell(option);
 
   // Make new vertices in parallel
   const auto [new_vertex_map, new_vertex_coords, xshape]
       = create_new_vertices(neighbor_comm, shared_edges, mesh, marked_edges);
 
-  std::vector<std::int32_t> parent_cell;
-  std::vector<std::int8_t> parent_facet;
+  std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
+  if (compute_parent_cell)
+    parent_cell.emplace();
+
+  std::optional<std::vector<std::int8_t>> parent_facet(std::nullopt);
+  if (compute_facets)
+    parent_facet.emplace();
+
   std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
   std::vector<std::int32_t> simplex_set;
 
@@ -375,14 +368,14 @@ compute_refinement(MPI_Comm neighbor_comm,
         cell_topology.push_back(global_indices[v]);
 
       if (compute_parent_cell)
-        parent_cell.push_back(c);
+        parent_cell->push_back(c);
 
       if (compute_facets)
       {
         if (tdim == 3)
-          parent_facet.insert(parent_facet.end(), {0, 1, 2, 3});
+          parent_facet->insert(parent_facet->end(), {0, 1, 2, 3});
         else
-          parent_facet.insert(parent_facet.end(), {0, 1, 2});
+          parent_facet->insert(parent_facet->end(), {0, 1, 2});
       }
     }
     else
@@ -417,7 +410,7 @@ compute_refinement(MPI_Comm neighbor_comm,
       if (compute_parent_cell)
       {
         for (std::int32_t i = 0; i < ncells; ++i)
-          parent_cell.push_back(c);
+          parent_cell->push_back(c);
       }
 
       if (compute_facets)
@@ -425,14 +418,14 @@ compute_refinement(MPI_Comm neighbor_comm,
         if (tdim == 3)
         {
           auto npf = compute_parent_facets<3>(simplex_set);
-          parent_facet.insert(parent_facet.end(), npf.begin(),
-                              std::next(npf.begin(), simplex_set.size()));
+          parent_facet->insert(parent_facet->end(), npf.begin(),
+                               std::next(npf.begin(), simplex_set.size()));
         }
         else
         {
           auto npf = compute_parent_facets<2>(simplex_set);
-          parent_facet.insert(parent_facet.end(), npf.begin(),
-                              std::next(npf.begin(), simplex_set.size()));
+          parent_facet->insert(parent_facet->end(), npf.begin(),
+                               std::next(npf.begin(), simplex_set.size()));
         }
       }
 
@@ -454,164 +447,6 @@ compute_refinement(MPI_Comm neighbor_comm,
 }
 } // namespace impl
 
-/// @brief Uniform refine, optionally redistributing and optionally
-/// calculating the parent-child relationships`.
-///
-/// @param[in] mesh Input mesh to be refined
-/// @param[in] redistribute Flag to call the mesh partitioner to
-/// redistribute after refinement
-/// @param[in] option Control the computation of parent facets, parent
-/// cells. If an option is unselected, an empty list is returned.
-/// @return Refined mesh and optional parent cell index, parent facet
-/// indices
-template <std::floating_point T>
-std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
-refine(const mesh::Mesh<T>& mesh, bool redistribute, Option option)
-{
-  auto [cell_adj, new_coords, xshape, parent_cell, parent_facet]
-      = compute_refinement_data(mesh, option);
-
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                              mesh.geometry().cmap(), new_coords, xshape,
-                              mesh::GhostMode::none),
-            std::move(parent_cell), std::move(parent_facet)};
-  }
-  else
-  {
-    std::shared_ptr<const common::IndexMap> map_c
-        = mesh.topology()->index_map(mesh.topology()->dim());
-    const int num_ghost_cells = map_c->num_ghosts();
-    // Check if mesh has ghost cells on any rank
-    // FIXME: this is not a robust test. Should be user option.
-    int max_ghost_cells = 0;
-    MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
-                  mesh.comm());
-
-    // Build mesh
-    const mesh::GhostMode ghost_mode = max_ghost_cells == 0
-                                           ? mesh::GhostMode::none
-                                           : mesh::GhostMode::shared_facet;
-    return {partition<T>(mesh, cell_adj, std::span(new_coords), xshape,
-                         redistribute, ghost_mode),
-            std::move(parent_cell), std::move(parent_facet)};
-  }
-}
-
-/// @brief Refine with markers, optionally redistributing, and
-/// optionally calculating the parent-child relationships.
-///
-/// @param[in] mesh Input mesh to be refined
-/// @param[in] edges Indices of the edges that should be split by this
-/// refinement
-/// @param[in] redistribute Flag to call the Mesh Partitioner to
-/// redistribute after refinement
-/// @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 indices
-template <std::floating_point T>
-std::tuple<mesh::Mesh<T>, std::vector<std::int32_t>, std::vector<std::int8_t>>
-refine(const mesh::Mesh<T>& mesh, std::span<const std::int32_t> edges,
-       bool redistribute, Option option)
-{
-  auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
-      = compute_refinement_data(mesh, edges, option);
-
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    return {mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                              mesh.geometry().cmap(), new_vertex_coords, xshape,
-                              mesh::GhostMode::none),
-            std::move(parent_cell), std::move(parent_facet)};
-  }
-  else
-  {
-    std::shared_ptr<const common::IndexMap> map_c
-        = mesh.topology()->index_map(mesh.topology()->dim());
-    const int num_ghost_cells = map_c->num_ghosts();
-    // Check if mesh has ghost cells on any rank
-    // FIXME: this is not a robust test. Should be user option.
-    int max_ghost_cells = 0;
-    MPI_Allreduce(&num_ghost_cells, &max_ghost_cells, 1, MPI_INT, MPI_MAX,
-                  mesh.comm());
-
-    // Build mesh
-    const mesh::GhostMode ghost_mode = max_ghost_cells == 0
-                                           ? mesh::GhostMode::none
-                                           : mesh::GhostMode::shared_facet;
-
-    return {partition<T>(mesh, cell_adj, new_vertex_coords, xshape,
-                         redistribute, ghost_mode),
-            std::move(parent_cell), std::move(parent_facet)};
-  }
-}
-
-/// @brief Refine mesh returning new mesh data.
-///
-/// @param[in] mesh Input mesh to be refined
-/// @param[in] option Control computation of parent facets and parent
-/// cells. If an option is unselected, an empty list is returned.
-/// @return New mesh data: cell topology, vertex coordinates, vertex
-/// coordinates shape, and optional parent cell index, and parent facet
-/// indices.
-template <std::floating_point T>
-std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
-           std::array<std::size_t, 2>, std::vector<std::int32_t>,
-           std::vector<std::int8_t>>
-compute_refinement_data(const mesh::Mesh<T>& mesh, Option option)
-{
-  common::Timer t0("PLAZA: refine");
-  auto topology = mesh.topology();
-  assert(topology);
-
-  if (topology->cell_type() != mesh::CellType::triangle
-      and topology->cell_type() != mesh::CellType::tetrahedron)
-  {
-    throw std::runtime_error("Cell type not supported");
-  }
-
-  auto map_e = topology->index_map(1);
-  if (!map_e)
-    throw std::runtime_error("Edges must be initialised");
-
-  // Get sharing ranks for each edge
-  graph::AdjacencyList<int> edge_ranks = map_e->index_to_dest_ranks();
-
-  // Create unique list of ranks that share edges (owners of ghosts
-  // plus ranks that ghost owned indices)
-  std::vector<int> ranks(edge_ranks.array().begin(), edge_ranks.array().end());
-  std::ranges::sort(ranks);
-  auto [unique_end, range_end] = std::ranges::unique(ranks);
-  ranks.erase(unique_end, range_end);
-
-  // Convert edge_ranks from global rank to to neighbourhood ranks
-  std::ranges::transform(edge_ranks.array(), edge_ranks.array().begin(),
-                         [&ranks](auto r)
-                         {
-                           auto it = std::ranges::lower_bound(ranks, r);
-                           assert(it != ranks.end() and *it == r);
-                           return std::distance(ranks.begin(), it);
-                         });
-
-  MPI_Comm comm;
-  MPI_Dist_graph_create_adjacent(mesh.comm(), ranks.size(), ranks.data(),
-                                 MPI_UNWEIGHTED, ranks.size(), ranks.data(),
-                                 MPI_UNWEIGHTED, MPI_INFO_NULL, false, &comm);
-
-  const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(mesh);
-  auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
-      = impl::compute_refinement(
-          comm,
-          std::vector<std::int8_t>(map_e->size_local() + map_e->num_ghosts(),
-                                   true),
-          edge_ranks, mesh, long_edge, edge_ratio_ok, option);
-  MPI_Comm_free(&comm);
-
-  return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
-          std::move(parent_cell), std::move(parent_facet)};
-}
-
 /// Refine with markers returning new mesh data.
 ///
 /// @param[in] mesh Input mesh to be refined
@@ -623,10 +458,11 @@ compute_refinement_data(const mesh::Mesh<T>& mesh, Option option)
 /// cell index, and stored parent facet indices (if requested).
 template <std::floating_point T>
 std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
-           std::array<std::size_t, 2>, std::vector<std::int32_t>,
-           std::vector<std::int8_t>>
+           std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
+           std::optional<std::vector<std::int8_t>>>
 compute_refinement_data(const mesh::Mesh<T>& mesh,
-                        std::span<const std::int32_t> edges, Option option)
+                        std::optional<std::span<const std::int32_t>> edges,
+                        Option option)
 {
   common::Timer t0("PLAZA: refine");
   auto topology = mesh.topology();
@@ -663,17 +499,21 @@ compute_refinement_data(const mesh::Mesh<T>& mesh,
 
   // Get number of neighbors
   std::vector<std::int8_t> marked_edges(
-      map_e->size_local() + map_e->num_ghosts(), false);
+      map_e->size_local() + map_e->num_ghosts(), !edges.has_value());
   std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
-  for (auto edge : edges)
+
+  if (edges.has_value())
   {
-    if (!marked_edges[edge])
+    for (auto edge : edges.value())
     {
-      marked_edges[edge] = true;
+      if (!marked_edges[edge])
+      {
+        marked_edges[edge] = true;
 
-      // If it is a shared edge, add all sharing neighbors to update set
-      for (int rank : edge_ranks.links(edge))
-        marked_for_update[rank].push_back(edge);
+        // If it is a shared edge, add all sharing neighbors to update set
+        for (int rank : edge_ranks.links(edge))
+          marked_for_update[rank].push_back(edge);
+      }
     }
   }
 
@@ -698,4 +538,5 @@ compute_refinement_data(const mesh::Mesh<T>& mesh,
   return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
           std::move(parent_cell), std::move(parent_facet)};
 }
+
 } // namespace dolfinx::refinement::plaza
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 5b243fc4f29..32f18d89827 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -6,74 +6,85 @@
 
 #pragma once
 
+#include "dolfinx/mesh/Mesh.h"
+#include "dolfinx/mesh/Topology.h"
+#include "dolfinx/mesh/cell_types.h"
+#include "dolfinx/mesh/utils.h"
+#include "interval.h"
 #include "plaza.h"
-#include "refine.h"
-#include <dolfinx/common/IndexMap.h>
-#include <dolfinx/common/log.h>
-#include <dolfinx/mesh/Mesh.h>
-#include <dolfinx/mesh/Topology.h>
+#include <algorithm>
+#include <concepts>
+#include <optional>
+#include <utility>
 
 namespace dolfinx::refinement
 {
-/// @brief Create a uniformly refined mesh.
-///
-/// @param[in] mesh Mesh to create a new, refined mesh from
-/// @param[in] redistribute If `true` refined mesh is re-partitioned
-/// across MPI ranks.
-/// @return Refined mesh
+namespace impl
+{
+/// @brief create the refined mesh by optionally redistributing it
 template <std::floating_point T>
-mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh, bool redistribute = true)
+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)
 {
-  auto topology = mesh.topology();
-  assert(topology);
-
-  if (topology->cell_type() != mesh::CellType::triangle
-      and topology->cell_type() != mesh::CellType::tetrahedron)
+  if (dolfinx::MPI::size(mesh.comm()) == 1)
   {
-    throw std::runtime_error("Refinement only defined for simplices");
+    // 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);
   }
-
-  auto [refined_mesh, parent_cell, parent_facet]
-      = plaza::refine(mesh, redistribute, plaza::Option::none);
-
-  // 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();
-  spdlog::info(
-      "Number of cells increased from {} to {} ({}% increase).", n0, n1,
-      100.0 * (static_cast<double>(n1) / static_cast<double>(n0) - 1.0));
-
-  return refined_mesh;
 }
+} // namespace impl
 
-/// @brief Create a locally refined mesh.
+/// @brief Refine with markers, optionally redistributing, and
+/// optionally calculating the parent-child relationships.
 ///
-/// @param[in] mesh Mesh to create a new, refined mesh from.
-/// @param[in] edges Indices of the edges that should be split during
-/// refinement. mesh::compute_incident_entities can be used to compute
-/// the edges that are incident to other entities, e.g. incident to
-/// cells.
-/// @param[in] redistribute If `true` refined mesh is re-partitioned
-/// across MPI ranks.
-/// @return Refined mesh.
+/// @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.
+/// @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
+/// indices.
 template <std::floating_point T>
-mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh,
-                     std::span<const std::int32_t> edges,
-                     bool redistribute = true)
+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,
+       Option option = Option::none)
 {
   auto topology = mesh.topology();
   assert(topology);
-  if (topology->cell_type() != mesh::CellType::triangle
-      and topology->cell_type() != mesh::CellType::tetrahedron)
-  {
+
+  mesh::CellType cell_t = topology->cell_type();
+  if (!mesh::is_simplex(cell_t))
     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);
 
-  auto [refined_mesh, parent_cell, parent_facet]
-      = plaza::refine(mesh, edges, redistribute, plaza::Option::none);
+  mesh::Mesh<T> refined_mesh = impl::create_refined_mesh(
+      mesh, cell_adj, std::span<const T>(new_vertex_coords), xshape,
+      redistribute, ghost_mode);
 
-  // Report the number of refined cells
+  // Report the number of refined cellse
   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();
@@ -81,7 +92,8 @@ mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh,
       "Number of cells increased from {} to {} ({}% increase).", n0, n1,
       100.0 * (static_cast<double>(n1) / static_cast<double>(n0) - 1.0));
 
-  return refined_mesh;
+  return {std::move(refined_mesh), 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 3d4cb017ae2..275b40d7cc2 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -35,23 +35,21 @@ class IndexMap;
 
 namespace dolfinx::refinement
 {
-
 namespace impl
 {
-
 /// @brief  Compute global index
 std::int64_t local_to_global(std::int32_t local_index,
                              const common::IndexMap& map);
 
-/// Create geometric points of new Mesh, from current Mesh and a
+/// @brief Create geometric points of new Mesh, from current Mesh and a
 /// edge_to_vertex map listing the new local points (midpoints of those
-/// edges)
+/// edges).
 ///
-/// @param mesh Current mesh
+/// @param mesh Current mesh.
 /// @param local_edge_to_new_vertex A map from a local edge to the new
-/// global vertex index that will be inserted at its midpoint
+/// global vertex index that will be inserted at its midpoint.
 /// @return (1) Array of new (flattened) mesh geometry and (2) its
-/// multi-dimensional shape
+/// multi-dimensional shape.
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 2>> create_new_geometry(
     const mesh::Mesh<T>& mesh,
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index a12aac898e4..95ccb34e614 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -45,6 +45,7 @@ add_executable(
   mesh/generation.cpp
   common/CIFailure.cpp
   refinement/interval.cpp
+  refinement/option.cpp
   ${CMAKE_CURRENT_BINARY_DIR}/poisson.c
 )
 target_link_libraries(unittests PRIVATE Catch2::Catch2WithMain dolfinx)
diff --git a/cpp/test/refinement/interval.cpp b/cpp/test/refinement/interval.cpp
index 7e3b90e5262..90ab09c4a7e 100644
--- a/cpp/test/refinement/interval.cpp
+++ b/cpp/test/refinement/interval.cpp
@@ -54,8 +54,10 @@ TEST_CASE("Interval uniform refinement", "refinement,interval,uniform")
   mesh::Mesh<double> mesh = create_3_vertex_interval_mesh();
   mesh.topology()->create_connectivity(1, 0);
 
-  auto [refined_mesh, parent_edge]
-      = refinement::refine_interval(mesh, std::nullopt, false);
+  // TODO: parent_facet
+  auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
+      mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
+      refinement::Option::parent_cell);
 
   // Check geometry
   {
@@ -113,12 +115,14 @@ TEST_CASE("Interval uniform refinement", "refinement,interval,uniform")
 
   // Check parent edges
   {
-    CHECK(parent_edge.size() == 4);
+    CHECK(parent_edge.has_value());
 
-    CHECK(parent_edge[0] == 0);
-    CHECK(parent_edge[1] == 0);
-    CHECK(parent_edge[2] == 1);
-    CHECK(parent_edge[3] == 1);
+    CHECK(parent_edge.value().size() == 4);
+
+    CHECK(parent_edge.value()[0] == 0);
+    CHECK(parent_edge.value()[1] == 0);
+    CHECK(parent_edge.value()[2] == 1);
+    CHECK(parent_edge.value()[3] == 1);
   }
 }
 
@@ -131,8 +135,10 @@ TEST_CASE("Interval adaptive refinement", "refinement,interval,adaptive")
   mesh.topology()->create_connectivity(1, 0);
 
   std::vector<std::int32_t> edges{1};
-  auto [refined_mesh, parent_edge]
-      = refinement::refine_interval(mesh, std::span(edges), false);
+  // TODO: parent_facet
+  auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
+      mesh, std::span(edges), false, mesh::GhostMode::shared_facet,
+      refinement::Option::parent_cell);
 
   // Check geometry
   {
@@ -181,11 +187,13 @@ TEST_CASE("Interval adaptive refinement", "refinement,interval,adaptive")
 
   // Check parent edges
   {
-    CHECK(parent_edge.size() == 3);
+    CHECK(parent_edge.has_value());
+
+    CHECK(parent_edge.value().size() == 3);
 
-    CHECK(parent_edge[0] == 0);
-    CHECK(parent_edge[1] == 1);
-    CHECK(parent_edge[2] == 1);
+    CHECK(parent_edge.value()[0] == 0);
+    CHECK(parent_edge.value()[1] == 1);
+    CHECK(parent_edge.value()[2] == 1);
   }
 }
 
@@ -239,8 +247,10 @@ TEST_CASE("Interval Refinement (parallel)", "refinement,interval,paralle")
 
   // complete refinement
   {
-    auto [refined_mesh, parent_edges]
-        = refinement::refine_interval(mesh, std::nullopt, false);
+    // TODO: parent_facet
+    auto [refined_mesh, parent_edges, parent_facet] = refinement::refine(
+        mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
+        refinement::Option::parent_cell);
 
     // Check geometry
     {
diff --git a/cpp/test/refinement/option.cpp b/cpp/test/refinement/option.cpp
new file mode 100644
index 00000000000..54219f83759
--- /dev/null
+++ b/cpp/test/refinement/option.cpp
@@ -0,0 +1,49 @@
+// Copyright (C) 2024 Paul Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include <catch2/catch_test_macros.hpp>
+#include <cstdint>
+#include <dolfinx/refinement/option.h>
+
+using namespace dolfinx::refinement;
+
+TEST_CASE("Refinement Option", "refinement,option")
+{
+  // Internal binary structure
+  CHECK(static_cast<std::uint8_t>(Option::none) == 0);
+  CHECK(static_cast<std::uint8_t>(Option::parent_facet) == 1);
+  CHECK(static_cast<std::uint8_t>(Option::parent_cell) == 2);
+  CHECK(static_cast<std::uint8_t>(Option::parent_cell_and_facet) == 3);
+
+  // Extraction of flags from possible multiple set flags
+  CHECK(option_parent_cell(Option::none) == false);
+  CHECK(option_parent_facet(Option::none) == false);
+
+  CHECK(option_parent_cell(Option::parent_facet) == false);
+  CHECK(option_parent_facet(Option::parent_facet) == true);
+
+  CHECK(option_parent_cell(Option::parent_cell) == true);
+  CHECK(option_parent_facet(Option::parent_cell) == false);
+
+  CHECK(option_parent_cell(Option::parent_cell_and_facet) == true);
+  CHECK(option_parent_facet(Option::parent_cell_and_facet) == true);
+
+  // Logical combination of options
+  CHECK((Option::none | Option::none) == Option::none);
+  CHECK((Option::none | Option::parent_facet) == Option::parent_facet);
+  CHECK((Option::none | Option::parent_cell) == Option::parent_cell);
+
+  CHECK((Option::parent_cell_and_facet | Option::none)
+        == Option::parent_cell_and_facet);
+  CHECK((Option::parent_cell_and_facet | Option::parent_facet)
+        == Option::parent_cell_and_facet);
+  CHECK((Option::parent_cell_and_facet | Option::parent_cell)
+        == Option::parent_cell_and_facet);
+
+  CHECK((Option::parent_cell | Option::none) == Option::parent_cell);
+  CHECK((Option::parent_cell | Option::parent_facet)
+        == Option::parent_cell_and_facet);
+}
\ No newline at end of file
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index cf04e746fe6..3d58b531bdd 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -60,8 +60,6 @@
     "create_unit_cube",
     "to_type",
     "to_string",
-    "refine_interval",
-    "refine_plaza",
     "transfer_meshtag",
     "entities_to_geometry",
 ]
@@ -313,65 +311,12 @@ def transfer_meshtag(
 
 
 def refine(
-    mesh: Mesh, edges: typing.Optional[np.ndarray] = None, redistribute: bool = True
-) -> Mesh:
-    """Refine a mesh.
-
-    Args:
-        mesh: 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``.
-
-    Returns:
-       Refined mesh.
-    """
-    if edges is None:
-        mesh1 = _cpp.refinement.refine(mesh._cpp_object, redistribute)
-    else:
-        mesh1 = _cpp.refinement.refine(mesh._cpp_object, edges, redistribute)
-    ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element())  # type: ignore
-    return Mesh(mesh1, ufl_domain)
-
-
-def refine_interval(
-    mesh: Mesh,
-    cells: typing.Optional[np.ndarray] = None,
-    redistribute: bool = True,
-    ghost_mode: GhostMode = GhostMode.shared_facet,
-) -> tuple[Mesh, npt.NDArray[np.int32]]:
-    """Refine a (topologically) one dimensional mesh.
-
-    Args:
-        mesh: Mesh to refine
-        cells: Indices of cells, i.e. edges, to split druing refinement. If ``None``, mesh
-            refinement is uniform.
-        redistribute: Refined mesh is re-partitioned if ``True``.
-        ghost_mode: ghost mode of the refined mesh
-
-    Returns:
-        Refined mesh and parent cells
-    """
-
-    if cells is None:
-        refined_mesh, parent_cells = _cpp.refinement.refine_interval(
-            mesh._cpp_object, redistribute, ghost_mode
-        )
-    else:
-        refined_mesh, parent_cells = _cpp.refinement.refine_interval(
-            mesh._cpp_object, cells, redistribute, ghost_mode
-        )
-
-    return Mesh(refined_mesh, mesh._ufl_domain), parent_cells
-
-
-def refine_plaza(
     mesh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
     redistribute: bool = True,
+    ghost_mode: GhostMode = GhostMode.shared_facet,
     option: RefinementOption = RefinementOption.none,
-) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int32]]:
+) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
 
     Args:
@@ -380,19 +325,14 @@ def refine_plaza(
             mesh refinement is uniform.
         redistribute:
             Refined mesh is re-partitioned if ``True``.
-        option:
-            Control computation of the parent-refined mesh data.
 
     Returns:
-       Refined mesh, list of parent cell for each refine cell, and list
+       Refined mesh.
     """
-    if edges is None:
-        mesh1, cells, facets = _cpp.refinement.refine_plaza(mesh._cpp_object, redistribute, option)
-    else:
-        mesh1, cells, facets = _cpp.refinement.refine_plaza(
-            mesh._cpp_object, edges, redistribute, option
-        )
-    return Mesh(mesh1, mesh._ufl_domain), cells, facets
+    mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
+        mesh._cpp_object, edges, redistribute, ghost_mode, option
+    )
+    return Mesh(mesh1, mesh._ufl_domain), parent_cell, parent_facet
 
 
 def create_mesh(
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index af5337a5684..37ca3937992 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -6,18 +6,22 @@
 
 #include "array.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 <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
+#include <nanobind/stl/optional.h>
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/tuple.h>
 #include <nanobind/stl/vector.h>
 #include <optional>
+#include <span>
 
 namespace nb = nanobind;
 
@@ -27,76 +31,40 @@ namespace dolfinx_wrappers
 template <std::floating_point T>
 void export_refinement_with_variable_mesh_type(nb::module_& m)
 {
-  m.def("refine",
-        nb::overload_cast<const dolfinx::mesh::Mesh<T>&, bool>(
-            &dolfinx::refinement::refine<T>),
-        nb::arg("mesh"), nb::arg("redistribute") = true);
   m.def(
       "refine",
       [](const dolfinx::mesh::Mesh<T>& mesh,
-         nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> edges,
-         bool redistribute)
+         std::optional<
+             nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
+             edges,
+         bool redistribute, dolfinx::mesh::GhostMode ghost_mode,
+         dolfinx::refinement::Option option)
       {
-        return dolfinx::refinement::refine(
-            mesh, std::span(edges.data(), edges.size()), redistribute);
-      },
-      nb::arg("mesh"), nb::arg("edges"), nb::arg("redistribute") = true);
+        std::optional<std::span<const std::int32_t>> cpp_edges(std::nullopt);
+        if (edges.has_value())
+        {
+          cpp_edges.emplace(
+              std::span(edges.value().data(), edges.value().size()));
+        }
 
-  m.def(
-      "refine_interval",
-      [](const dolfinx::mesh::Mesh<T>& mesh, bool redistribute,
-         dolfinx::mesh::GhostMode ghost_mode)
-      {
-        auto [mesh_refined, parent_cells]
-            = dolfinx::refinement::refine_interval(mesh, std::nullopt,
-                                                   redistribute, ghost_mode);
-        return std::tuple(std::move(mesh_refined),
-                          as_nbarray(std::move(parent_cells)));
-      },
-      nb::arg("mesh"), nb::arg("redistribute"), nb::arg("ghost_mode"));
+        auto [mesh1, cell, facet] = dolfinx::refinement::refine(
+            mesh, cpp_edges, redistribute, ghost_mode, option);
 
-  m.def(
-      "refine_interval",
-      [](const dolfinx::mesh::Mesh<T>& mesh,
-         nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> cells,
-         bool redistribute, dolfinx::mesh::GhostMode ghost_mode)
-      {
-        auto [mesh_refined, parent_cells]
-            = dolfinx::refinement::refine_interval(
-                mesh, std::make_optional(std::span(cells.data(), cells.size())),
-                redistribute, ghost_mode);
-        return std::tuple(std::move(mesh_refined),
-                          as_nbarray(std::move(parent_cells)));
-      },
-      nb::arg("mesh"), nb::arg("edges"), nb::arg("redistribute"),
-      nb::arg("ghost_mode"));
+        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())));
 
-  m.def(
-      "refine_plaza",
-      [](const dolfinx::mesh::Mesh<T>& mesh0, bool redistribute,
-         dolfinx::refinement::plaza::Option option)
-      {
-        auto [mesh1, cell, facet]
-            = dolfinx::refinement::plaza::refine(mesh0, redistribute, option);
-        return std::tuple{std::move(mesh1), as_nbarray(std::move(cell)),
-                          as_nbarray(std::move(facet))};
-      },
-      nb::arg("mesh"), nb::arg("redistribute"), nb::arg("option"));
+        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())));
 
-  m.def(
-      "refine_plaza",
-      [](const dolfinx::mesh::Mesh<T>& mesh0,
-         nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> edges,
-         bool redistribute, dolfinx::refinement::plaza::Option option)
-      {
-        auto [mesh1, cell, facet] = dolfinx::refinement::plaza::refine(
-            mesh0, std::span<const std::int32_t>(edges.data(), edges.size()),
-            redistribute, option);
-        return std::tuple{std::move(mesh1), as_nbarray(std::move(cell)),
-                          as_nbarray(std::move(facet))};
+        return std::tuple{std::move(mesh1), std::move(python_cell),
+                          std::move(python_facet)};
       },
-      nb::arg("mesh"), nb::arg("edges"), nb::arg("redistribute"),
-      nb::arg("option"));
+      nb::arg("mesh"), nb::arg("edges") = nb::none(), nb::arg("redistribute"),
+      nb::arg("ghost_mode"), nb::arg("option"));
 }
 
 void refinement(nb::module_& m)
@@ -104,12 +72,12 @@ void refinement(nb::module_& m)
   export_refinement_with_variable_mesh_type<float>(m);
   export_refinement_with_variable_mesh_type<double>(m);
 
-  nb::enum_<dolfinx::refinement::plaza::Option>(m, "RefinementOption")
-      .value("none", dolfinx::refinement::plaza::Option::none)
-      .value("parent_facet", dolfinx::refinement::plaza::Option::parent_facet)
-      .value("parent_cell", dolfinx::refinement::plaza::Option::parent_cell)
+  nb::enum_<dolfinx::refinement::Option>(m, "RefinementOption")
+      .value("none", dolfinx::refinement::Option::none)
+      .value("parent_facet", dolfinx::refinement::Option::parent_facet)
+      .value("parent_cell", dolfinx::refinement::Option::parent_cell)
       .value("parent_cell_and_facet",
-             dolfinx::refinement::plaza::Option::parent_cell_and_facet);
+             dolfinx::refinement::Option::parent_cell_and_facet);
   m.def(
       "transfer_facet_meshtag",
       [](const dolfinx::mesh::MeshTags<std::int32_t>& parent_meshtag,
diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py
index 670a37ca46f..3b82e8d3899 100644
--- a/python/test/unit/refinement/test_interval.py
+++ b/python/test/unit/refinement/test_interval.py
@@ -21,20 +21,27 @@
     [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet],
 )
 @pytest.mark.parametrize("redistribute", [True, False])
-def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined):
+@pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell])
+def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
-    msh_refined, edges = mesh.refine_interval(
-        msh, redistribute=redistribute, ghost_mode=ghost_mode_refined
+    msh_refined, edges, vertices = mesh.refine(
+        msh,
+        redistribute=redistribute,
+        ghost_mode=ghost_mode_refined,
+        option=option,
     )
 
     # vertex count
     assert msh_refined.topology.index_map(0).size_global == 2 * n + 1
 
-    # edge count
-    edge_count = np.array([len(edges)], dtype=np.int64)
-    edge_count = MPI.COMM_WORLD.allreduce(edge_count)
+    if option is mesh.RefinementOption.parent_cell:
+        # edge count
+        edge_count = np.array([len(edges)], dtype=np.int64)
+        edge_count = MPI.COMM_WORLD.allreduce(edge_count)
 
-    assert edge_count == msh_refined.topology.index_map(1).size_global == 2 * n
+        assert edge_count == msh_refined.topology.index_map(1).size_global == 2 * n
+    else:
+        assert edges is None
 
 
 @pytest.mark.parametrize("n", [50, 100])
@@ -48,8 +55,12 @@ def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined):
 @pytest.mark.parametrize("redistribute", [True, False])
 def test_refine_interval_adaptive(n, ghost_mode, redistribute, ghost_mode_refined):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
-    msh_refined, edges = mesh.refine_interval(
-        msh, np.arange(10, dtype=np.int32), redistribute=redistribute, ghost_mode=ghost_mode_refined
+    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,
     )
 
     # vertex count
diff --git a/python/test/unit/mesh/test_refinement.py b/python/test/unit/refinement/test_refinement.py
similarity index 76%
rename from python/test/unit/mesh/test_refinement.py
rename to python/test/unit/refinement/test_refinement.py
index f9415fe2c01..440057ce265 100644
--- a/python/test/unit/mesh/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -24,7 +24,6 @@
     locate_entities_boundary,
     meshtags,
     refine,
-    refine_plaza,
     transfer_meshtag,
 )
 
@@ -33,36 +32,24 @@ 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 = refine(mesh, redistribute=False)
-    assert mesh.topology.index_map(0).size_global == 165
-    assert mesh.topology.index_map(2).size_global == 280
+    mesh_refined, _, _ = refine(mesh, redistribute=False)
+    assert mesh_refined.topology.index_map(0).size_global == 165
+    assert mesh_refined.topology.index_map(2).size_global == 280
 
-
-def test_refine_create_unit_cube_repartition():
-    """Refine mesh of unit cube."""
-    mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=GhostMode.none)
-    mesh.topology.create_entities(1)
-    mesh = refine(mesh, redistribute=True)
-    assert mesh.topology.index_map(0).size_global == 3135
-    assert mesh.topology.index_map(3).size_global == 15120
-
-    mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=GhostMode.shared_facet)
-    mesh.topology.create_entities(1)
-    mesh = refine(mesh, redistribute=True)
-    assert mesh.topology.index_map(0).size_global == 3135
-    assert mesh.topology.index_map(3).size_global == 15120
-
-    Q = functionspace(mesh, ("Lagrange", 1))
-    assert Q
+    # Test that 2D refinement is still 2D
+    assert mesh.geometry.dim == mesh_refined.geometry.dim
 
 
-def test_refine_create_unit_cube_keep_partition():
+@pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet])
+@pytest.mark.parametrize("redistribute", [True, False])
+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=GhostMode.none)
+    mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode)
     mesh.topology.create_entities(1)
-    mesh = refine(mesh, redistribute=False)
+    mesh, _, _ = refine(mesh, redistribute=redistribute)
     assert mesh.topology.index_map(0).size_global == 3135
     assert mesh.topology.index_map(3).size_global == 15120
+
     Q = functionspace(mesh, ("Lagrange", 1))
     assert Q
 
@@ -71,7 +58,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, redistribute=True)
 
     V = functionspace(mesh, ("Lagrange", 1))
 
@@ -82,14 +69,6 @@ def test_refine_create_form():
     assemble_matrix(a)
 
 
-def test_refinement_gdim():
-    """Test that 2D refinement is still 2D"""
-    mesh = create_unit_square(MPI.COMM_WORLD, 3, 4, ghost_mode=GhostMode.none)
-    mesh.topology.create_entities(1)
-    mesh2 = refine(mesh, redistribute=True)
-    assert mesh.geometry.dim == mesh2.geometry.dim
-
-
 def test_sub_refine():
     """Test that refinement of a subset of edges works"""
     mesh = create_unit_square(
@@ -104,7 +83,7 @@ def left_corner_edge(x, tol=1e-7):
     if MPI.COMM_WORLD.size == 0:
         assert edges == 1
 
-    mesh2 = refine(mesh, edges, redistribute=False)
+    mesh2, _, _ = refine(mesh, edges, redistribute=False)
     assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global
 
 
@@ -126,7 +105,7 @@ def left_side(x, tol=1e-16):
     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)
+    mesh2, _, _ = refine(msh, edges, redistribute=True)
 
     num_cells_global = mesh2.topology.index_map(2).size_global
     actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny
@@ -137,12 +116,14 @@ def left_side(x, tol=1e-16):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine_plaza(mesh, None, False, RefinementOption.parent_cell_and_facet),
-        lambda mesh: refine_plaza(
+        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),
-            False,
-            RefinementOption.parent_cell_and_facet,
+            redistribute=False,
+            option=RefinementOption.parent_cell_and_facet,
         ),
     ],
 )
@@ -197,12 +178,14 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine_plaza(mesh, None, False, RefinementOption.parent_cell_and_facet),
-        lambda mesh: refine_plaza(
+        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),
-            False,
-            RefinementOption.parent_cell_and_facet,
+            redistribute=False,
+            option=RefinementOption.parent_cell_and_facet,
         ),
     ],
 )