From e1ed42ad55951be975cdf5481ba86f0ee0289e61 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Wed, 21 Aug 2024 13:13:07 +0200
Subject: [PATCH 01/31] Remove refine option `distribute` and `ghost_mode` in
 favor of an explicit partitioner

---
 cpp/dolfinx/refinement/refine.h               |  67 ++--
 cpp/dolfinx/refinement/utils.h                |  45 ---
 cpp/test/mesh/refinement/interval.cpp         |   7 +-
 cpp/test/mesh/refinement/rectangle.cpp        |   6 +-
 cpp/test/transfer/transfer_matrix.cpp         | 315 ++++++++++++++++++
 python/dolfinx/mesh.py                        |  12 +-
 python/dolfinx/wrappers/refinement.cpp        |  78 ++++-
 python/test/unit/refinement/test_interval.py  |  12 +-
 .../test/unit/refinement/test_refinement.py   |  20 +-
 9 files changed, 434 insertions(+), 128 deletions(-)
 create mode 100644 cpp/test/transfer/transfer_matrix.cpp

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 32f18d89827..50d2181207c 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -17,46 +17,40 @@
 #include <optional>
 #include <utility>
 
+#include "dolfinx/graph/AdjacencyList.h"
+#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"
+
 namespace dolfinx::refinement
 {
-namespace impl
-{
-/// @brief create the refined mesh by optionally redistributing it
-template <std::floating_point T>
-mesh::Mesh<T>
-create_refined_mesh(const mesh::Mesh<T>& mesh,
-                    const graph::AdjacencyList<std::int64_t>& cell_adj,
-                    std::span<const T> new_vertex_coords,
-                    std::array<std::size_t, 2> xshape, bool redistribute,
-                    mesh::GhostMode ghost_mode)
+
+// TODO: move to cpp?
+inline graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
+    MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
+    const std::vector<std::span<const std::int64_t>>& cell_topology)
 {
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    // No parallel construction necessary, i.e. redistribute also has no
-    // effect
-    return mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                             mesh.geometry().cmap(), new_vertex_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    // Let partition handle the parallel construction of the mesh
-    return partition<T>(mesh, cell_adj, new_vertex_coords, xshape, redistribute,
-                        ghost_mode);
-  }
+  std::int32_t mpi_rank = MPI::rank(comm);
+  std::int32_t num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
+  std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
+  std::vector<std::int32_t> destinations(num_cells, mpi_rank);
+  std::vector<std::int32_t> dest_offsets(num_cells + 1);
+  std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
+  return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
 }
-} // namespace impl
 
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
 ///
 /// @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] 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] partitioner partitioner to be used for 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
@@ -65,8 +59,8 @@ template <std::floating_point T>
 std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
            std::optional<std::vector<std::int8_t>>>
 refine(const mesh::Mesh<T>& mesh,
-       std::optional<std::span<const std::int32_t>> edges, bool redistribute,
-       mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet,
+       std::optional<std::span<const std::int32_t>> edges,
+       mesh::CellPartitionFunction partitioner = maintain_coarse_partitioner,
        Option option = Option::none)
 {
   auto topology = mesh.topology();
@@ -80,9 +74,10 @@ refine(const mesh::Mesh<T>& mesh,
       = oned ? interval::compute_refinement_data(mesh, edges, option)
              : plaza::compute_refinement_data(mesh, edges, option);
 
-  mesh::Mesh<T> refined_mesh = impl::create_refined_mesh(
-      mesh, cell_adj, std::span<const T>(new_vertex_coords), xshape,
-      redistribute, ghost_mode);
+  mesh::Mesh<T> refined_mesh
+      = mesh::create_mesh(mesh.comm(), mesh.comm(), std::move(cell_adj).array(),
+                          mesh.geometry().cmap(), mesh.comm(),
+                          std::move(new_vertex_coords), xshape, partitioner);
 
   // Report the number of refined cellse
   const int D = topology->dim();
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index 275b40d7cc2..1c6af077551 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -264,51 +264,6 @@ create_new_vertices(MPI_Comm comm,
           xshape};
 }
 
-/// Use vertex and topology data to partition new mesh across
-/// processes
-/// @param[in] old_mesh
-/// @param[in] cell_topology Topology of cells, (vertex indices)
-/// @param[in] new_coords New coordinates, row-major storage
-/// @param[in] xshape The shape of `new_coords`
-/// @param[in] redistribute Call graph partitioner if true
-/// @param[in] ghost_mode None or shared_facet
-/// @return New mesh
-template <std::floating_point T>
-mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
-                        const graph::AdjacencyList<std::int64_t>& cell_topology,
-                        std::span<const T> new_coords,
-                        std::array<std::size_t, 2> xshape, bool redistribute,
-                        mesh::GhostMode ghost_mode)
-{
-  if (redistribute)
-  {
-    return mesh::create_mesh(old_mesh.comm(), cell_topology.array(),
-                             old_mesh.geometry().cmap(), new_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    auto partitioner
-        = [](MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-             const std::vector<std::span<const std::int64_t>>& cell_topology)
-    {
-      std::int32_t mpi_rank = MPI::rank(comm);
-      std::int32_t num_cell_vertices
-          = mesh::num_cell_vertices(cell_types.front());
-      std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
-      std::vector<std::int32_t> destinations(num_cells, mpi_rank);
-      std::vector<std::int32_t> dest_offsets(num_cells + 1);
-      std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
-      return graph::AdjacencyList(std::move(destinations),
-                                  std::move(dest_offsets));
-    };
-
-    return mesh::create_mesh(old_mesh.comm(), old_mesh.comm(),
-                             cell_topology.array(), old_mesh.geometry().cmap(),
-                             old_mesh.comm(), new_coords, xshape, partitioner);
-  }
-}
-
 /// @brief Given an index map, add "n" extra indices at the end of local range
 ///
 /// @note The returned global indices (local and ghosts) are adjust for the
diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp
index 763418ac4b3..74b3f434627 100644
--- a/cpp/test/mesh/refinement/interval.cpp
+++ b/cpp/test/mesh/refinement/interval.cpp
@@ -73,7 +73,7 @@ TEMPLATE_TEST_CASE("Interval uniform refinement",
 
   // TODO: parent_facet
   auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
-      mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
+      mesh, std::nullopt, &refinement::maintain_coarse_partitioner,
       refinement::Option::parent_cell);
 
   std::vector<T> expected_x = {
@@ -114,7 +114,8 @@ TEMPLATE_TEST_CASE("Interval adaptive refinement",
   std::vector<std::int32_t> edges{1};
   // TODO: parent_facet
   auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
-      mesh, std::span(edges), false, mesh::GhostMode::shared_facet,
+      mesh, std::span(edges),
+      mesh::create_cell_partitioner(mesh::GhostMode::shared_facet),
       refinement::Option::parent_cell);
 
   std::vector<T> expected_x = {
@@ -192,7 +193,7 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)",
 
   // TODO: parent_facet
   auto [refined_mesh, parent_edges, parent_facet] = refinement::refine(
-      mesh, std::nullopt, false, mesh::GhostMode::shared_facet,
+      mesh, std::nullopt, &refinement::maintain_coarse_partitioner,
       refinement::Option::parent_cell);
 
   T rank_d = static_cast<T>(rank);
diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp
index 954879abb1d..ff2a67979df 100644
--- a/cpp/test/mesh/refinement/rectangle.cpp
+++ b/cpp/test/mesh/refinement/rectangle.cpp
@@ -99,9 +99,9 @@ plotter.show()
   // plaza requires the edges to be pre initialized!
   mesh.topology()->create_entities(1);
 
-  auto [mesh_fine, parent_cell, parent_facet]
-      = refinement::refine(mesh, std::nullopt, false, mesh::GhostMode::none,
-                           refinement::Option::parent_cell_and_facet);
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      mesh, std::nullopt, mesh::create_cell_partitioner(mesh::GhostMode::none),
+      refinement::Option::parent_cell_and_facet);
 
   // vertex layout:
   // 8---7---5
diff --git a/cpp/test/transfer/transfer_matrix.cpp b/cpp/test/transfer/transfer_matrix.cpp
new file mode 100644
index 00000000000..9747e776fcb
--- /dev/null
+++ b/cpp/test/transfer/transfer_matrix.cpp
@@ -0,0 +1,315 @@
+// 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 <concepts>
+#include <cstdint>
+#include <memory>
+#include <optional>
+
+#include <mpi.h>
+
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
+
+#include <dolfinx/fem/FunctionSpace.h>
+#include <dolfinx/fem/utils.h>
+#include <dolfinx/la/MatrixCSR.h>
+#include <dolfinx/la/SparsityPattern.h>
+#include <dolfinx/mesh/cell_types.h>
+#include <dolfinx/mesh/generation.h>
+#include <dolfinx/mesh/utils.h>
+#include <dolfinx/refinement/plaza.h>
+#include <dolfinx/refinement/refine.h>
+#include <dolfinx/transfer/transfer_matrix.h>
+
+#include "../mesh/util.h"
+
+using namespace dolfinx;
+using namespace Catch::Matchers;
+
+template <std::floating_point T>
+constexpr auto weight = [](std::int32_t d) -> T { return d == 0 ? 1. : .5; };
+
+TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]",
+                   double) // TODO: float
+{
+  using T = TestType;
+
+  auto mesh_coarse
+      = std::make_shared<mesh::Mesh<T>>(dolfinx::mesh::create_interval<T>(
+          MPI_COMM_SELF, 2, {0.0, 1.0}, mesh::GhostMode::none));
+
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
+      refinement::Option::parent_cell);
+
+  auto element = basix::create_element<T>(
+      basix::element::family::P, basix::cell::type::interval, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
+  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
+      fem::create_functionspace<T>(mesh_coarse, element, {}));
+  auto V_fine
+      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
+
+  mesh_fine.topology()->create_connectivity(1, 0);
+  mesh_fine.topology()->create_connectivity(0, 1);
+
+  std::vector<std::int64_t> from_to_map{0, 2, 4}; // TODO: general computation!
+
+  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix<T>(
+      *V_coarse, *V_fine, from_to_map, weight<T>);
+
+  std::vector<T> expected{1.0, .5, 0, 0, 0, 0, .5, 1, .5, 0, 0, 0, 0, .5, 1};
+  CHECK_THAT(transfer_matrix.to_dense(),
+             RangeEquals(expected, [](auto a, auto b)
+                         { return std::abs(a - b) <= EPS<T>; }));
+}
+
+TEMPLATE_TEST_CASE("Transfer Matrix 1D (parallel)", "[transfer_matrix]",
+                   double) // TODO: float
+{
+  using T = TestType;
+
+  int comm_size = dolfinx::MPI::size(MPI_COMM_WORLD);
+  int rank = dolfinx::MPI::rank(MPI_COMM_WORLD);
+
+  // TODO: see https://github.com/FEniCS/dolfinx/issues/3358
+  // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
+  mesh::CellPartitionFunction part
+      = [&](MPI_Comm /* comm */, int /* nparts */,
+            const std::vector<mesh::CellType>& /* cell_types */,
+            const std::vector<std::span<const std::int64_t>>& /* cells */)
+  {
+    std::vector<std::vector<std::int32_t>> data;
+    if (comm_size == 1)
+      data = {{0}, {0}, {0}, {0}};
+    else if (comm_size == 2)
+      data = {{0}, {0}, {0}, {0}, {0, 1}, {1, 0}, {1}, {1}, {1}};
+    else if (comm_size == 3)
+      data = {{1}, {1},    {1},    {1}, {1, 2}, {2, 1}, {2},
+              {2}, {2, 0}, {0, 2}, {0}, {0},    {0},    {0}};
+    else
+      FAIL("Test only supports <= 3 processes");
+    return graph::AdjacencyList<std::int32_t>(std::move(data));
+  };
+
+  auto mesh_coarse = std::make_shared<mesh::Mesh<T>>(
+      mesh::create_interval<T>(MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.},
+                               mesh::GhostMode::shared_facet, part));
+
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      *mesh_coarse, std::nullopt, part, mesh::GhostMode::shared_facet,
+      refinement::Option::none);
+
+  auto element = basix::create_element<T>(
+      basix::element::family::P, basix::cell::type::interval, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
+  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
+      fem::create_functionspace<T>(mesh_coarse, element, {}));
+  auto V_fine
+      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
+
+  mesh_fine.topology()->create_connectivity(1, 0);
+  mesh_fine.topology()->create_connectivity(0, 1);
+
+  std::vector<int64_t> from_to_map;
+  switch (comm_size)
+  {
+  case 1:
+    from_to_map = {0, 2, 4, 6, 8};
+    break;
+  case 2:
+    from_to_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18};
+    break;
+  case 3:
+    from_to_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 28, 26, 24, 22};
+    break;
+  }
+
+  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
+      *V_coarse, *V_fine, from_to_map, weight<T>);
+
+  std::array<std::vector<T>, 3> expected;
+  if (comm_size == 1)
+  {
+    // clang-format off
+    expected[0] = {
+        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0
+    };
+    // clang-format on
+  }
+  else if (comm_size == 2)
+  {
+    // clang-format off
+    expected[0] = {
+        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+    };
+    expected[1] = {
+        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+    };
+    // clang-format on
+  }
+  else if (comm_size == 3)
+  {
+    // clang-format off
+    expected[0] = {
+        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    };
+    expected[1] = {
+        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5,
+        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    };
+    expected[2] = {
+        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0,
+        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
+        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
+        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    };
+    // clang-format on
+  }
+  else
+  {
+    CHECK(false);
+  }
+
+  CHECK_THAT(transfer_matrix.to_dense(),
+             RangeEquals(expected[rank], [](auto a, auto b)
+                         { return std::abs(a - b) <= EPS<T>; }));
+}
+
+TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double)
+{
+  using T = TestType;
+
+  auto mesh_coarse
+      = std::make_shared<mesh::Mesh<T>>(dolfinx::mesh::create_rectangle<T>(
+          MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle));
+  mesh_coarse->topology()->create_entities(1);
+
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
+      refinement::Option::parent_cell);
+
+  auto element = basix::create_element<T>(
+      basix::element::family::P, basix::cell::type::triangle, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
+  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
+      fem::create_functionspace<T>(mesh_coarse, element, {}));
+  auto V_fine
+      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
+
+  mesh_fine.topology()->create_connectivity(1, 0);
+  mesh_fine.topology()->create_connectivity(0, 1);
+
+  std::vector<std::int64_t> from_to_map{4, 1, 5,
+                                        8}; // TODO: general computation!
+
+  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
+      *V_coarse, *V_fine, from_to_map, weight<T>);
+
+  std::vector<T> expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0,
+                          0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
+                          0.5, 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0,
+                          0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0};
+  CHECK_THAT(transfer_matrix.to_dense(),
+             RangeEquals(expected, [](auto a, auto b)
+                         { return std::abs(a - b) <= EPS<T>; }));
+}
+
+TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double)
+{
+  using T = TestType;
+
+  auto mesh_coarse = std::make_shared<mesh::Mesh<T>>(
+      dolfinx::mesh::create_box<T>(MPI_COMM_SELF, {{{0, 0, 0}, {1, 1, 1}}},
+                                   {1, 1, 1}, mesh::CellType::tetrahedron));
+  mesh_coarse->topology()->create_entities(1);
+
+  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
+      refinement::Option::parent_cell);
+
+  auto element = basix::create_element<T>(
+      basix::element::family::P, basix::cell::type::tetrahedron, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
+  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
+      fem::create_functionspace<T>(mesh_coarse, element, {}));
+  auto V_fine
+      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
+
+  mesh_fine.topology()->create_connectivity(1, 0);
+  mesh_fine.topology()->create_connectivity(0, 1);
+
+  std::vector<std::int64_t> from_to_map{
+      0, 6, 15, 25, 17, 9, 11, 22}; // TODO: general computation!
+
+  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
+      *V_coarse, *V_fine, from_to_map, weight<T>);
+
+  std::vector<T> expected{
+      1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5,
+      0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0,
+      0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5,
+      0.5, 1.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.5, 0.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0,
+      0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0,
+      0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
+      0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
+      0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5,
+      0.5, 1.0, 0.0, 0.0, 0.0, 0.5};
+  CHECK_THAT(transfer_matrix.to_dense(),
+             RangeEquals(expected, [](auto a, auto b)
+                         { return std::abs(a - b) <= EPS<T>; }));
+}
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index b9cefce09b3..46fec0bed34 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -370,8 +370,7 @@ def transfer_meshtag(
 def refine(
     mesh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
-    redistribute: bool = True,
-    ghost_mode: GhostMode = GhostMode.shared_facet,
+    partitioner: typing.Optional[typing.Callable] = None,
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
@@ -380,17 +379,16 @@ def refine(
         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``.
-        ghost_mode: Ghost mode to use for the refined mesh.
+        partitioner:
+            partitioner to use for the refined mesh, If ``None`` no redistribution is performed,
+            i.e. previous local mesh is equally parallelized now with new vertices.s
         option: Controls whether parent cells and/or parent facets are computed.
 
-
     Returns:
        Refined mesh, (optional) parent cells, (optional) parent facets
     """
     mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
-        mesh._cpp_object, edges, redistribute, ghost_mode, option
+        mesh._cpp_object, edges, partitioner, option
     )
     return Mesh(mesh1, mesh._ufl_domain), parent_cell, parent_facet
 
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 37ca3937992..7921ffb22e5 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -4,9 +4,20 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
-#include "array.h"
 #include <concepts>
 #include <cstdint>
+#include <functional>
+#include <optional>
+#include <span>
+
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/array.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>
@@ -14,30 +25,67 @@
 #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>
+
+#include "MPICommWrapper.h"
+#include "array.h"
 
 namespace nb = nanobind;
 
 namespace dolfinx_wrappers
 {
 
+namespace
+{
+
+using PythonCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        dolfinx_wrappers::MPICommWrapper, int,
+        const std::vector<dolfinx::mesh::CellType>&,
+        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
+
+using CppCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
+        const std::vector<std::span<const std::int64_t>>&)>;
+
+/// Wrap a Python graph partitioning function as a C++ function
+CppCellPartitionFunction
+create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
+{
+  if (p)
+  {
+    return [p](MPI_Comm comm, int n,
+               const std::vector<dolfinx::mesh::CellType>& cell_types,
+               const std::vector<std::span<const std::int64_t>>& cells)
+    {
+      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
+      std::ranges::transform(
+          cells, std::back_inserter(cells_nb),
+          [](auto c)
+          {
+            return nb::ndarray<const std::int64_t, nb::numpy>(
+                c.data(), {c.size()}, nb::handle());
+          });
+
+      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
+    };
+  }
+  else
+    return nullptr;
+}
+} // namespace
+
 template <std::floating_point T>
 void export_refinement_with_variable_mesh_type(nb::module_& m)
 {
+
   m.def(
       "refine",
       [](const dolfinx::mesh::Mesh<T>& mesh,
          std::optional<
              nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
              edges,
-         bool redistribute, dolfinx::mesh::GhostMode ghost_mode,
+         std::optional<PythonCellPartitionFunction> partitioner,
          dolfinx::refinement::Option option)
       {
         std::optional<std::span<const std::int32_t>> cpp_edges(std::nullopt);
@@ -47,8 +95,12 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
               std::span(edges.value().data(), edges.value().size()));
         }
 
+        CppCellPartitionFunction cpp_partitioner
+            = partitioner.has_value()
+                  ? create_cell_partitioner_cpp(partitioner.value())
+                  : dolfinx::refinement::maintain_coarse_partitioner;
         auto [mesh1, cell, facet] = dolfinx::refinement::refine(
-            mesh, cpp_edges, redistribute, ghost_mode, option);
+            mesh, cpp_edges, cpp_partitioner, option);
 
         std::optional<nb::ndarray<std::int32_t, nb::numpy>> python_cell(
             std::nullopt);
@@ -63,8 +115,8 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
         return std::tuple{std::move(mesh1), std::move(python_cell),
                           std::move(python_facet)};
       },
-      nb::arg("mesh"), nb::arg("edges") = nb::none(), nb::arg("redistribute"),
-      nb::arg("ghost_mode"), nb::arg("option"));
+      nb::arg("mesh"), nb::arg("edges") = nb::none(),
+      nb::arg("partitioner") = nb::none(), nb::arg("option"));
 }
 
 void refinement(nb::module_& m)
diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py
index 28779644840..27040fb0b90 100644
--- a/python/test/unit/refinement/test_interval.py
+++ b/python/test/unit/refinement/test_interval.py
@@ -20,16 +20,14 @@
     "ghost_mode_refined",
     [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet],
 )
-@pytest.mark.parametrize("redistribute", [True, False])
 @pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell])
-def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option):
+def test_refine_interval(n, ghost_mode, ghost_mode_refined, option):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
     msh_refined, edges, vertices = mesh.refine(
         msh,
-        redistribute=redistribute,
-        ghost_mode=ghost_mode_refined,
         option=option,
     )
+    # TODO: add create_cell_partitioner(ghost_mode) when works
 
     # vertex count
     assert msh_refined.topology.index_map(0).size_global == 2 * n + 1
@@ -52,16 +50,14 @@ def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option
     "ghost_mode_refined",
     [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet],
 )
-@pytest.mark.parametrize("redistribute", [True, False])
-def test_refine_interval_adaptive(n, ghost_mode, redistribute, ghost_mode_refined):
+def test_refine_interval_adaptive(n, ghost_mode, ghost_mode_refined):
     msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode)
     msh_refined, edges, vertices = mesh.refine(
         msh,
         np.arange(10, dtype=np.int32),
-        redistribute=redistribute,
-        ghost_mode=ghost_mode_refined,
         option=mesh.RefinementOption.parent_cell,
     )
+    # TODO: add create_cell_partitioner(ghost_mode) when works
 
     # vertex count
     assert msh_refined.topology.index_map(0).size_global == n + 1 + 10 * MPI.COMM_WORLD.size
diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 440057ce265..eb5188d5bf1 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -32,7 +32,7 @@ def test_refine_create_unit_square():
     """Refine mesh of unit square."""
     mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none)
     mesh.topology.create_entities(1)
-    mesh_refined, _, _ = refine(mesh, redistribute=False)
+    mesh_refined, _, _ = refine(mesh)
     assert mesh_refined.topology.index_map(0).size_global == 165
     assert mesh_refined.topology.index_map(2).size_global == 280
 
@@ -46,7 +46,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute):
     """Refine mesh of unit cube."""
     mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode)
     mesh.topology.create_entities(1)
-    mesh, _, _ = refine(mesh, redistribute=redistribute)
+    mesh, _, _ = refine(mesh)  # , partitioner=create_cell_partitioner(ghost_mode))
     assert mesh.topology.index_map(0).size_global == 3135
     assert mesh.topology.index_map(3).size_global == 15120
 
@@ -58,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)  # TODO: recover redistribute=True behavior
 
     V = functionspace(mesh, ("Lagrange", 1))
 
@@ -83,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)
     assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global
 
 
@@ -105,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)  # TODO: redistribute=True
 
     num_cells_global = mesh2.topology.index_map(2).size_global
     actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny
@@ -116,13 +116,10 @@ def left_side(x, tol=1e-16):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(
-            mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet
-        ),
+        lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet),
         lambda mesh: refine(
             mesh,
             np.arange(mesh.topology.index_map(1).size_local),
-            redistribute=False,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],
@@ -178,13 +175,10 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(
-            mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet
-        ),
+        lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet),
         lambda mesh: refine(
             mesh,
             np.arange(mesh.topology.index_map(1).size_local),
-            redistribute=False,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],

From 32fca7cb63bbf64732a24541066f61ca254f20df Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 17 Sep 2024 22:43:48 +0200
Subject: [PATCH 02/31] Remove wrongly added transfer_matrix -> later PR

---
 cpp/test/transfer/transfer_matrix.cpp | 315 --------------------------
 1 file changed, 315 deletions(-)
 delete mode 100644 cpp/test/transfer/transfer_matrix.cpp

diff --git a/cpp/test/transfer/transfer_matrix.cpp b/cpp/test/transfer/transfer_matrix.cpp
deleted file mode 100644
index 9747e776fcb..00000000000
--- a/cpp/test/transfer/transfer_matrix.cpp
+++ /dev/null
@@ -1,315 +0,0 @@
-// 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 <concepts>
-#include <cstdint>
-#include <memory>
-#include <optional>
-
-#include <mpi.h>
-
-#include <catch2/catch_template_test_macros.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers.hpp>
-#include <catch2/matchers/catch_matchers_floating_point.hpp>
-#include <catch2/matchers/catch_matchers_range_equals.hpp>
-
-#include <dolfinx/fem/FunctionSpace.h>
-#include <dolfinx/fem/utils.h>
-#include <dolfinx/la/MatrixCSR.h>
-#include <dolfinx/la/SparsityPattern.h>
-#include <dolfinx/mesh/cell_types.h>
-#include <dolfinx/mesh/generation.h>
-#include <dolfinx/mesh/utils.h>
-#include <dolfinx/refinement/plaza.h>
-#include <dolfinx/refinement/refine.h>
-#include <dolfinx/transfer/transfer_matrix.h>
-
-#include "../mesh/util.h"
-
-using namespace dolfinx;
-using namespace Catch::Matchers;
-
-template <std::floating_point T>
-constexpr auto weight = [](std::int32_t d) -> T { return d == 0 ? 1. : .5; };
-
-TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]",
-                   double) // TODO: float
-{
-  using T = TestType;
-
-  auto mesh_coarse
-      = std::make_shared<mesh::Mesh<T>>(dolfinx::mesh::create_interval<T>(
-          MPI_COMM_SELF, 2, {0.0, 1.0}, mesh::GhostMode::none));
-
-  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
-      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
-      refinement::Option::parent_cell);
-
-  auto element = basix::create_element<T>(
-      basix::element::family::P, basix::cell::type::interval, 1,
-      basix::element::lagrange_variant::unset,
-      basix::element::dpc_variant::unset, false);
-
-  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace<T>(mesh_coarse, element, {}));
-  auto V_fine
-      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
-          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
-
-  mesh_fine.topology()->create_connectivity(1, 0);
-  mesh_fine.topology()->create_connectivity(0, 1);
-
-  std::vector<std::int64_t> from_to_map{0, 2, 4}; // TODO: general computation!
-
-  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix<T>(
-      *V_coarse, *V_fine, from_to_map, weight<T>);
-
-  std::vector<T> expected{1.0, .5, 0, 0, 0, 0, .5, 1, .5, 0, 0, 0, 0, .5, 1};
-  CHECK_THAT(transfer_matrix.to_dense(),
-             RangeEquals(expected, [](auto a, auto b)
-                         { return std::abs(a - b) <= EPS<T>; }));
-}
-
-TEMPLATE_TEST_CASE("Transfer Matrix 1D (parallel)", "[transfer_matrix]",
-                   double) // TODO: float
-{
-  using T = TestType;
-
-  int comm_size = dolfinx::MPI::size(MPI_COMM_WORLD);
-  int rank = dolfinx::MPI::rank(MPI_COMM_WORLD);
-
-  // TODO: see https://github.com/FEniCS/dolfinx/issues/3358
-  // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
-  mesh::CellPartitionFunction part
-      = [&](MPI_Comm /* comm */, int /* nparts */,
-            const std::vector<mesh::CellType>& /* cell_types */,
-            const std::vector<std::span<const std::int64_t>>& /* cells */)
-  {
-    std::vector<std::vector<std::int32_t>> data;
-    if (comm_size == 1)
-      data = {{0}, {0}, {0}, {0}};
-    else if (comm_size == 2)
-      data = {{0}, {0}, {0}, {0}, {0, 1}, {1, 0}, {1}, {1}, {1}};
-    else if (comm_size == 3)
-      data = {{1}, {1},    {1},    {1}, {1, 2}, {2, 1}, {2},
-              {2}, {2, 0}, {0, 2}, {0}, {0},    {0},    {0}};
-    else
-      FAIL("Test only supports <= 3 processes");
-    return graph::AdjacencyList<std::int32_t>(std::move(data));
-  };
-
-  auto mesh_coarse = std::make_shared<mesh::Mesh<T>>(
-      mesh::create_interval<T>(MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.},
-                               mesh::GhostMode::shared_facet, part));
-
-  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
-      *mesh_coarse, std::nullopt, part, mesh::GhostMode::shared_facet,
-      refinement::Option::none);
-
-  auto element = basix::create_element<T>(
-      basix::element::family::P, basix::cell::type::interval, 1,
-      basix::element::lagrange_variant::unset,
-      basix::element::dpc_variant::unset, false);
-
-  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace<T>(mesh_coarse, element, {}));
-  auto V_fine
-      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
-          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
-
-  mesh_fine.topology()->create_connectivity(1, 0);
-  mesh_fine.topology()->create_connectivity(0, 1);
-
-  std::vector<int64_t> from_to_map;
-  switch (comm_size)
-  {
-  case 1:
-    from_to_map = {0, 2, 4, 6, 8};
-    break;
-  case 2:
-    from_to_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18};
-    break;
-  case 3:
-    from_to_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 28, 26, 24, 22};
-    break;
-  }
-
-  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
-      *V_coarse, *V_fine, from_to_map, weight<T>);
-
-  std::array<std::vector<T>, 3> expected;
-  if (comm_size == 1)
-  {
-    // clang-format off
-    expected[0] = {
-        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0
-    };
-    // clang-format on
-  }
-  else if (comm_size == 2)
-  {
-    // clang-format off
-    expected[0] = {
-        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
-    };
-    expected[1] = {
-        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
-    };
-    // clang-format on
-  }
-  else if (comm_size == 3)
-  {
-    // clang-format off
-    expected[0] = {
-        /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-    };
-    expected[1] = {
-        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5,
-        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-    };
-    expected[2] = {
-        /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0,
-        /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0,
-        /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0,
-        /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-    };
-    // clang-format on
-  }
-  else
-  {
-    CHECK(false);
-  }
-
-  CHECK_THAT(transfer_matrix.to_dense(),
-             RangeEquals(expected[rank], [](auto a, auto b)
-                         { return std::abs(a - b) <= EPS<T>; }));
-}
-
-TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double)
-{
-  using T = TestType;
-
-  auto mesh_coarse
-      = std::make_shared<mesh::Mesh<T>>(dolfinx::mesh::create_rectangle<T>(
-          MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle));
-  mesh_coarse->topology()->create_entities(1);
-
-  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
-      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
-      refinement::Option::parent_cell);
-
-  auto element = basix::create_element<T>(
-      basix::element::family::P, basix::cell::type::triangle, 1,
-      basix::element::lagrange_variant::unset,
-      basix::element::dpc_variant::unset, false);
-
-  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace<T>(mesh_coarse, element, {}));
-  auto V_fine
-      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
-          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
-
-  mesh_fine.topology()->create_connectivity(1, 0);
-  mesh_fine.topology()->create_connectivity(0, 1);
-
-  std::vector<std::int64_t> from_to_map{4, 1, 5,
-                                        8}; // TODO: general computation!
-
-  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
-      *V_coarse, *V_fine, from_to_map, weight<T>);
-
-  std::vector<T> expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0,
-                          0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
-                          0.5, 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0,
-                          0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0};
-  CHECK_THAT(transfer_matrix.to_dense(),
-             RangeEquals(expected, [](auto a, auto b)
-                         { return std::abs(a - b) <= EPS<T>; }));
-}
-
-TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double)
-{
-  using T = TestType;
-
-  auto mesh_coarse = std::make_shared<mesh::Mesh<T>>(
-      dolfinx::mesh::create_box<T>(MPI_COMM_SELF, {{{0, 0, 0}, {1, 1, 1}}},
-                                   {1, 1, 1}, mesh::CellType::tetrahedron));
-  mesh_coarse->topology()->create_entities(1);
-
-  auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
-      *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none,
-      refinement::Option::parent_cell);
-
-  auto element = basix::create_element<T>(
-      basix::element::family::P, basix::cell::type::tetrahedron, 1,
-      basix::element::lagrange_variant::unset,
-      basix::element::dpc_variant::unset, false);
-
-  auto V_coarse = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace<T>(mesh_coarse, element, {}));
-  auto V_fine
-      = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
-          std::make_shared<mesh::Mesh<T>>(mesh_fine), element, {}));
-
-  mesh_fine.topology()->create_connectivity(1, 0);
-  mesh_fine.topology()->create_connectivity(0, 1);
-
-  std::vector<std::int64_t> from_to_map{
-      0, 6, 15, 25, 17, 9, 11, 22}; // TODO: general computation!
-
-  la::MatrixCSR<T> transfer_matrix = transfer::create_transfer_matrix(
-      *V_coarse, *V_fine, from_to_map, weight<T>);
-
-  std::vector<T> expected{
-      1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5,
-      0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0,
-      0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5,
-      0.5, 1.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-      0.5, 0.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0,
-      0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0,
-      0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
-      0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
-      0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5,
-      0.5, 1.0, 0.0, 0.0, 0.0, 0.5};
-  CHECK_THAT(transfer_matrix.to_dense(),
-             RangeEquals(expected, [](auto a, auto b)
-                         { return std::abs(a - b) <= EPS<T>; }));
-}

From 5bd8ba5387c894bd993c2de87747baca345394cf Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 17 Sep 2024 22:46:02 +0200
Subject: [PATCH 03/31] Add author name

---
 cpp/dolfinx/refinement/refine.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 50d2181207c..12d5ccce196 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2010-2023 Garth N. Wells
+// Copyright (C) 2010-2024 Garth N. Wells and Paul T. Kühner
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //

From 124630aea480a1823c43e1dcc51371bba1c5efec Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Wed, 18 Sep 2024 11:49:28 +0100
Subject: [PATCH 04/31] Doc fix

---
 cpp/dolfinx/refinement/refine.h | 29 ++++++++++-------------------
 1 file changed, 10 insertions(+), 19 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 12d5ccce196..4e72c35837c 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -6,6 +6,7 @@
 
 #pragma once
 
+#include "dolfinx/graph/AdjacencyList.h"
 #include "dolfinx/mesh/Mesh.h"
 #include "dolfinx/mesh/Topology.h"
 #include "dolfinx/mesh/cell_types.h"
@@ -17,25 +18,16 @@
 #include <optional>
 #include <utility>
 
-#include "dolfinx/graph/AdjacencyList.h"
-#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"
-
 namespace dolfinx::refinement
 {
 
-// TODO: move to cpp?
+/// TODO: Document function
 inline graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
     MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
     const std::vector<std::span<const std::int64_t>>& cell_topology)
 {
-  std::int32_t mpi_rank = MPI::rank(comm);
-  std::int32_t num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
+  int mpi_rank = MPI::rank(comm);
+  int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
   std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
   std::vector<std::int32_t> destinations(num_cells, mpi_rank);
   std::vector<std::int32_t> dest_offsets(num_cells + 1);
@@ -46,10 +38,10 @@ inline graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
 ///
-/// @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] 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] partitioner partitioner to be used for the refined mesh
 /// @param[in] option Control the computation of parent facets, parent
 /// cells. If an option is unselected, an empty list is returned.
@@ -65,10 +57,9 @@ refine(const mesh::Mesh<T>& mesh,
 {
   auto topology = mesh.topology();
   assert(topology);
-
-  mesh::CellType cell_t = topology->cell_type();
-  if (!mesh::is_simplex(cell_t))
+  if (!mesh::is_simplex(topology->cell_type()))
     throw std::runtime_error("Refinement only defined for simplices");
+
   bool oned = topology->cell_type() == mesh::CellType::interval;
   auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
       = oned ? interval::compute_refinement_data(mesh, edges, option)

From 8a706621c4c064a6b2d66b74e1bab206ac8c8764 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Wed, 18 Sep 2024 11:49:58 +0100
Subject: [PATCH 05/31] Fix typi

---
 cpp/dolfinx/refinement/refine.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 4e72c35837c..1473a4db4e0 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -70,7 +70,7 @@ refine(const mesh::Mesh<T>& mesh,
                           mesh.geometry().cmap(), mesh.comm(),
                           std::move(new_vertex_coords), xshape, partitioner);
 
-  // Report the number of refined cellse
+  // Report the number of refined cells
   const int D = topology->dim();
   const std::int64_t n0 = topology->index_map(D)->size_global();
   const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global();

From 973f16dca8705a87c23cbbbc5d4e0401520ad139 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Sep 2024 07:44:53 +0200
Subject: [PATCH 06/31] Move maintain_coarse_partitioner to refine.cpp

---
 cpp/dolfinx/refinement/CMakeLists.txt |  1 +
 cpp/dolfinx/refinement/refine.h       | 13 ++-----------
 2 files changed, 3 insertions(+), 11 deletions(-)

diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt
index 9560e45c823..b11546361db 100644
--- a/cpp/dolfinx/refinement/CMakeLists.txt
+++ b/cpp/dolfinx/refinement/CMakeLists.txt
@@ -11,4 +11,5 @@ set(HEADERS_refinement
 target_sources(
   dolfinx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/plaza.cpp
                   ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
+                  ${CMAKE_CURRENT_SOURCE_DIR}/refine.cpp
 )
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 1473a4db4e0..e55245d4750 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -22,18 +22,9 @@ namespace dolfinx::refinement
 {
 
 /// TODO: Document function
-inline graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
+graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
     MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-    const std::vector<std::span<const std::int64_t>>& cell_topology)
-{
-  int mpi_rank = MPI::rank(comm);
-  int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
-  std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
-  std::vector<std::int32_t> destinations(num_cells, mpi_rank);
-  std::vector<std::int32_t> dest_offsets(num_cells + 1);
-  std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
-  return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
-}
+    const std::vector<std::span<const std::int64_t>>& cell_topology);
 
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.

From face987aba59853b04a65bf65459da9057a1f76c Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Sep 2024 07:47:41 +0200
Subject: [PATCH 07/31] Works better with the actual file

---
 cpp/dolfinx/refinement/refine.cpp | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 cpp/dolfinx/refinement/refine.cpp

diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp
new file mode 100644
index 00000000000..961f5cb98fc
--- /dev/null
+++ b/cpp/dolfinx/refinement/refine.cpp
@@ -0,0 +1,25 @@
+// Copyright (C) 2024 Paul T. Kühner
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include "refine.h"
+
+namespace dolfinx::refinement
+{
+
+graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
+    MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
+    const std::vector<std::span<const std::int64_t>>& cell_topology)
+{
+  int mpi_rank = MPI::rank(comm);
+  int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
+  std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
+  std::vector<std::int32_t> destinations(num_cells, mpi_rank);
+  std::vector<std::int32_t> dest_offsets(num_cells + 1);
+  std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
+  return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
+}
+
+} // namespace dolfinx::refinement

From 965c28bdcb59166da33a36a2cb51efffe044c39c Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Sep 2024 07:53:48 +0200
Subject: [PATCH 08/31] Add doc string to partitioner

---
 cpp/dolfinx/refinement/refine.h | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index e55245d4750..4d1420e7460 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -21,7 +21,15 @@
 namespace dolfinx::refinement
 {
 
-/// TODO: Document function
+/// @brief Partitioner that maintains the regional destribution for a refined
+/// mesh. Meaning, process local data is just maintained and no redistribution
+/// happens.
+///
+/// @param[in] comm MPI Communicator
+/// @param[in] nparts Number of partitions (input has no effect)
+/// @param[in] cell_types Cell types in the mesh
+/// @param[in] cells Lists of cells of each cell type.
+/// @return Destination ranks for each cell on this process
 graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
     MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
     const std::vector<std::span<const std::int64_t>>& cell_topology);

From ba75f6afd86dd897bcfcd6a7aa595d0be789ee0b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Sep 2024 08:01:37 +0200
Subject: [PATCH 09/31] Make doxygen happy

---
 cpp/dolfinx/refinement/refine.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 4d1420e7460..dddfa908104 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -31,7 +31,7 @@ namespace dolfinx::refinement
 /// @param[in] cells Lists of cells of each cell type.
 /// @return Destination ranks for each cell on this process
 graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
-    MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
+    MPI_Comm comm, int nparts, const std::vector<mesh::CellType>& cell_types,
     const std::vector<std::span<const std::int64_t>>& cell_topology);
 
 /// @brief Refine with markers, optionally redistributing, and

From 1f001cae4e81949fe6400a0e97957d8b542ace4f Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Sep 2024 08:13:11 +0200
Subject: [PATCH 10/31] Switch to general partitioner naming

---
 cpp/dolfinx/refinement/refine.cpp | 4 ++--
 cpp/dolfinx/refinement/refine.h   | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp
index 961f5cb98fc..2b5ce85d3c3 100644
--- a/cpp/dolfinx/refinement/refine.cpp
+++ b/cpp/dolfinx/refinement/refine.cpp
@@ -11,11 +11,11 @@ namespace dolfinx::refinement
 
 graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
     MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-    const std::vector<std::span<const std::int64_t>>& cell_topology)
+    const std::vector<std::span<const std::int64_t>>& cells)
 {
   int mpi_rank = MPI::rank(comm);
   int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
-  std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
+  std::int32_t num_cells = cells.front().size() / num_cell_vertices;
   std::vector<std::int32_t> destinations(num_cells, mpi_rank);
   std::vector<std::int32_t> dest_offsets(num_cells + 1);
   std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index dddfa908104..c23bc0a3c30 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -32,7 +32,7 @@ namespace dolfinx::refinement
 /// @return Destination ranks for each cell on this process
 graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
     MPI_Comm comm, int nparts, const std::vector<mesh::CellType>& cell_types,
-    const std::vector<std::span<const std::int64_t>>& cell_topology);
+    const std::vector<std::span<const std::int64_t>>& cells);
 
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.

From 22127f7248d495fd71ed485c5490e5f60b665d9b Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Mon, 23 Sep 2024 20:03:05 +0100
Subject: [PATCH 11/31] Small updates

---
 cpp/dolfinx/refinement/refine.cpp      | 15 ++++++--------
 cpp/dolfinx/refinement/refine.h        | 14 ++++++-------
 python/dolfinx/wrappers/refinement.cpp | 27 +++++++++++---------------
 3 files changed, 24 insertions(+), 32 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp
index 2b5ce85d3c3..3fba0e2410e 100644
--- a/cpp/dolfinx/refinement/refine.cpp
+++ b/cpp/dolfinx/refinement/refine.cpp
@@ -5,21 +5,18 @@
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
 #include "refine.h"
+#include <numeric>
 
-namespace dolfinx::refinement
-{
+using namespace dolfinx;
 
-graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
-    MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-    const std::vector<std::span<const std::int64_t>>& cells)
+graph::AdjacencyList<std::int32_t> refinement::maintain_coarse_partitioner(
+    MPI_Comm comm, int /*nparts*/, std::vector<mesh::CellType> cell_types,
+    std::vector<std::span<const std::int64_t>> cells)
 {
-  int mpi_rank = MPI::rank(comm);
   int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
   std::int32_t num_cells = cells.front().size() / num_cell_vertices;
-  std::vector<std::int32_t> destinations(num_cells, mpi_rank);
+  std::vector<std::int32_t> destinations(num_cells, dolfinx::MPI::rank(comm));
   std::vector<std::int32_t> dest_offsets(num_cells + 1);
   std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
   return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
 }
-
-} // namespace dolfinx::refinement
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index c23bc0a3c30..07a5b411a86 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -20,19 +20,19 @@
 
 namespace dolfinx::refinement
 {
-
-/// @brief Partitioner that maintains the regional destribution for a refined
-/// mesh. Meaning, process local data is just maintained and no redistribution
-/// happens.
+/// @brief Partitioner that maintains the regional distribution for a
+/// refined mesh. Meaning, process local data is just maintained and no
+/// redistribution happens.
 ///
 /// @param[in] comm MPI Communicator
 /// @param[in] nparts Number of partitions (input has no effect)
 /// @param[in] cell_types Cell types in the mesh
 /// @param[in] cells Lists of cells of each cell type.
 /// @return Destination ranks for each cell on this process
-graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
-    MPI_Comm comm, int nparts, const std::vector<mesh::CellType>& cell_types,
-    const std::vector<std::span<const std::int64_t>>& cells);
+graph::AdjacencyList<std::int32_t>
+maintain_coarse_partitioner(MPI_Comm comm, int nparts,
+                            std::vector<mesh::CellType> cell_types,
+                            std::vector<std::span<const std::int64_t>> cells);
 
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 7921ffb22e5..8bb89b7260d 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -4,12 +4,18 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
+#include "MPICommWrapper.h"
+#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 <functional>
-#include <optional>
-#include <span>
-
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/array.h>
@@ -17,26 +23,15 @@
 #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 "MPICommWrapper.h"
-#include "array.h"
+#include <optional>
+#include <span>
 
 namespace nb = nanobind;
 
 namespace dolfinx_wrappers
 {
-
 namespace
 {
-
 using PythonCellPartitionFunction
     = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
         dolfinx_wrappers::MPICommWrapper, int,

From 3b74e494b7fd1bf2bbd6e4f1ca9d01a8da4aed2b Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Fri, 27 Sep 2024 09:38:33 +0100
Subject: [PATCH 12/31] Fix test failure

---
 python/test/unit/refinement/test_refinement.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 25b1d460bbc..202fed23393 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -211,5 +211,5 @@ def test_refine_cell_meshtag(tdim, refine_plaza_wrapper):
 def test_refine_ufl_cargo():
     mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3)
     mesh.topology.create_entities(1)
-    refined_mesh, _, _ = refine(mesh, redistribute=False)
+    refined_mesh, _, _ = refine(mesh)
     assert refined_mesh.ufl_domain().ufl_cargo() != mesh.ufl_domain().ufl_cargo()

From 2e8a65ca8108b34a494aee76f47a694b703d61b3 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 10:14:02 +0100
Subject: [PATCH 13/31] Simplifications

---
 cpp/dolfinx/refinement/CMakeLists.txt |  1 -
 cpp/dolfinx/refinement/refine.cpp     | 22 --------------------
 cpp/dolfinx/refinement/refine.h       | 30 +++++++--------------------
 3 files changed, 8 insertions(+), 45 deletions(-)
 delete mode 100644 cpp/dolfinx/refinement/refine.cpp

diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt
index b11546361db..9560e45c823 100644
--- a/cpp/dolfinx/refinement/CMakeLists.txt
+++ b/cpp/dolfinx/refinement/CMakeLists.txt
@@ -11,5 +11,4 @@ set(HEADERS_refinement
 target_sources(
   dolfinx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/plaza.cpp
                   ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
-                  ${CMAKE_CURRENT_SOURCE_DIR}/refine.cpp
 )
diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp
deleted file mode 100644
index 3fba0e2410e..00000000000
--- a/cpp/dolfinx/refinement/refine.cpp
+++ /dev/null
@@ -1,22 +0,0 @@
-// Copyright (C) 2024 Paul T. Kühner
-//
-// This file is part of DOLFINx (https://www.fenicsproject.org)
-//
-// SPDX-License-Identifier:    LGPL-3.0-or-later
-
-#include "refine.h"
-#include <numeric>
-
-using namespace dolfinx;
-
-graph::AdjacencyList<std::int32_t> refinement::maintain_coarse_partitioner(
-    MPI_Comm comm, int /*nparts*/, std::vector<mesh::CellType> cell_types,
-    std::vector<std::span<const std::int64_t>> cells)
-{
-  int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
-  std::int32_t num_cells = cells.front().size() / num_cell_vertices;
-  std::vector<std::int32_t> destinations(num_cells, dolfinx::MPI::rank(comm));
-  std::vector<std::int32_t> dest_offsets(num_cells + 1);
-  std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
-  return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
-}
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 07a5b411a86..17346a80180 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -20,38 +20,24 @@
 
 namespace dolfinx::refinement
 {
-/// @brief Partitioner that maintains the regional distribution for a
-/// refined mesh. Meaning, process local data is just maintained and no
-/// redistribution happens.
-///
-/// @param[in] comm MPI Communicator
-/// @param[in] nparts Number of partitions (input has no effect)
-/// @param[in] cell_types Cell types in the mesh
-/// @param[in] cells Lists of cells of each cell type.
-/// @return Destination ranks for each cell on this process
-graph::AdjacencyList<std::int32_t>
-maintain_coarse_partitioner(MPI_Comm comm, int nparts,
-                            std::vector<mesh::CellType> cell_types,
-                            std::vector<std::span<const std::int64_t>> cells);
-
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
 ///
 /// @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] edges Indices of the edges that should be split by this
+/// refinement. If not provides, a uniform refinement is performed.
 /// @param[in] partitioner partitioner to be used for 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
+/// @return New mesh and optional parent cell index, parent facet
 /// indices.
 template <std::floating_point T>
 std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
            std::optional<std::vector<std::int8_t>>>
 refine(const mesh::Mesh<T>& mesh,
        std::optional<std::span<const std::int32_t>> edges,
-       mesh::CellPartitionFunction partitioner = maintain_coarse_partitioner,
+       const mesh::CellPartitionFunction& partitioner
+       = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet),
        Option option = Option::none)
 {
   auto topology = mesh.topology();
@@ -59,10 +45,10 @@ refine(const mesh::Mesh<T>& mesh,
   if (!mesh::is_simplex(topology->cell_type()))
     throw std::runtime_error("Refinement only defined for simplices");
 
-  bool oned = topology->cell_type() == mesh::CellType::interval;
   auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
-      = oned ? interval::compute_refinement_data(mesh, edges, option)
-             : plaza::compute_refinement_data(mesh, edges, option);
+      = (topology->cell_type() == mesh::CellType::interval)
+            ? interval::compute_refinement_data(mesh, edges, option)
+            : plaza::compute_refinement_data(mesh, edges, option);
 
   mesh::Mesh<T> refined_mesh
       = mesh::create_mesh(mesh.comm(), mesh.comm(), std::move(cell_adj).array(),

From 905497f1073eacf6affaad8a969749debc78b13c Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 10:16:56 +0100
Subject: [PATCH 14/31] Small fix

---
 cpp/test/mesh/refinement/interval.cpp | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp
index 74b3f434627..a1ed0d7b9ab 100644
--- a/cpp/test/mesh/refinement/interval.cpp
+++ b/cpp/test/mesh/refinement/interval.cpp
@@ -73,8 +73,7 @@ TEMPLATE_TEST_CASE("Interval uniform refinement",
 
   // TODO: parent_facet
   auto [refined_mesh, parent_edge, parent_facet] = refinement::refine(
-      mesh, std::nullopt, &refinement::maintain_coarse_partitioner,
-      refinement::Option::parent_cell);
+      mesh, std::nullopt, nullptr, refinement::Option::parent_cell);
 
   std::vector<T> expected_x = {
       /* v_0 */ 0.0, 0.0, 0.0,
@@ -193,8 +192,7 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)",
 
   // TODO: parent_facet
   auto [refined_mesh, parent_edges, parent_facet] = refinement::refine(
-      mesh, std::nullopt, &refinement::maintain_coarse_partitioner,
-      refinement::Option::parent_cell);
+      mesh, std::nullopt, nullptr, refinement::Option::parent_cell);
 
   T rank_d = static_cast<T>(rank);
   T comm_size_d = static_cast<T>(comm_size);

From d701654385c4663ce2b588a220a66aad0fb9fdd2 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 10:23:38 +0100
Subject: [PATCH 15/31] Small fix

---
 python/dolfinx/wrappers/refinement.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 8bb89b7260d..df39d4fb241 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -93,7 +93,7 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
         CppCellPartitionFunction cpp_partitioner
             = partitioner.has_value()
                   ? create_cell_partitioner_cpp(partitioner.value())
-                  : dolfinx::refinement::maintain_coarse_partitioner;
+                  : nullptr;
         auto [mesh1, cell, facet] = dolfinx::refinement::refine(
             mesh, cpp_edges, cpp_partitioner, option);
 

From 9a651e41559070300411004a9671d42261b0fe51 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 11:03:37 +0100
Subject: [PATCH 16/31] Work on docs

---
 cpp/dolfinx/refinement/refine.h | 13 ++++++++++++-
 python/dolfinx/mesh.py          | 20 +++++++++++++++-----
 2 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 17346a80180..5bd3a80703e 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -23,10 +23,21 @@ namespace dolfinx::refinement
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
 ///
+/// @note Using the default partitioner for a refined mesh, the refined
+/// mesh will include ghosts cells (connected by facet).
+///
+/// @warning Passing `nullptr` for `partitioner`, refined cells will be
+/// on the same process as the parent cell but the refined mesh will
+/// **not** have ghosts cells. The possibility to not re-partition the
+/// refined mesh and include ghost cells in the refined mesh will be
+/// added in a future release.
+///
 /// @param[in] mesh Input mesh to be refined.
 /// @param[in] edges Indices of the edges that should be split by this
 /// refinement. If not provides, a uniform refinement is performed.
-/// @param[in] partitioner partitioner to be used for the refined mesh
+/// @param[in] partitioner Partitioner to be used for the refined mesh.
+/// If not callable, refined cells will be on the same process as the
+/// parent cell.
 /// @param[in] option Control the computation of parent facets, parent
 /// cells. If an option is unselected, an empty list is returned.
 /// @return New mesh and optional parent cell index, parent facet
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index bac5232cbb1..ccb0f272835 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -535,19 +535,29 @@ def transfer_meshtag(
 def refine(
     mesh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
-    partitioner: typing.Optional[typing.Callable] = None,
+    partitioner: typing.Optional[typing.Callable] = _cpp.mesh.create_cell_partitioner(
+        GhostMode.shared_facet
+    ),
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
 
+    Warning:
+        Passing ``None`` for ``partitioner``, refined cells will be on
+        the same process as the parent cell but the refined mesh will
+        **not** have ghosts cells. The possibility to not re-partition
+        the refined mesh and include ghost cells in the refined will be
+        added in a future release.
+
     Args:
         mesh: Mesh from which to create the refined mesh.
         edges: Indices of edges to split during refinement. If ``None``,
             mesh refinement is uniform.
-        partitioner:
-            partitioner to use for the refined mesh, If ``None`` no redistribution is performed,
-            i.e. previous local mesh is equally parallelized now with new vertices.s
-        option: Controls whether parent cells and/or parent facets are computed.
+        partitioner: Partitioner to distribute the refined mesh. If
+            ``None`` no redistribution is performed, i.e. refined cells remain on the
+            same process as the parent cell.
+        option: Controls whether parent cells and/or parent facets are
+            computed.
 
     Returns:
        Refined mesh, (optional) parent cells, (optional) parent facets

From f52c12feb2e57c468e153bdb5444a1fa4988d570 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 11:24:52 +0100
Subject: [PATCH 17/31] Simplify

---
 python/dolfinx/wrappers/mesh.cpp       | 111 +++++--------------------
 python/dolfinx/wrappers/mesh.h         |  81 ++++++++++++++++++
 python/dolfinx/wrappers/refinement.cpp |  47 +----------
 3 files changed, 108 insertions(+), 131 deletions(-)
 create mode 100644 python/dolfinx/wrappers/mesh.h

diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index f8c3679974d..913ca654592 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -4,11 +4,11 @@
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
+#include "mesh.h"
 #include "MPICommWrapper.h"
 #include "array.h"
 #include "caster_mpi.h"
 #include "numpy_dtype.h"
-#include <cfloat>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/fem/CoordinateElement.h>
 #include <dolfinx/fem/ElementDofLayout.h>
@@ -21,7 +21,6 @@
 #include <dolfinx/mesh/graphbuild.h>
 #include <dolfinx/mesh/topologycomputation.h>
 #include <dolfinx/mesh/utils.h>
-#include <iostream>
 #include <memory>
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
@@ -36,75 +35,6 @@
 
 namespace nb = nanobind;
 
-namespace
-{
-/// Wrap a Python graph partitioning function as a C++ function
-template <typename Functor>
-auto create_partitioner_cpp(Functor p)
-{
-  return [p](MPI_Comm comm, int nparts,
-             const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
-             bool ghosting)
-  {
-    return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph,
-             ghosting);
-  };
-}
-
-/// Wrap a C++ cell partitioning function as a Python function
-template <typename Functor>
-auto create_cell_partitioner_py(Functor p)
-{
-  return [p](dolfinx_wrappers::MPICommWrapper comm, int n,
-             const std::vector<dolfinx::mesh::CellType>& cell_types,
-             std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb)
-  {
-    std::vector<std::span<const std::int64_t>> cells;
-    std::ranges::transform(
-        cells_nb, std::back_inserter(cells), [](auto c)
-        { return std::span<const std::int64_t>(c.data(), c.size()); });
-    return p(comm.get(), n, cell_types, cells);
-  };
-}
-
-using PythonCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        dolfinx_wrappers::MPICommWrapper, int,
-        const std::vector<dolfinx::mesh::CellType>&,
-        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
-
-using CppCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
-        const std::vector<std::span<const std::int64_t>>&)>;
-
-/// Wrap a Python cell graph partitioning function as a C++ function
-CppCellPartitionFunction
-create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
-{
-  if (p)
-  {
-    return [p](MPI_Comm comm, int n,
-               const std::vector<dolfinx::mesh::CellType>& cell_types,
-               const std::vector<std::span<const std::int64_t>>& cells)
-    {
-      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
-      std::ranges::transform(
-          cells, std::back_inserter(cells_nb),
-          [](auto c)
-          {
-            return nb::ndarray<const std::int64_t, nb::numpy>(
-                c.data(), {c.size()}, nb::handle());
-          });
-
-      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
-    };
-  }
-  else
-    return nullptr;
-}
-} // namespace
-
 namespace dolfinx_wrappers
 {
 template <typename T>
@@ -260,7 +190,8 @@ void declare_mesh(nb::module_& m, std::string type)
           "__init__",
           [](dolfinx::mesh::Mesh<T>* mesh, MPICommWrapper comm,
              std::shared_ptr<dolfinx::mesh::Topology> topology,
-             dolfinx::mesh::Geometry<T>& geometry) {
+             dolfinx::mesh::Geometry<T>& geometry)
+          {
             new (mesh) dolfinx::mesh::Mesh<T>(comm.get(), topology, geometry);
           },
           nb::arg("comm"), nb::arg("topology"), nb::arg("geometry"))
@@ -280,10 +211,11 @@ void declare_mesh(nb::module_& m, std::string type)
       create_interval.c_str(),
       [](MPICommWrapper comm, std::int64_t n, std::array<T, 2> p,
          dolfinx::mesh::GhostMode ghost_mode,
-         const PythonCellPartitionFunction& part)
+         const part::impl::PythonCellPartitionFunction& part)
       {
         return dolfinx::mesh::create_interval<T>(
-            comm.get(), n, p, ghost_mode, create_cell_partitioner_cpp(part));
+            comm.get(), n, p, ghost_mode,
+            part::impl::create_cell_partitioner_cpp(part));
       },
       nb::arg("comm"), nb::arg("n"), nb::arg("p"), nb::arg("ghost_mode"),
       nb::arg("partitioner").none());
@@ -293,12 +225,12 @@ void declare_mesh(nb::module_& m, std::string type)
       create_rectangle.c_str(),
       [](MPICommWrapper comm, std::array<std::array<T, 2>, 2> p,
          std::array<std::int64_t, 2> n, dolfinx::mesh::CellType celltype,
-         const PythonCellPartitionFunction& part,
+         const part::impl::PythonCellPartitionFunction& part,
          dolfinx::mesh::DiagonalType diagonal)
       {
         return dolfinx::mesh::create_rectangle<T>(
-            comm.get(), p, n, celltype, create_cell_partitioner_cpp(part),
-            diagonal);
+            comm.get(), p, n, celltype,
+            part::impl::create_cell_partitioner_cpp(part), diagonal);
       },
       nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"),
       nb::arg("partitioner").none(), nb::arg("diagonal"));
@@ -308,11 +240,12 @@ void declare_mesh(nb::module_& m, std::string type)
       create_box.c_str(),
       [](MPICommWrapper comm, std::array<std::array<T, 3>, 2> p,
          std::array<std::int64_t, 3> n, dolfinx::mesh::CellType celltype,
-         const PythonCellPartitionFunction& part)
+         const part::impl::PythonCellPartitionFunction& part)
       {
         MPI_Comm _comm = comm.get();
-        return dolfinx::mesh::create_box<T>(_comm, _comm, p, n, celltype,
-                                            create_cell_partitioner_cpp(part));
+        return dolfinx::mesh::create_box<T>(
+            _comm, _comm, p, n, celltype,
+            part::impl::create_cell_partitioner_cpp(part));
       },
       nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"),
       nb::arg("partitioner").none());
@@ -323,7 +256,7 @@ void declare_mesh(nb::module_& m, std::string type)
                                          nb::c_contig>>& cells_nb,
            const std::vector<dolfinx::fem::CoordinateElement<T>>& elements,
            nb::ndarray<const T, nb::c_contig> x,
-           const PythonCellPartitionFunction& p)
+           const part::impl::PythonCellPartitionFunction& p)
         {
           std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1);
 
@@ -365,7 +298,7 @@ void declare_mesh(nb::module_& m, std::string type)
          nb::ndarray<const std::int64_t, nb::ndim<2>, nb::c_contig> cells,
          const dolfinx::fem::CoordinateElement<T>& element,
          nb::ndarray<const T, nb::c_contig> x,
-         const PythonCellPartitionFunction& p)
+         const part::impl::PythonCellPartitionFunction& p)
       {
         std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1);
         if (p)
@@ -599,7 +532,8 @@ void mesh(nb::module_& m)
   m.def(
       "compute_entities",
       [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology, int dim,
-         int index) {
+         int index)
+      {
         return dolfinx::mesh::compute_entities(comm.get(), topology, dim,
                                                index);
       },
@@ -733,9 +667,9 @@ void mesh(nb::module_& m)
 
   m.def(
       "create_cell_partitioner",
-      [](dolfinx::mesh::GhostMode gm) -> PythonCellPartitionFunction
+      [](dolfinx::mesh::GhostMode gm) -> part::impl::PythonCellPartitionFunction
       {
-        return create_cell_partitioner_py(
+        return part::impl::create_cell_partitioner_py(
             dolfinx::mesh::create_cell_partitioner(gm));
       },
       "Create default cell partitioner.");
@@ -746,11 +680,12 @@ void mesh(nb::module_& m)
              const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
              bool ghosting)>
              part,
-         dolfinx::mesh::GhostMode ghost_mode) -> PythonCellPartitionFunction
+         dolfinx::mesh::GhostMode ghost_mode)
+          -> part::impl::PythonCellPartitionFunction
       {
-        return create_cell_partitioner_py(
+        return part::impl::create_cell_partitioner_py(
             dolfinx::mesh::create_cell_partitioner(
-                ghost_mode, create_partitioner_cpp(part)));
+                ghost_mode, part::impl::create_partitioner_cpp(part)));
       },
       nb::arg("part"), nb::arg("ghost_mode") = dolfinx::mesh::GhostMode::none,
       "Create a cell partitioner from a graph partitioning function.");
diff --git a/python/dolfinx/wrappers/mesh.h b/python/dolfinx/wrappers/mesh.h
new file mode 100644
index 00000000000..cc8df9d82e9
--- /dev/null
+++ b/python/dolfinx/wrappers/mesh.h
@@ -0,0 +1,81 @@
+// Copyright (C) 2017-2024 Chris N. Richardson and Garth N. Wells
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include "MPICommWrapper.h"
+#include <dolfinx/mesh/cell_types.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/vector.h>
+
+namespace nb = nanobind;
+
+namespace dolfinx_wrappers::part::impl
+{
+/// Wrap a Python graph partitioning function as a C++ function
+template <typename Functor>
+auto create_partitioner_cpp(Functor p)
+{
+  return [p](MPI_Comm comm, int nparts,
+             const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
+             bool ghosting)
+  {
+    return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph,
+             ghosting);
+  };
+}
+
+/// Wrap a C++ cell partitioning function as a Python function
+template <typename Functor>
+auto create_cell_partitioner_py(Functor p)
+{
+  return [p](dolfinx_wrappers::MPICommWrapper comm, int n,
+             const std::vector<dolfinx::mesh::CellType>& cell_types,
+             std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb)
+  {
+    std::vector<std::span<const std::int64_t>> cells;
+    std::ranges::transform(
+        cells_nb, std::back_inserter(cells), [](auto c)
+        { return std::span<const std::int64_t>(c.data(), c.size()); });
+    return p(comm.get(), n, cell_types, cells);
+  };
+}
+
+using PythonCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        dolfinx_wrappers::MPICommWrapper, int,
+        const std::vector<dolfinx::mesh::CellType>&,
+        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
+
+using CppCellPartitionFunction
+    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
+        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
+        const std::vector<std::span<const std::int64_t>>&)>;
+
+/// Wrap a Python cell graph partitioning function as a C++ function
+inline CppCellPartitionFunction
+create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
+{
+  if (p)
+  {
+    return [p](MPI_Comm comm, int n,
+               const std::vector<dolfinx::mesh::CellType>& cell_types,
+               const std::vector<std::span<const std::int64_t>>& cells)
+    {
+      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
+      std::ranges::transform(
+          cells, std::back_inserter(cells_nb),
+          [](auto c)
+          {
+            return nb::ndarray<const std::int64_t, nb::numpy>(
+                c.data(), {c.size()}, nb::handle());
+          });
+
+      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
+    };
+  }
+  else
+    return nullptr;
+}
+} // namespace dolfinx_wrappers::part::impl
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index df39d4fb241..399f5d2305c 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -6,6 +6,7 @@
 
 #include "MPICommWrapper.h"
 #include "array.h"
+#include "mesh.h"
 #include <concepts>
 #include <cstdint>
 #include <dolfinx/mesh/Mesh.h>
@@ -30,46 +31,6 @@ namespace nb = nanobind;
 
 namespace dolfinx_wrappers
 {
-namespace
-{
-using PythonCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        dolfinx_wrappers::MPICommWrapper, int,
-        const std::vector<dolfinx::mesh::CellType>&,
-        std::vector<nb::ndarray<const std::int64_t, nb::numpy>>)>;
-
-using CppCellPartitionFunction
-    = std::function<dolfinx::graph::AdjacencyList<std::int32_t>(
-        MPI_Comm, int, const std::vector<dolfinx::mesh::CellType>& q,
-        const std::vector<std::span<const std::int64_t>>&)>;
-
-/// Wrap a Python graph partitioning function as a C++ function
-CppCellPartitionFunction
-create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
-{
-  if (p)
-  {
-    return [p](MPI_Comm comm, int n,
-               const std::vector<dolfinx::mesh::CellType>& cell_types,
-               const std::vector<std::span<const std::int64_t>>& cells)
-    {
-      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
-      std::ranges::transform(
-          cells, std::back_inserter(cells_nb),
-          [](auto c)
-          {
-            return nb::ndarray<const std::int64_t, nb::numpy>(
-                c.data(), {c.size()}, nb::handle());
-          });
-
-      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
-    };
-  }
-  else
-    return nullptr;
-}
-} // namespace
-
 template <std::floating_point T>
 void export_refinement_with_variable_mesh_type(nb::module_& m)
 {
@@ -80,7 +41,7 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
          std::optional<
              nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
              edges,
-         std::optional<PythonCellPartitionFunction> partitioner,
+         std::optional<part::impl::PythonCellPartitionFunction> partitioner,
          dolfinx::refinement::Option option)
       {
         std::optional<std::span<const std::int32_t>> cpp_edges(std::nullopt);
@@ -90,9 +51,9 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
               std::span(edges.value().data(), edges.value().size()));
         }
 
-        CppCellPartitionFunction cpp_partitioner
+        part::impl::CppCellPartitionFunction cpp_partitioner
             = partitioner.has_value()
-                  ? create_cell_partitioner_cpp(partitioner.value())
+                  ? part::impl::create_cell_partitioner_cpp(partitioner.value())
                   : nullptr;
         auto [mesh1, cell, facet] = dolfinx::refinement::refine(
             mesh, cpp_edges, cpp_partitioner, option);

From fb8edd1d5b598f2384f7079868e19ae8e4feda12 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 11:27:48 +0100
Subject: [PATCH 18/31] Improve function name

---
 python/dolfinx/wrappers/refinement.cpp | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 399f5d2305c..c644cb76d53 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -32,9 +32,8 @@ namespace nb = nanobind;
 namespace dolfinx_wrappers
 {
 template <std::floating_point T>
-void export_refinement_with_variable_mesh_type(nb::module_& m)
+void export_refinement(nb::module_& m)
 {
-
   m.def(
       "refine",
       [](const dolfinx::mesh::Mesh<T>& mesh,
@@ -77,8 +76,8 @@ void export_refinement_with_variable_mesh_type(nb::module_& m)
 
 void refinement(nb::module_& m)
 {
-  export_refinement_with_variable_mesh_type<float>(m);
-  export_refinement_with_variable_mesh_type<double>(m);
+  export_refinement<float>(m);
+  export_refinement<double>(m);
 
   nb::enum_<dolfinx::refinement::Option>(m, "RefinementOption")
       .value("none", dolfinx::refinement::Option::none)

From 6989537bbd9891b252db402ab0a7d828aed4306b Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 11:28:21 +0100
Subject: [PATCH 19/31] Namespace fix

---
 python/dolfinx/wrappers/refinement.cpp | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index c644cb76d53..7c4ac5dc98e 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -29,7 +29,7 @@
 
 namespace nb = nanobind;
 
-namespace dolfinx_wrappers
+namespace
 {
 template <std::floating_point T>
 void export_refinement(nb::module_& m)
@@ -73,6 +73,10 @@ void export_refinement(nb::module_& m)
       nb::arg("mesh"), nb::arg("edges") = nb::none(),
       nb::arg("partitioner") = nb::none(), nb::arg("option"));
 }
+} // namespace
+
+namespace dolfinx_wrappers
+{
 
 void refinement(nb::module_& m)
 {

From c45574e68c4219b61f42ef993be8c2e8baa17e87 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 11:30:40 +0100
Subject: [PATCH 20/31] Compile fix

---
 python/dolfinx/wrappers/refinement.cpp | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 7c4ac5dc98e..756c008c102 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -40,7 +40,9 @@ void export_refinement(nb::module_& m)
          std::optional<
              nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
              edges,
-         std::optional<part::impl::PythonCellPartitionFunction> partitioner,
+         std::optional<
+             dolfinx_wrappers::part::impl::PythonCellPartitionFunction>
+             partitioner,
          dolfinx::refinement::Option option)
       {
         std::optional<std::span<const std::int32_t>> cpp_edges(std::nullopt);
@@ -50,9 +52,10 @@ void export_refinement(nb::module_& m)
               std::span(edges.value().data(), edges.value().size()));
         }
 
-        part::impl::CppCellPartitionFunction cpp_partitioner
+        dolfinx_wrappers::part::impl::CppCellPartitionFunction cpp_partitioner
             = partitioner.has_value()
-                  ? part::impl::create_cell_partitioner_cpp(partitioner.value())
+                  ? dolfinx_wrappers::part::impl::create_cell_partitioner_cpp(
+                        partitioner.value())
                   : nullptr;
         auto [mesh1, cell, facet] = dolfinx::refinement::refine(
             mesh, cpp_edges, cpp_partitioner, option);
@@ -60,12 +63,18 @@ void export_refinement(nb::module_& m)
         std::optional<nb::ndarray<std::int32_t, nb::numpy>> python_cell(
             std::nullopt);
         if (cell.has_value())
-          python_cell.emplace(as_nbarray(std::move(cell.value())));
+        {
+          python_cell.emplace(
+              dolfinx_wrappers::as_nbarray(std::move(cell.value())));
+        }
 
         std::optional<nb::ndarray<std::int8_t, nb::numpy>> python_facet(
             std::nullopt);
         if (facet.has_value())
-          python_facet.emplace(as_nbarray(std::move(facet.value())));
+        {
+          python_facet.emplace(
+              dolfinx_wrappers::as_nbarray(std::move(facet.value())));
+        }
 
         return std::tuple{std::move(mesh1), std::move(python_cell),
                           std::move(python_facet)};

From 35df167971fa9147dd26a710cfb6e3f5a63d9e1d Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 12:19:51 +0100
Subject: [PATCH 21/31] Add nanobind include

---
 cpp/dolfinx/refinement/refine.h        |  3 ++-
 python/dolfinx/mesh.py                 | 18 +++++++++++++-----
 python/dolfinx/wrappers/refinement.cpp |  9 ++++-----
 3 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 5bd3a80703e..b625bdfc510 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -24,7 +24,8 @@ namespace dolfinx::refinement
 /// optionally calculating the parent-child relationships.
 ///
 /// @note Using the default partitioner for a refined mesh, the refined
-/// mesh will include ghosts cells (connected by facet).
+/// mesh will include ghosts cells (connected by facet) even if the
+/// parent mesh is not ghosted.
 ///
 /// @warning Passing `nullptr` for `partitioner`, refined cells will be
 /// on the same process as the parent cell but the refined mesh will
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index ccb0f272835..b073d068e7b 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -542,12 +542,20 @@ def refine(
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
 
+    Passing ``None`` for ``partitioner``, refined cells will be on the
+    same process as the parent cell.
+
+    Note:
+        Using the default partitioner for a refined mesh, the refined
+        mesh will include ghosts cells (connected by facet) even if the
+        parent mesh does not include ghost cells.
+
     Warning:
-        Passing ``None`` for ``partitioner``, refined cells will be on
-        the same process as the parent cell but the refined mesh will
-        **not** have ghosts cells. The possibility to not re-partition
-        the refined mesh and include ghost cells in the refined will be
-        added in a future release.
+        Passing ``None`` for ``partitioner``, the refined mesh will
+        **not** have ghosts cells even if the parent mesh has ghost
+        cells. The possibility to not re-partition the refined mesh and
+        include ghost cells in the refined mesh will be added in a
+        future release.
 
     Args:
         mesh: Mesh from which to create the refined mesh.
diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp
index 756c008c102..d9b86fcf781 100644
--- a/python/dolfinx/wrappers/refinement.cpp
+++ b/python/dolfinx/wrappers/refinement.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2018 Chris N. Richardson and Garth N. Wells
+// Copyright (C) 2018-2024 Chris N. Richardson and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -13,13 +13,12 @@
 #include <dolfinx/mesh/MeshTags.h>
 #include <dolfinx/refinement/interval.h>
 #include <dolfinx/refinement/option.h>
-#include <dolfinx/refinement/plaza.h>
 #include <dolfinx/refinement/refine.h>
 #include <dolfinx/refinement/utils.h>
 #include <functional>
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
-#include <nanobind/stl/array.h>
+#include <nanobind/stl/function.h>
 #include <nanobind/stl/optional.h>
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/tuple.h>
@@ -79,8 +78,8 @@ void export_refinement(nb::module_& m)
         return std::tuple{std::move(mesh1), std::move(python_cell),
                           std::move(python_facet)};
       },
-      nb::arg("mesh"), nb::arg("edges") = nb::none(),
-      nb::arg("partitioner") = nb::none(), nb::arg("option"));
+      nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(),
+      nb::arg("option"));
 }
 } // namespace
 

From c1e543035a5f16a475b421160d59214f590e387a Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 12:23:06 +0100
Subject: [PATCH 22/31] Small re-factor

---
 python/dolfinx/wrappers/mesh.cpp | 28 ++++++++++++++++++++++++++++
 python/dolfinx/wrappers/mesh.h   | 26 ++------------------------
 2 files changed, 30 insertions(+), 24 deletions(-)

diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index 913ca654592..e2cae585563 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -37,6 +37,34 @@ namespace nb = nanobind;
 
 namespace dolfinx_wrappers
 {
+namespace part::impl
+{
+CppCellPartitionFunction
+create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
+{
+  if (p)
+  {
+    return [p](MPI_Comm comm, int n,
+               const std::vector<dolfinx::mesh::CellType>& cell_types,
+               const std::vector<std::span<const std::int64_t>>& cells)
+    {
+      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
+      std::ranges::transform(
+          cells, std::back_inserter(cells_nb),
+          [](auto c)
+          {
+            return nb::ndarray<const std::int64_t, nb::numpy>(
+                c.data(), {c.size()}, nb::handle());
+          });
+
+      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
+    };
+  }
+  else
+    return nullptr;
+}
+} // namespace part::impl
+
 template <typename T>
 void declare_meshtags(nb::module_& m, std::string type)
 {
diff --git a/python/dolfinx/wrappers/mesh.h b/python/dolfinx/wrappers/mesh.h
index cc8df9d82e9..5471dc1ab24 100644
--- a/python/dolfinx/wrappers/mesh.h
+++ b/python/dolfinx/wrappers/mesh.h
@@ -54,28 +54,6 @@ using CppCellPartitionFunction
         const std::vector<std::span<const std::int64_t>>&)>;
 
 /// Wrap a Python cell graph partitioning function as a C++ function
-inline CppCellPartitionFunction
-create_cell_partitioner_cpp(const PythonCellPartitionFunction& p)
-{
-  if (p)
-  {
-    return [p](MPI_Comm comm, int n,
-               const std::vector<dolfinx::mesh::CellType>& cell_types,
-               const std::vector<std::span<const std::int64_t>>& cells)
-    {
-      std::vector<nb::ndarray<const std::int64_t, nb::numpy>> cells_nb;
-      std::ranges::transform(
-          cells, std::back_inserter(cells_nb),
-          [](auto c)
-          {
-            return nb::ndarray<const std::int64_t, nb::numpy>(
-                c.data(), {c.size()}, nb::handle());
-          });
-
-      return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb);
-    };
-  }
-  else
-    return nullptr;
-}
+CppCellPartitionFunction
+create_cell_partitioner_cpp(const PythonCellPartitionFunction& p);
 } // namespace dolfinx_wrappers::part::impl

From 198c4d9d9fae94114b44d3efba753be02689594f Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 18:40:54 +0100
Subject: [PATCH 23/31] Fixes

---
 cpp/dolfinx/mesh/MeshTags.h                    |  4 ++--
 cpp/dolfinx/refinement/refine.h                |  9 ++++-----
 cpp/dolfinx/refinement/utils.h                 | 12 ++++++++----
 python/dolfinx/mesh.py                         |  2 +-
 python/test/unit/refinement/test_refinement.py | 11 ++++++++---
 5 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/cpp/dolfinx/mesh/MeshTags.h b/cpp/dolfinx/mesh/MeshTags.h
index 32cb2de6a31..3026c258158 100644
--- a/cpp/dolfinx/mesh/MeshTags.h
+++ b/cpp/dolfinx/mesh/MeshTags.h
@@ -127,8 +127,8 @@ class MeshTags
 };
 
 /// @brief Create MeshTags from arrays
-/// @param[in] topology Mesh topology that the tags are associated with
-/// @param[in] dim Topological dimension of tagged entities
+/// @param[in] topology Mesh topology that the tags are associated with.
+/// @param[in] dim Topological dimension of tagged entities.
 /// @param[in] entities Local vertex indices for tagged entities.
 /// @param[in] values Tag values for each entity in `entities`. The
 /// length of `values` must be equal to number of rows in `entities`.
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index b625bdfc510..d92ba24c4ff 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -49,7 +49,7 @@ std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
 refine(const mesh::Mesh<T>& mesh,
        std::optional<std::span<const std::int32_t>> edges,
        const mesh::CellPartitionFunction& partitioner
-       = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet),
+       = mesh::create_cell_partitioner(mesh::GhostMode::none),
        Option option = Option::none)
 {
   auto topology = mesh.topology();
@@ -62,10 +62,9 @@ refine(const mesh::Mesh<T>& mesh,
             ? interval::compute_refinement_data(mesh, edges, option)
             : plaza::compute_refinement_data(mesh, edges, option);
 
-  mesh::Mesh<T> refined_mesh
-      = mesh::create_mesh(mesh.comm(), mesh.comm(), std::move(cell_adj).array(),
-                          mesh.geometry().cmap(), mesh.comm(),
-                          std::move(new_vertex_coords), xshape, partitioner);
+  mesh::Mesh<T> refined_mesh = mesh::create_mesh(
+      mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(),
+      mesh.comm(), new_vertex_coords, xshape, partitioner);
 
   // Report the number of refined cells
   const int D = topology->dim();
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index 1c6af077551..d0f6d158e4b 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -276,9 +276,12 @@ std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
                                          std::int32_t n);
 
 /// @brief Transfer facet MeshTags from coarse mesh to refined mesh
-/// @note The refined mesh must not have been redistributed during
-/// refinement
-/// @note GhostMode must be GhostMode.none
+///
+/// @warning The refined mesh must not have been redistributed during
+/// refinement.
+///
+/// @warning GhostMode must be GhostMode.none.
+///
 /// @param[in] tags0 Tags on the parent mesh
 /// @param[in] topology1 Refined mesh topology
 /// @param[in] cell Parent cell of each cell in refined mesh
@@ -292,7 +295,8 @@ std::array<std::vector<std::int32_t>, 2> transfer_facet_meshtag(
 ///
 /// @note The refined mesh must not have been redistributed during
 /// refinement.
-/// @note GhostMode must be GhostMode.none
+///
+/// @warning GhostMode must be GhostMode.none.
 ///
 /// @param[in] tags0 Tags on the parent mesh
 /// @param[in] topology1 Refined mesh topology
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index b073d068e7b..4cc7d001859 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -536,7 +536,7 @@ def refine(
     mesh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
     partitioner: typing.Optional[typing.Callable] = _cpp.mesh.create_cell_partitioner(
-        GhostMode.shared_facet
+        GhostMode.none
     ),
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 202fed23393..719acc8e366 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -112,14 +112,18 @@ def left_side(x, tol=1e-16):
     assert num_cells_global == actual_cells
 
 
-@pytest.mark.parametrize("tdim", [2, 3])
+@pytest.mark.parametrize(
+    "tdim",
+    [2, 3],
+)
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet),
+        lambda mesh: refine(mesh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
         lambda mesh: refine(
             mesh,
-            np.arange(mesh.topology.index_map(1).size_local),
+            edges=np.arange(mesh.topology.index_map(1).size_local),
+            partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],
@@ -149,6 +153,7 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
     )
 
     fine_mesh, parent_cell, parent_facet = refine_plaza_wrapper(mesh)
+
     fine_mesh.topology.create_entities(tdim - 1)
     new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet)
     assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices)

From de0e55b095234394c570e11343b49228a447d82c Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 19:05:01 +0100
Subject: [PATCH 24/31] More updates

---
 python/test/unit/refinement/test_refinement.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 719acc8e366..b6fac5f5903 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -180,10 +180,11 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet),
+        lambda mesh: refine(mesh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
         lambda mesh: refine(
             mesh,
             np.arange(mesh.topology.index_map(1).size_local),
+            partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
     ],

From b2feb52d9a436aa06eb9b84a4bcc54ac2af80f89 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 19:23:07 +0100
Subject: [PATCH 25/31] Fix broken tests

---
 python/test/unit/refinement/test_refinement.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index b6fac5f5903..c8b9333a956 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -80,7 +80,7 @@ def left_corner_edge(x, tol=1e-7):
         return isclose(x[0], 0) & (x[1] < 1 / 4 + tol)
 
     edges = locate_entities_boundary(mesh, 1, left_corner_edge)
-    if MPI.COMM_WORLD.size == 0:
+    if MPI.COMM_WORLD.size == 1:
         assert edges == 1
 
     mesh2, _, _ = refine(mesh, edges)
@@ -100,12 +100,12 @@ def left_side(x, tol=1e-16):
         return x[0] <= 0.5 + tol
 
     cells = locate_entities(msh, msh.topology.dim, left_side)
-    if MPI.COMM_WORLD.size == 0:
-        assert cells.__len__() == Nx * Ny
+    if MPI.COMM_WORLD.size == 1:
+        assert len(cells) == Nx * Ny
     edges = compute_incident_entities(msh.topology, cells, 2, 1)
-    if MPI.COMM_WORLD.size == 0:
-        assert edges.__len__() == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny
-    mesh2, _, _ = refine(msh, edges)  # TODO: redistribute=True
+    if MPI.COMM_WORLD.size == 1:
+        assert len(edges) == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny
+    mesh2, _, _ = refine(msh, edges)
 
     num_cells_global = mesh2.topology.index_map(2).size_global
     actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny

From 2854936bd9f7831fa75b4d5c33711195b96e7f64 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 19:33:33 +0100
Subject: [PATCH 26/31] Doc update

---
 cpp/dolfinx/refinement/refine.h |  4 ++--
 cpp/dolfinx/refinement/utils.h  | 38 +++++++++++++++++----------------
 2 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index d92ba24c4ff..dd904820af6 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -24,8 +24,8 @@ namespace dolfinx::refinement
 /// optionally calculating the parent-child relationships.
 ///
 /// @note Using the default partitioner for a refined mesh, the refined
-/// mesh will include ghosts cells (connected by facet) even if the
-/// parent mesh is not ghosted.
+/// mesh will **not** include ghosts cells (connected by facet) even if
+/// the parent mesh is ghosted.
 ///
 /// @warning Passing `nullptr` for `partitioner`, refined cells will be
 /// on the same process as the parent cell but the refined mesh will
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index d0f6d158e4b..9ac12932b79 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -264,14 +264,16 @@ create_new_vertices(MPI_Comm comm,
           xshape};
 }
 
-/// @brief Given an index map, add "n" extra indices at the end of local range
+/// @brief Given an index map, add "n" extra indices at the end of local
+/// range.
 ///
-/// @note The returned global indices (local and ghosts) are adjust for the
-/// addition of new local indices.
-/// @param[in] map Index map
-/// @param[in] n Number of new local indices
-/// @return New global indices for both owned and ghosted indices in input
-/// index map.
+/// @note The returned global indices (local and ghosts) are adjust for
+/// the addition of new local indices.
+///
+/// @param[in] map Index map.
+/// @param[in] n Number of new local indices.
+/// @return New global indices for both owned and ghosted indices in
+/// input index map.
 std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
                                          std::int32_t n);
 
@@ -280,28 +282,28 @@ std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
 /// @warning The refined mesh must not have been redistributed during
 /// refinement.
 ///
-/// @warning GhostMode must be GhostMode.none.
+/// @warning `GhostMode` must be mesh::GhostMode::none.
 ///
-/// @param[in] tags0 Tags on the parent mesh
-/// @param[in] topology1 Refined mesh topology
+/// @param[in] tags0 Tags on the parent mesh.
+/// @param[in] topology1 Refined mesh topology.
 /// @param[in] cell Parent cell of each cell in refined mesh
-/// @param[in] facet Local facets of parent in each cell in refined mesh
-/// @return (0) entities and (1) values on the refined topology
+/// @param[in] facet Local facets of parent in each cell in refined mesh.
+/// @return (0) entities and (1) values on the refined topology.
 std::array<std::vector<std::int32_t>, 2> transfer_facet_meshtag(
     const mesh::MeshTags<std::int32_t>& tags0, const mesh::Topology& topology1,
     std::span<const std::int32_t> cell, std::span<const std::int8_t> facet);
 
 /// @brief Transfer cell MeshTags from coarse mesh to refined mesh.
 ///
-/// @note The refined mesh must not have been redistributed during
+/// @warning The refined mesh must not have been redistributed during
 /// refinement.
 ///
-/// @warning GhostMode must be GhostMode.none.
+/// @warning `GhostMode` must be mesh::GhostMode::none.
 ///
-/// @param[in] tags0 Tags on the parent mesh
-/// @param[in] topology1 Refined mesh topology
-/// @param[in] parent_cell Parent cell of each cell in refined mesh
-/// @return (0) entities and (1) values on the refined topology
+/// @param[in] tags0 Tags on the parent mesh.
+/// @param[in] topology1 Refined mesh topology.
+/// @param[in] parent_cell Parent cell of each cell in refined mesh.
+/// @return (0) entities and (1) values on the refined topology.
 std::array<std::vector<std::int32_t>, 2>
 transfer_cell_meshtag(const mesh::MeshTags<std::int32_t>& tags0,
                       const mesh::Topology& topology1,

From b65857b0effa9fd305a4e29a267342c43e691eeb Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 19:49:54 +0100
Subject: [PATCH 27/31] Update test

---
 python/dolfinx/mesh.py                         | 6 ++----
 python/test/unit/refinement/test_refinement.py | 3 ++-
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 4cc7d001859..9740e3837d6 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -535,9 +535,7 @@ def transfer_meshtag(
 def refine(
     mesh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
-    partitioner: typing.Optional[typing.Callable] = _cpp.mesh.create_cell_partitioner(
-        GhostMode.none
-    ),
+    partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none),
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
@@ -609,7 +607,7 @@ def create_mesh(
         A mesh.
     """
     if partitioner is None and comm.size > 1:
-        partitioner = _cpp.mesh.create_cell_partitioner(GhostMode.none)
+        partitioner = create_cell_partitioner(GhostMode.none)
 
     x = np.asarray(x, order="C")
     if x.ndim == 1:
diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index c8b9333a956..19401e31f4c 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -11,6 +11,7 @@
 from numpy import isclose
 
 import ufl
+from dolfinx.cpp.mesh import create_cell_partitioner
 from dolfinx.fem import assemble_matrix, form, functionspace
 from dolfinx.mesh import (
     CellType,
@@ -46,7 +47,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute):
     """Refine mesh of unit cube."""
     mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode)
     mesh.topology.create_entities(1)
-    mesh, _, _ = refine(mesh)  # , partitioner=create_cell_partitioner(ghost_mode))
+    mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode))
     assert mesh.topology.index_map(0).size_global == 3135
     assert mesh.topology.index_map(3).size_global == 15120
 

From 885ab2d935ec165c4a00aee04c91752f7b4a8a1d Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sat, 28 Sep 2024 19:51:15 +0100
Subject: [PATCH 28/31] Test updates

---
 python/test/unit/refinement/test_refinement.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index 19401e31f4c..b767ba7c58a 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -59,8 +59,7 @@ def test_refine_create_form():
     """Check that forms can be assembled on refined mesh"""
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)
     mesh.topology.create_entities(1)
-    mesh, _, _ = refine(mesh)  # TODO: recover redistribute=True behavior
-
+    mesh, _, _ = refine(mesh)
     V = functionspace(mesh, ("Lagrange", 1))
 
     # Define variational problem

From 7d724286c2196d9b16b9ff93a4f79ee7ba9e176f Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sun, 29 Sep 2024 10:56:42 +0100
Subject: [PATCH 29/31] Work on docs

---
 cpp/dolfinx/refinement/refine.h               |  48 +++++----
 cpp/dolfinx/refinement/utils.h                |   4 +-
 cpp/test/mesh/refinement/interval.cpp         |   3 +-
 python/dolfinx/mesh.py                        | 100 +++++++++---------
 python/dolfinx/wrappers/mesh.cpp              |  15 ++-
 .../test/unit/refinement/test_refinement.py   |  89 +++++++---------
 6 files changed, 131 insertions(+), 128 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index dd904820af6..e784558ce68 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -16,31 +16,39 @@
 #include <algorithm>
 #include <concepts>
 #include <optional>
+#include <spdlog/spdlog.h>
 #include <utility>
 
 namespace dolfinx::refinement
 {
-/// @brief Refine with markers, optionally redistributing, and
-/// optionally calculating the parent-child relationships.
+/// @brief Refine a mesh with markers.
 ///
-/// @note Using the default partitioner for a refined mesh, the refined
-/// mesh will **not** include ghosts cells (connected by facet) even if
-/// the parent mesh is ghosted.
+/// The refined mesh can be optionally re-partitioned across processes.
+/// Passing `nullptr` for `partitioner`, refined cells will be on the
+/// same process as the parent cell.
 ///
-/// @warning Passing `nullptr` for `partitioner`, refined cells will be
-/// on the same process as the parent cell but the refined mesh will
-/// **not** have ghosts cells. The possibility to not re-partition the
-/// refined mesh and include ghost cells in the refined mesh will be
-/// added in a future release.
+/// Parent-child relationships can be optionally computed. Parent-child
+/// relationships can be used to create MeshTags on the refined mesh
+/// from MeshTags on the parent mesh.
+///
+/// @warning Using the default partitioner for a refined mesh, the
+/// refined mesh will **not** include ghosts cells (cells connected by
+/// facet to an owned cell) even if the parent mesh is ghosted.
+///
+/// @warning Passing `nullptr` for `partitioner`, the refined mesh will
+/// **not** have ghosts cells even if the parent mesh is ghosted. The
+/// possibility to not re-partition the refined mesh and include ghost
+/// cells in the refined mesh will be added in a future release.
 ///
 /// @param[in] mesh Input mesh to be refined.
-/// @param[in] edges Indices of the edges that should be split by this
-/// refinement. If not provides, a uniform refinement is performed.
-/// @param[in] partitioner Partitioner to be used for the refined mesh.
-/// If not callable, refined cells will be on the same process as the
-/// parent cell.
+/// @param[in] edges Indices of the edges that should be split in the
+/// refinement. If not provided (`std::nullopt`), uniform refinement is
+/// performed.
+/// @param[in] partitioner Partitioner to be used to distribute the
+/// refined mesh. If not callable, refined cells will be on the same
+/// process as the parent cell.
 /// @param[in] option Control the computation of parent facets, parent
-/// cells. If an option is unselected, an empty list is returned.
+/// cells. If an option is not selected, an empty list is returned.
 /// @return New mesh and optional parent cell index, parent facet
 /// indices.
 template <std::floating_point T>
@@ -62,20 +70,18 @@ refine(const mesh::Mesh<T>& mesh,
             ? interval::compute_refinement_data(mesh, edges, option)
             : plaza::compute_refinement_data(mesh, edges, option);
 
-  mesh::Mesh<T> refined_mesh = mesh::create_mesh(
+  mesh::Mesh<T> mesh1 = mesh::create_mesh(
       mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(),
       mesh.comm(), new_vertex_coords, xshape, partitioner);
 
   // Report the number of refined cells
   const int D = topology->dim();
   const std::int64_t n0 = topology->index_map(D)->size_global();
-  const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global();
+  const std::int64_t n1 = mesh1.topology()->index_map(D)->size_global();
   spdlog::info(
       "Number of cells increased from {} to {} ({}% increase).", n0, n1,
       100.0 * (static_cast<double>(n1) / static_cast<double>(n0) - 1.0));
 
-  return {std::move(refined_mesh), std::move(parent_cell),
-          std::move(parent_facet)};
+  return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)};
 }
-
 } // namespace dolfinx::refinement
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index 9ac12932b79..3b708e9459d 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -264,6 +264,8 @@ create_new_vertices(MPI_Comm comm,
           xshape};
 }
 
+/// @todo Improve docstring.
+///
 /// @brief Given an index map, add "n" extra indices at the end of local
 /// range.
 ///
@@ -277,7 +279,7 @@ create_new_vertices(MPI_Comm comm,
 std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
                                          std::int32_t n);
 
-/// @brief Transfer facet MeshTags from coarse mesh to refined mesh
+/// @brief Transfer facet MeshTags from coarse mesh to refined mesh.
 ///
 /// @warning The refined mesh must not have been redistributed during
 /// refinement.
diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp
index a1ed0d7b9ab..1253fe00b04 100644
--- a/cpp/test/mesh/refinement/interval.cpp
+++ b/cpp/test/mesh/refinement/interval.cpp
@@ -190,14 +190,13 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)",
   mesh::Mesh<T> mesh = create_mesh();
   mesh.topology()->create_connectivity(1, 0);
 
-  // TODO: parent_facet
   auto [refined_mesh, parent_edges, parent_facet] = refinement::refine(
       mesh, std::nullopt, nullptr, refinement::Option::parent_cell);
 
   T rank_d = static_cast<T>(rank);
   T comm_size_d = static_cast<T>(comm_size);
 
-  auto x = refined_mesh.geometry().x();
+  std::span x = refined_mesh.geometry().x();
   std::ranges::sort(x);
   std::vector<T> expected_x
       = {rank_d / comm_size_d,
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 9740e3837d6..3b0ef96b0d6 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -276,11 +276,11 @@ class Mesh:
     _geometry: Geometry
     _ufl_domain: typing.Optional[ufl.Mesh]
 
-    def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
+    def __init__(self, msh, domain: typing.Optional[ufl.Mesh]):
         """Initialize mesh from a C++ mesh.
 
         Args:
-            mesh: A C++ mesh object.
+            msh: A C++ mesh object.
             domain: A UFL domain.
 
         Note:
@@ -289,7 +289,7 @@ def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
             This class is combined with different base classes that
             depend on the scalar type used in the Mesh.
         """
-        self._cpp_object = mesh
+        self._cpp_object = msh
         self._topology = Topology(self._cpp_object.topology)
         self._geometry = Geometry(self._cpp_object.geometry)
         self._ufl_domain = domain
@@ -439,25 +439,25 @@ def compute_incident_entities(
     return _cpp.mesh.compute_incident_entities(topology._cpp_object, entities, d0, d1)
 
 
-def compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]):
+def compute_midpoints(msh: Mesh, dim: int, entities: npt.NDArray[np.int32]):
     """Compute the midpoints of a set of mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entities to consider.
         entities: Indices of entities in ``mesh`` to consider.
 
     Returns:
         Midpoints of the entities, shape ``(num_entities, 3)``.
     """
-    return _cpp.mesh.compute_midpoints(mesh._cpp_object, dim, entities)
+    return _cpp.mesh.compute_midpoints(msh._cpp_object, dim, entities)
 
 
-def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
+def locate_entities(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
     """Compute mesh entities satisfying a geometric marking function.
 
     Args:
-        mesh: Mesh to locate entities on.
+        msh: Mesh to locate entities on.
         dim: Topological dimension of the mesh entities to consider.
         marker: A function that takes an array of points ``x`` with
             shape ``(gdim, num_points)`` and returns an array of
@@ -470,7 +470,7 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray
     return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker)
 
 
-def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
+def locate_entities_boundary(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
     """Compute mesh entities that are connected to an owned boundary
     facet and satisfy a geometric marking function.
 
@@ -484,7 +484,7 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n
     parallel communication.
 
     Args:
-        mesh: Mesh to locate boundary entities on.
+        msh: Mesh to locate boundary entities on.
         dim: Topological dimension of the mesh entities to consider
         marker: Function that takes an array of points ``x`` with shape
             ``(gdim, num_points)`` and returns an array of booleans of
@@ -494,12 +494,12 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n
     Returns:
         Indices (local to the process) of marked mesh entities.
     """
-    return _cpp.mesh.locate_entities_boundary(mesh._cpp_object, dim, marker)
+    return _cpp.mesh.locate_entities_boundary(msh._cpp_object, dim, marker)
 
 
 def transfer_meshtag(
     meshtag: MeshTags,
-    mesh1: Mesh,
+    msh1: Mesh,
     parent_cell: npt.NDArray[np.int32],
     parent_facet: typing.Optional[npt.NDArray[np.int8]] = None,
 ) -> MeshTags:
@@ -507,7 +507,7 @@ def transfer_meshtag(
 
     Args:
         meshtag: Mesh tags on the coarse, parent mesh.
-        mesh1: The refined mesh.
+        msh1: The refined mesh.
         parent_cell: Index of the parent cell for each cell in the
             refined mesh.
         parent_facet: Index of the local parent facet for each cell
@@ -519,13 +519,13 @@ def transfer_meshtag(
     """
     if meshtag.dim == meshtag.topology.dim:
         mt = _cpp.refinement.transfer_cell_meshtag(
-            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell
+            meshtag._cpp_object, msh1.topology._cpp_object, parent_cell
         )
         return MeshTags(mt)
     elif meshtag.dim == meshtag.topology.dim - 1:
         assert parent_facet is not None
         mt = _cpp.refinement.transfer_facet_meshtag(
-            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell, parent_facet
+            meshtag._cpp_object, msh1.topology._cpp_object, parent_cell, parent_facet
         )
         return MeshTags(mt)
     else:
@@ -533,7 +533,7 @@ def transfer_meshtag(
 
 
 def refine(
-    mesh: Mesh,
+    msh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
     partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none),
     option: RefinementOption = RefinementOption.none,
@@ -544,9 +544,9 @@ def refine(
     same process as the parent cell.
 
     Note:
-        Using the default partitioner for a refined mesh, the refined
-        mesh will include ghosts cells (connected by facet) even if the
-        parent mesh does not include ghost cells.
+        Using the default partitioner for the refined mesh, the refined
+        mesh will **not** include ghosts cells (cells connected by facet
+        to an owned cells) even if the parent mesh is ghosted.
 
     Warning:
         Passing ``None`` for ``partitioner``, the refined mesh will
@@ -556,12 +556,12 @@ def refine(
         future release.
 
     Args:
-        mesh: Mesh from which to create the refined mesh.
+        msh: Mesh from which to create the refined mesh.
         edges: Indices of edges to split during refinement. If ``None``,
             mesh refinement is uniform.
         partitioner: Partitioner to distribute the refined mesh. If
-            ``None`` no redistribution is performed, i.e. refined cells remain on the
-            same process as the parent cell.
+            ``None`` no redistribution is performed, i.e. refined cells
+            remain on the same process as the parent cell.
         option: Controls whether parent cells and/or parent facets are
             computed.
 
@@ -569,10 +569,11 @@ def refine(
        Refined mesh, (optional) parent cells, (optional) parent facets
     """
     mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
-        mesh._cpp_object, edges, partitioner, option
+        msh._cpp_object, edges, partitioner, option
     )
-    # Create new ufl domain as it will carry a reference to the C++ mesh in the ufl_cargo
-    ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element())  # type: ignore
+    # Create new ufl domain as it will carry a reference to the C++ mesh
+    # in the ufl_cargo
+    ufl_domain = ufl.Mesh(msh._ufl_domain.ufl_coordinate_element())  # type: ignore
     return Mesh(mesh1, ufl_domain), parent_cell, parent_facet
 
 
@@ -596,8 +597,8 @@ def create_mesh(
         x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, gdim)``.
         e: UFL mesh. The mesh scalar type is determined by the scalar
             type of ``e``.
-        partitioner: Function that computes the parallel distribution of
-            cells across MPI ranks.
+        partitioner: Function that determines the parallel distribution
+            of cells across MPI ranks.
 
     Note:
         If required, the coordinates ``x`` will be cast to the same
@@ -649,9 +650,9 @@ def create_mesh(
 
     x = np.asarray(x, dtype=dtype, order="C")
     cells = np.asarray(cells, dtype=np.int64, order="C")
-    mesh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner)
+    msh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner)
 
-    return Mesh(mesh, domain)
+    return Mesh(msh, domain)
 
 
 def create_submesh(
@@ -660,7 +661,7 @@ def create_submesh(
     """Create a mesh with specified entities from another mesh.
 
     Args:
-        mesh: Mesh to create the sub-mesh from.
+        msh: Mesh to create the sub-mesh from.
         dim: Topological dimension of the entities in ``msh`` to include
             in the sub-mesh.
         entities: Indices of entities in ``msh`` to include in the
@@ -689,7 +690,7 @@ def create_submesh(
 
 
 def meshtags(
-    mesh: Mesh,
+    msh: Mesh,
     dim: int,
     entities: npt.NDArray[np.int32],
     values: typing.Union[np.ndarray, int, float],
@@ -697,7 +698,7 @@ def meshtags(
     """Create a MeshTags object that associates data with a subset of mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entity.
         entities: Indices(local to process) of entities to associate
             values with . The array must be sorted and must not contain
@@ -730,19 +731,19 @@ def meshtags(
         raise NotImplementedError(f"Type {values.dtype} not supported.")
 
     return MeshTags(
-        ftype(mesh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values)
+        ftype(msh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values)
     )
 
 
 def meshtags_from_entities(
-    mesh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any]
+    msh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any]
 ):
     """Create a :class:dolfinx.mesh.MeshTags` object that associates
     data with a subset of mesh entities, where the entities are defined
     by their vertices.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the mesh entity.
         entities: Entities to associated values with, with entities
             defined by their vertices.
@@ -762,7 +763,7 @@ def meshtags_from_entities(
     elif isinstance(values, float):
         values = np.full(entities.num_nodes, values, dtype=np.double)
     values = np.asarray(values)
-    return MeshTags(_cpp.mesh.create_meshtags(mesh.topology._cpp_object, dim, entities, values))
+    return MeshTags(_cpp.mesh.create_meshtags(msh.topology._cpp_object, dim, entities, values))
 
 
 def create_interval(
@@ -793,12 +794,13 @@ def create_interval(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", "interval", 1, shape=(1,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner)
+        msh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner)
+        msh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_interval(
@@ -861,12 +863,13 @@ def create_rectangle(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(2,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal)
+        msh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal)
+        msh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_square(
@@ -939,12 +942,13 @@ def create_box(
         partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)
     domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(3,), dtype=dtype))  # type: ignore
     if np.issubdtype(dtype, np.float32):
-        mesh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner)
+        msh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner)
     elif np.issubdtype(dtype, np.float64):
-        mesh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner)
+        msh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner)
     else:
         raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
-    return Mesh(mesh, domain)
+
+    return Mesh(msh, domain)
 
 
 def create_unit_cube(
@@ -987,12 +991,12 @@ def create_unit_cube(
 
 
 def entities_to_geometry(
-    mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False
+    msh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False
 ) -> npt.NDArray[np.int32]:
     """Compute the geometric DOFs associated with the closure of the given mesh entities.
 
     Args:
-        mesh: The mesh.
+        msh: The mesh.
         dim: Topological dimension of the entities of interest.
         entities: Entity indices (local to the process).
         permute: Permute the DOFs such that they are consistent with the
@@ -1003,7 +1007,7 @@ def entities_to_geometry(
         The geometric DOFs associated with the closure of the entities
         in `entities`.
     """
-    return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, permute)
+    return _cpp.mesh.entities_to_geometry(msh._cpp_object, dim, entities, permute)
 
 
 def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]:
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index e2cae585563..5cf66c9c6ab 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -238,11 +238,11 @@ void declare_mesh(nb::module_& m, std::string type)
   m.def(
       create_interval.c_str(),
       [](MPICommWrapper comm, std::int64_t n, std::array<T, 2> p,
-         dolfinx::mesh::GhostMode ghost_mode,
+         dolfinx::mesh::GhostMode mode,
          const part::impl::PythonCellPartitionFunction& part)
       {
         return dolfinx::mesh::create_interval<T>(
-            comm.get(), n, p, ghost_mode,
+            comm.get(), n, p, mode,
             part::impl::create_cell_partitioner_cpp(part));
       },
       nb::arg("comm"), nb::arg("n"), nb::arg("p"), nb::arg("ghost_mode"),
@@ -683,8 +683,6 @@ void mesh(nb::module_& m)
               ghost_owners_span, boundary_vertices_span);
         });
 
-  // dolfinx::mesh::MeshTags
-
   declare_meshtags<std::int8_t>(m, "int8");
   declare_meshtags<std::int32_t>(m, "int32");
   declare_meshtags<std::int64_t>(m, "int64");
@@ -695,10 +693,11 @@ void mesh(nb::module_& m)
 
   m.def(
       "create_cell_partitioner",
-      [](dolfinx::mesh::GhostMode gm) -> part::impl::PythonCellPartitionFunction
+      [](dolfinx::mesh::GhostMode mode)
+          -> part::impl::PythonCellPartitionFunction
       {
         return part::impl::create_cell_partitioner_py(
-            dolfinx::mesh::create_cell_partitioner(gm));
+            dolfinx::mesh::create_cell_partitioner(mode));
       },
       "Create default cell partitioner.");
   m.def(
@@ -708,12 +707,12 @@ void mesh(nb::module_& m)
              const dolfinx::graph::AdjacencyList<std::int64_t>& local_graph,
              bool ghosting)>
              part,
-         dolfinx::mesh::GhostMode ghost_mode)
+         dolfinx::mesh::GhostMode mode)
           -> part::impl::PythonCellPartitionFunction
       {
         return part::impl::create_cell_partitioner_py(
             dolfinx::mesh::create_cell_partitioner(
-                ghost_mode, part::impl::create_partitioner_cpp(part)));
+                mode, part::impl::create_partitioner_cpp(part)));
       },
       nb::arg("part"), nb::arg("ghost_mode") = dolfinx::mesh::GhostMode::none,
       "Create a cell partitioner from a graph partitioning function.");
diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py
index b767ba7c58a..0a8ee67b79a 100644
--- a/python/test/unit/refinement/test_refinement.py
+++ b/python/test/unit/refinement/test_refinement.py
@@ -71,20 +71,20 @@ def test_refine_create_form():
 
 def test_sub_refine():
     """Test that refinement of a subset of edges works"""
-    mesh = create_unit_square(
+    msh = create_unit_square(
         MPI.COMM_WORLD, 3, 4, diagonal=DiagonalType.left, ghost_mode=GhostMode.none
     )
-    mesh.topology.create_entities(1)
+    msh.topology.create_entities(1)
 
     def left_corner_edge(x, tol=1e-7):
         return isclose(x[0], 0) & (x[1] < 1 / 4 + tol)
 
-    edges = locate_entities_boundary(mesh, 1, left_corner_edge)
+    edges = locate_entities_boundary(msh, 1, left_corner_edge)
     if MPI.COMM_WORLD.size == 1:
         assert edges == 1
 
-    mesh2, _, _ = refine(mesh, edges)
-    assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global
+    msh1, _, _ = refine(msh, edges)
+    assert msh.topology.index_map(2).size_global + 3 == msh1.topology.index_map(2).size_global
 
 
 def test_refine_from_cells():
@@ -112,17 +112,14 @@ def left_side(x, tol=1e-16):
     assert num_cells_global == actual_cells
 
 
-@pytest.mark.parametrize(
-    "tdim",
-    [2, 3],
-)
+@pytest.mark.parametrize("tdim", [2, 3])
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(mesh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
-        lambda mesh: refine(
-            mesh,
-            edges=np.arange(mesh.topology.index_map(1).size_local),
+        lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
+        lambda msh: refine(
+            msh,
+            edges=np.arange(msh.topology.index_map(1).size_local),
             partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
@@ -130,49 +127,47 @@ def left_side(x, tol=1e-16):
 )
 def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
     if tdim == 3:
-        mesh = create_unit_cube(
+        msh = create_unit_cube(
             MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none
         )
     else:
-        mesh = create_unit_square(
-            MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none
-        )
-    mesh.topology.create_entities(tdim - 1)
-    mesh.topology.create_connectivity(tdim - 1, tdim)
-    mesh.topology.create_entities(1)
-    f_to_c = mesh.topology.connectivity(tdim - 1, tdim)
+        msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none)
+    msh.topology.create_entities(tdim - 1)
+    msh.topology.create_connectivity(tdim - 1, tdim)
+    msh.topology.create_entities(1)
+    f_to_c = msh.topology.connectivity(tdim - 1, tdim)
     facet_indices = []
-    for f in range(mesh.topology.index_map(tdim - 1).size_local):
+    for f in range(msh.topology.index_map(tdim - 1).size_local):
         if len(f_to_c.links(f)) == 1:
             facet_indices += [f]
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim - 1,
         np.array(facet_indices, dtype=np.int32),
         np.arange(len(facet_indices), dtype=np.int32),
     )
 
-    fine_mesh, parent_cell, parent_facet = refine_plaza_wrapper(mesh)
+    msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh)
 
-    fine_mesh.topology.create_entities(tdim - 1)
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet)
+    msh1.topology.create_entities(tdim - 1)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet)
     assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices)
 
     # New tags should be on facets with one cell (i.e. exterior)
-    fine_mesh.topology.create_connectivity(tdim - 1, tdim)
-    new_f_to_c = fine_mesh.topology.connectivity(tdim - 1, tdim)
+    msh1.topology.create_connectivity(tdim - 1, tdim)
+    new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim)
     for f in new_meshtag.indices:
         assert len(new_f_to_c.links(f)) == 1
 
     # Now mark all facets (including internal)
-    facet_indices = np.arange(mesh.topology.index_map(tdim - 1).size_local)
+    facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local)
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim - 1,
         np.array(facet_indices, dtype=np.int32),
         np.arange(len(facet_indices), dtype=np.int32),
     )
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet)
     assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices)
 
 
@@ -180,10 +175,10 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 @pytest.mark.parametrize(
     "refine_plaza_wrapper",
     [
-        lambda mesh: refine(mesh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
-        lambda mesh: refine(
-            mesh,
-            np.arange(mesh.topology.index_map(1).size_local),
+        lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet),
+        lambda msh: refine(
+            msh,
+            np.arange(msh.topology.index_map(1).size_local),
             partitioner=None,
             option=RefinementOption.parent_cell_and_facet,
         ),
@@ -191,31 +186,29 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper):
 )
 def test_refine_cell_meshtag(tdim, refine_plaza_wrapper):
     if tdim == 3:
-        mesh = create_unit_cube(
+        msh = create_unit_cube(
             MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none
         )
     else:
-        mesh = create_unit_square(
-            MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none
-        )
+        msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none)
 
-    mesh.topology.create_entities(1)
-    cell_indices = np.arange(mesh.topology.index_map(tdim).size_local)
+    msh.topology.create_entities(1)
+    cell_indices = np.arange(msh.topology.index_map(tdim).size_local)
     meshtag = meshtags(
-        mesh,
+        msh,
         tdim,
         np.array(cell_indices, dtype=np.int32),
         np.arange(len(cell_indices), dtype=np.int32),
     )
 
-    fine_mesh, parent_cell, _ = refine_plaza_wrapper(mesh)
-    new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell)
+    msh1, parent_cell, _ = refine_plaza_wrapper(msh)
+    new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell)
     assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values)
     assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices)
 
 
 def test_refine_ufl_cargo():
-    mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3)
-    mesh.topology.create_entities(1)
-    refined_mesh, _, _ = refine(mesh)
-    assert refined_mesh.ufl_domain().ufl_cargo() != mesh.ufl_domain().ufl_cargo()
+    msh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3)
+    msh.topology.create_entities(1)
+    msh1, _, _ = refine(msh)
+    assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo()

From bd6ca515119a68dee3b26afe51fed84b042acc9f Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Sun, 29 Sep 2024 12:38:38 +0100
Subject: [PATCH 30/31] Lint

---
 python/dolfinx/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 3b0ef96b0d6..2203ca8d961 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -467,7 +467,7 @@ def locate_entities(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
     Returns:
         Indices (local to the process) of marked mesh entities.
     """
-    return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker)
+    return _cpp.mesh.locate_entities(msh._cpp_object, dim, marker)
 
 
 def locate_entities_boundary(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:

From 9e58656d2427865917b221306d4f9714d4bd3fde Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Tue, 1 Oct 2024 07:37:11 +0100
Subject: [PATCH 31/31] Minor doc updates

---
 cpp/dolfinx/refinement/refine.h | 2 +-
 cpp/dolfinx/refinement/utils.h  | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index e784558ce68..430bcc495bc 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -49,7 +49,7 @@ namespace dolfinx::refinement
 /// process as the parent cell.
 /// @param[in] option Control the computation of parent facets, parent
 /// cells. If an option is not selected, an empty list is returned.
-/// @return New mesh and optional parent cell index, parent facet
+/// @return New mesh, and optional parent cell indices and parent facet
 /// indices.
 template <std::floating_point T>
 std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h
index 3b708e9459d..ec550cb3575 100644
--- a/cpp/dolfinx/refinement/utils.h
+++ b/cpp/dolfinx/refinement/utils.h
@@ -284,7 +284,7 @@ std::vector<std::int64_t> adjust_indices(const common::IndexMap& map,
 /// @warning The refined mesh must not have been redistributed during
 /// refinement.
 ///
-/// @warning `GhostMode` must be mesh::GhostMode::none.
+/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none.
 ///
 /// @param[in] tags0 Tags on the parent mesh.
 /// @param[in] topology1 Refined mesh topology.
@@ -300,7 +300,7 @@ std::array<std::vector<std::int32_t>, 2> transfer_facet_meshtag(
 /// @warning The refined mesh must not have been redistributed during
 /// refinement.
 ///
-/// @warning `GhostMode` must be mesh::GhostMode::none.
+/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none.
 ///
 /// @param[in] tags0 Tags on the parent mesh.
 /// @param[in] topology1 Refined mesh topology.