From a1d0027bdc72ba33af7a0f9cce325ec09ad1cc7b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Wed, 24 Jul 2024 13:37:59 +0200
Subject: [PATCH] Unify refine interface

---
 cpp/dolfinx/refinement/CMakeLists.txt         |   1 +
 cpp/dolfinx/refinement/interval.h             | 104 +++-----
 cpp/dolfinx/refinement/option.h               |  50 ++++
 cpp/dolfinx/refinement/plaza.h                | 252 ++++--------------
 cpp/dolfinx/refinement/refine.h               | 117 ++++----
 cpp/test/CMakeLists.txt                       |   1 +
 cpp/test/refinement/interval.cpp              |  40 +--
 cpp/test/refinement/option.cpp                |  51 ++++
 python/dolfinx/mesh.py                        |  74 +----
 python/dolfinx/wrappers/refinement.cpp        | 114 +++-----
 python/test/unit/refinement/test_interval.py  |  29 +-
 .../{mesh => refinement}/test_refinement.py   |  69 ++---
 12 files changed, 369 insertions(+), 533 deletions(-)
 create mode 100644 cpp/dolfinx/refinement/option.h
 create mode 100644 cpp/test/refinement/option.cpp
 rename python/test/unit/{mesh => refinement}/test_refinement.py (76%)

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 00a8b098d23..05b49e836ea 100644
--- a/cpp/dolfinx/refinement/interval.h
+++ b/cpp/dolfinx/refinement/interval.h
@@ -7,39 +7,44 @@
 #pragma once
 
 #include <algorithm>
-#include <cstddef>
-#include <mpi.h>
-
 #include <concepts>
+#include <cstddef>
 #include <cstdint>
 #include <optional>
 #include <stdexcept>
 #include <vector>
 
-#include "dolfinx/mesh/Mesh.h"
-#include "dolfinx/mesh/cell_types.h"
-#include "dolfinx/mesh/utils.h"
-#include "dolfinx/refinement/plaza.h"
+#include <mpi.h>
 
-namespace dolfinx::refinement
-{
+#include "dolfinx/mesh/Mesh.h"
+#include "dolfinx/refinement/option.h"
+#include "dolfinx/refinement/utils.h"
 
-namespace impl
+namespace dolfinx::refinement::interval
 {
 
 /// 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] 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);
@@ -127,8 +132,12 @@ 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)
   {
@@ -150,18 +159,22 @@ 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() == refined_cell_count * 2);
-  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++; });
@@ -169,56 +182,7 @@ compute_interval_refinement(const mesh::Mesh<T>& mesh,
   graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
 
   return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
-          std::move(parent_cell)};
-}
-
-} // namespace impl
-
-/// Refines 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)};
-  }
+          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..86ed5a30ccf
--- /dev/null
+++ b/cpp/dolfinx/refinement/option.h
@@ -0,0 +1,50 @@
+// 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,
+    parent_facet = 0b01,
+    parent_cell = 0b10,
+    parent_cell_and_facet = 0b11
+};
+
+/// @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;
+}
+}
diff --git a/cpp/dolfinx/refinement/plaza.h b/cpp/dolfinx/refinement/plaza.h
index 9231144a84a..c4e042622fd 100644
--- a/cpp/dolfinx/refinement/plaza.h
+++ b/cpp/dolfinx/refinement/plaza.h
@@ -4,21 +4,23 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
-#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>
 #include <vector>
 
+#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"
+
 #pragma once
 
 /// @brief Plaza mesh refinement.
@@ -29,18 +31,6 @@
 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
@@ -295,30 +285,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 +369,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 +411,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 +419,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 +448,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 +459,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 +500,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 +539,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..42753c434ab 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -6,74 +6,87 @@
 
 #pragma once
 
+#include <algorithm>
+#include <concepts>
+#include <optional>
+#include <utility>
+
+#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>
 
 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
-template <std::floating_point T>
-mesh::Mesh<T> refine(const mesh::Mesh<T>& mesh, bool redistribute = true)
-{
-  auto topology = mesh.topology();
-  assert(topology);
-
-  if (topology->cell_type() != mesh::CellType::triangle
-      and topology->cell_type() != mesh::CellType::tetrahedron)
-  {
-    throw std::runtime_error("Refinement only defined for simplices");
-  }
 
-  auto [refined_mesh, parent_cell, parent_facet]
-      = plaza::refine(mesh, redistribute, plaza::Option::none);
+namespace
+{
 
-  // 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));
+/// @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,
+                    const std::vector<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);
 
-  return refined_mesh;
+  else
+    // let partition handle the parallel construction of the mesh
+    return partition<T>(mesh, cell_adj, new_vertex_coords, xshape, redistribute,
+                        ghost_mode);
 }
+} // namespace
 
-/// @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 performend, 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 (cell_t != mesh::CellType::interval and cell_t != mesh::CellType::triangle
+      and cell_t != mesh::CellType::tetrahedron)
     throw std::runtime_error("Refinement only defined for simplices");
-  }
 
-  auto [refined_mesh, parent_cell, parent_facet]
-      = plaza::refine(mesh, edges, redistribute, plaza::Option::none);
+  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);
+
+  mesh::Mesh<T> refined_mesh = create_refined_mesh(
+      mesh, std::move(cell_adj), std::move(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 +94,7 @@ 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 {refined_mesh, std::move(parent_cell), std::move(parent_facet)};
 }
 
 } // namespace dolfinx::refinement
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index e652537d80c..fc91064cb0a 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -44,6 +44,7 @@ add_executable(
   mesh/distributed_mesh.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..fcaf2ab24e3
--- /dev/null
+++ b/cpp/test/refinement/option.cpp
@@ -0,0 +1,51 @@
+// 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 <cstdint>
+
+#include <catch2/catch_test_macros.hpp>
+
+#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 45770c1c476..f91907ad5f3 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -4,20 +4,27 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
-#include "array.h"
 #include <concepts>
+#include <cstdint>
+#include <optional>
+#include <span>
+
+#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 <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/shared_ptr.h>
-#include <nanobind/stl/tuple.h>
-#include <nanobind/stl/vector.h>
-#include <optional>
+
+#include "array.h"
 
 namespace nb = nanobind;
 
@@ -27,75 +34,38 @@ 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)
@@ -103,12 +73,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,
         ),
     ],
 )