diff --git a/cpp/doc/source/api.rst b/cpp/doc/source/api.rst
index f9b914c732b..8464594f7a2 100644
--- a/cpp/doc/source/api.rst
+++ b/cpp/doc/source/api.rst
@@ -14,6 +14,7 @@ generated documentation is `here <doxygen>`_.
    io
    la
    mesh
+   multigrid
    refinement
 
 * :ref:`genindex`
diff --git a/cpp/doc/source/multigrid.rst b/cpp/doc/source/multigrid.rst
new file mode 100644
index 00000000000..f0fe0e6ecb1
--- /dev/null
+++ b/cpp/doc/source/multigrid.rst
@@ -0,0 +1,5 @@
+Multigrid (``dolfinx::multigrid``)
+====================================
+
+.. doxygennamespace:: dolfinx::multigrid
+   :project: DOLFINx
diff --git a/cpp/dolfinx/CMakeLists.txt b/cpp/dolfinx/CMakeLists.txt
index 7adaaeb73cf..001cb2116be 100644
--- a/cpp/dolfinx/CMakeLists.txt
+++ b/cpp/dolfinx/CMakeLists.txt
@@ -17,6 +17,7 @@ set(DOLFINX_DIRS
     io
     la
     mesh
+    multigrid
     nls
     refinement
 )
diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp
index fade5e160cf..425ca4e2910 100644
--- a/cpp/dolfinx/graph/partitioners.cpp
+++ b/cpp/dolfinx/graph/partitioners.cpp
@@ -35,21 +35,8 @@ extern "C"
 
 using namespace dolfinx;
 
-namespace
-{
-/// @todo Is it un-documented that the owning rank must come first in
-/// reach list of edges?
-///
-/// @param[in] comm The communicator
-/// @param[in] graph Graph, using global indices for graph edges
-/// @param[in] node_disp The distribution of graph nodes across MPI
-/// ranks. The global index `gidx` of local index `lidx` is `lidx +
-/// node_disp[my_rank]`.
-/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]`
-/// is the destination of the node with local index `i`.
-/// @return Destination ranks for each local node.
 template <typename T>
-graph::AdjacencyList<int> compute_destination_ranks(
+graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
     MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
     const std::vector<T>& node_disp, const std::vector<T>& part)
 {
@@ -202,6 +189,18 @@ graph::AdjacencyList<int> compute_destination_ranks(
 
   return g;
 }
+
+/// @cond
+template graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<int>& node_disp, const std::vector<int>& part);
+
+template graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<unsigned long long>& node_disp,
+    const std::vector<unsigned long long>& part);
+/// @endcond
+
 //-----------------------------------------------------------------------------
 #ifdef HAS_PARMETIS
 template <typename T>
@@ -304,7 +303,6 @@ std::vector<int> refine(MPI_Comm comm, const graph::AdjacencyList<T>& adj_graph)
   //-----------------------------------------------------------------------------
 }
 #endif
-} // namespace
 
 //-----------------------------------------------------------------------------
 #ifdef HAS_PTSCOTCH
@@ -587,7 +585,7 @@ graph::partition_fn graph::parmetis::partitioner(double imbalance,
     {
       // FIXME: Is it implicit that the first entry is the owner?
       graph::AdjacencyList<int> dest
-          = compute_destination_ranks(pcomm, graph, node_disp, part);
+          = graph::compute_destination_ranks(pcomm, graph, node_disp, part);
       if (split_comm)
         MPI_Comm_free(&pcomm);
       return dest;
@@ -653,7 +651,7 @@ graph::partition_fn graph::kahip::partitioner(int mode, int seed,
     timer2.stop();
 
     if (ghosting)
-      return compute_destination_ranks(comm, graph, node_disp, part);
+      return graph::compute_destination_ranks(comm, graph, node_disp, part);
     else
     {
       return regular_adjacency_list(std::vector<int>(part.begin(), part.end()),
diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h
index aba51a8430a..26a78309790 100644
--- a/cpp/dolfinx/graph/partitioners.h
+++ b/cpp/dolfinx/graph/partitioners.h
@@ -11,6 +11,23 @@
 
 namespace dolfinx::graph
 {
+
+/// @todo Is it un-documented that the owning rank must come first in
+/// reach list of edges?
+///
+/// @param[in] comm The communicator
+/// @param[in] graph Graph, using global indices for graph edges
+/// @param[in] node_disp The distribution of graph nodes across MPI
+/// ranks. The global index `gidx` of local index `lidx` is `lidx +
+/// node_disp[my_rank]`.
+/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]`
+/// is the destination of the node with local index `i`.
+/// @return Destination ranks for each local node.
+template <typename T>
+graph::AdjacencyList<int> compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<T>& node_disp, const std::vector<T>& part);
+
 namespace scotch
 {
 #ifdef HAS_PTSCOTCH
diff --git a/cpp/dolfinx/multigrid/CMakeLists.txt b/cpp/dolfinx/multigrid/CMakeLists.txt
new file mode 100644
index 00000000000..55516dc1ef0
--- /dev/null
+++ b/cpp/dolfinx/multigrid/CMakeLists.txt
@@ -0,0 +1,5 @@
+set(HEADERS_multigrid
+    ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h
+    PARENT_SCOPE
+)
diff --git a/cpp/dolfinx/multigrid/dolfinx_multigrid.h b/cpp/dolfinx/multigrid/dolfinx_multigrid.h
new file mode 100644
index 00000000000..5dff5c3502f
--- /dev/null
+++ b/cpp/dolfinx/multigrid/dolfinx_multigrid.h
@@ -0,0 +1,12 @@
+#pragma once
+
+/// @brief Multigrid algorithms.
+///
+namespace dolfinx::multigrid
+{
+}
+
+// DOLFINx multigrid interface
+
+#include <dolfinx/multigrid/inclusion.h>
+#include <dolfinx/multigrid/transfer_matrix.h>
diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h
new file mode 100644
index 00000000000..69f543e8e95
--- /dev/null
+++ b/cpp/dolfinx/multigrid/inclusion.h
@@ -0,0 +1,186 @@
+// 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
+
+#pragma once
+
+#include <algorithm>
+#include <concepts>
+#include <cstdint>
+#include <iterator>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+
+#include <mpi.h>
+
+#include "dolfinx/common/IndexMap.h"
+#include "dolfinx/la/SparsityPattern.h"
+#include "dolfinx/mesh/Mesh.h"
+
+namespace dolfinx::multigrid
+{
+
+/**
+ * @brief Gathers a global vector from combination of local data.
+ * @note Performs an all-to-all communication.
+ *
+ * @param local local data
+ * @param global_size number of global data entries
+ * @param comm MPI communicator
+ * @return std::vector<T> on communicator gathered global data.
+ */
+template <std::floating_point T>
+std::vector<T> gather_global(std::span<const T> local, std::int64_t global_size,
+                             MPI_Comm comm)
+{
+  // 1) exchange local sizes
+  std::vector<std::int32_t> local_sizes(dolfinx::MPI::size(comm));
+  {
+    std::array<std::int32_t, 1> tmp{local.size()};
+    MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T,
+                  comm);
+  }
+
+  // 2) compute displacement vector
+  std::vector<std::int32_t> displacements(local_sizes.size() + 1, 0);
+  std::partial_sum(local_sizes.begin(), local_sizes.end(),
+                   displacements.begin() + 1);
+
+  // 3) Allgather global vector
+  std::vector<T> global(global_size);
+  MPI_Allgatherv(local.data(), local.size(), dolfinx::MPI::mpi_t<T>,
+                 global.data(), local_sizes.data(), displacements.data(),
+                 dolfinx::MPI::mpi_t<T>, comm);
+
+  return global;
+}
+
+/**
+ * @brief Computes an inclusion map, i.e. local list of global vertex indices of
+ * another mesh, between to meshes.
+ *
+ *
+ * @param mesh_from Coarser mesh (domain of the inclusion map)
+ * @param mesh_to Finer mesh (range of the inclusion map)
+ * @param allow_all_to_all If the vertices of `mesh_from` are not equally
+ * spatially parallelized as `mesh_to` an all-to-all gathering of all vertices
+ * in `mesh_to` is performed. If true, performs all-to-all gathering, otherwise
+ * throws an exception if this becomes necessary.
+ * @return std::vector<std::int64_t> Map from local vertex index in `mesh_from`
+ * to global vertex index in `mesh_to`, i.e. `mesh_from.geometry.x()[i:i+3] ==
+ * mesh_to.geometry.x()[map[i]:map[i]+3]`.
+ */
+template <std::floating_point T>
+std::vector<std::int64_t>
+inclusion_mapping(const dolfinx::mesh::Mesh<T>& mesh_from,
+                  const dolfinx::mesh::Mesh<T>& mesh_to,
+                  bool allow_all_to_all = false)
+{
+  {
+    // Check comms equal
+    int result;
+    MPI_Comm_compare(mesh_from.comm(), mesh_to.comm(), &result);
+    assert(result == MPI_CONGRUENT);
+  }
+
+  const common::IndexMap& im_from = *mesh_from.topology()->index_map(0);
+  const common::IndexMap& im_to = *mesh_to.topology()->index_map(0);
+
+  std::vector<std::int64_t> map(im_from.size_local() + im_from.num_ghosts(),
+                                -1);
+
+  std::span<const T> x_from = mesh_from.geometry().x();
+  std::span<const T> x_to = mesh_to.geometry().x();
+
+  auto build_global_to_local = [&](const auto& im)
+  {
+    return [&](std::int32_t idx)
+    {
+      std::array<std::int64_t, 1> tmp;
+      im.local_to_global(std::vector<std::int32_t>{idx}, tmp);
+      return tmp[0];
+    };
+  };
+
+  auto to_global_to = build_global_to_local(im_to);
+  auto to_global_from = build_global_to_local(im_from);
+
+  for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++)
+  {
+    std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i),
+                                      std::next(x_from.begin(), 3 * (i + 1)));
+    for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++)
+    {
+      std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j),
+                                      std::next(x_to.begin(), 3 * (j + 1)));
+
+      if (std::ranges::equal(
+              vertex_from, vertex_to, [](T a, T b)
+              { return std::abs(a - b) <= std::numeric_limits<T>::epsilon(); }))
+      {
+        assert(map[i] == -1);
+        map[i] = to_global_to(j);
+        break;
+      }
+    }
+  }
+
+  if (dolfinx::MPI::size(mesh_to.comm()) == 1)
+  {
+    // no communication required
+    assert(std::ranges::all_of(map, [](auto e) { return e >= 0; }));
+    return map;
+  }
+
+  bool all_found = std::ranges::all_of(map, [](auto e) { return e >= 0; });
+  MPI_Allreduce(MPI_IN_PLACE, &all_found, 1, MPI_CXX_BOOL, MPI_LAND,
+                mesh_to.comm());
+  if (all_found)
+    return map;
+
+  if (!allow_all_to_all)
+  {
+    throw std::runtime_error(
+        "Parallelization of mesh requires all to all communication to compute "
+        "inclusion map, but allow_all_to_all is not set.");
+  }
+
+  // Build global to vertex list
+  std::vector<T> global_x_to
+      = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3),
+                      im_to.size_global() * 3, mesh_to.comm());
+
+  // Recheck indices on global data structure
+  for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++)
+  {
+    // TODO:
+    // if (map[i] >= 0)
+    //   continue;
+
+    std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i),
+                                      std::next(x_from.begin(), 3 * (i + 1)));
+    for (std::int64_t j = 0; j < im_to.size_global(); j++)
+    {
+      std::ranges::subrange vertex_to(
+          std::next(global_x_to.begin(), 3 * j),
+          std::next(global_x_to.begin(), 3 * (j + 1)));
+
+      if (std::ranges::equal(
+              vertex_from, vertex_to, [](T a, T b)
+              { return std::abs(a - b) <= std::numeric_limits<T>::epsilon(); }))
+      {
+        map[i] = j;
+        break;
+      }
+    }
+  }
+
+  assert(std::ranges::all_of(map, [&](auto e)
+                             { return e >= 0 && e < im_to.size_global(); }));
+  return map;
+}
+
+} // namespace dolfinx::multigrid
diff --git a/cpp/dolfinx/multigrid/transfer_matrix.h b/cpp/dolfinx/multigrid/transfer_matrix.h
new file mode 100644
index 00000000000..fc828875666
--- /dev/null
+++ b/cpp/dolfinx/multigrid/transfer_matrix.h
@@ -0,0 +1,193 @@
+// 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
+
+#pragma once
+
+#include <algorithm>
+#include <concepts>
+#include <cstdint>
+#include <stdexcept>
+#include <vector>
+
+#include <basix/element-families.h>
+
+#include "dolfinx/common/IndexMap.h"
+#include "dolfinx/fem/FunctionSpace.h"
+#include "dolfinx/la/SparsityPattern.h"
+#include "dolfinx/la/utils.h"
+
+namespace dolfinx::multigrid
+{
+
+template <std::floating_point T>
+la::SparsityPattern
+create_sparsity_pattern(const dolfinx::fem::FunctionSpace<T>& V_from,
+                        const dolfinx::fem::FunctionSpace<T>& V_to,
+                        const std::vector<std::int64_t>& inclusion_map)
+{
+  if (*V_from.element() != *V_to.element())
+    throw std::runtime_error(
+        "Transfer between different element types not supported");
+
+  if (V_from.element()->basix_element().family() != basix::element::family::P)
+    throw std::runtime_error("Only Lagrange elements supported");
+
+  if (V_from.element()->basix_element().degree() != 1)
+    throw std::runtime_error("Only piecewise linear elements supported");
+
+  // TODO: mixed elements? value shapes? DG?
+
+  auto mesh_from = V_from.mesh();
+  auto mesh_to = V_to.mesh();
+  assert(mesh_from);
+  assert(mesh_to);
+
+  MPI_Comm comm = mesh_from->comm();
+  {
+    // Check comms equal
+    int result;
+    MPI_Comm_compare(comm, mesh_to->comm(), &result);
+    assert(result == MPI_CONGRUENT);
+  }
+  assert(mesh_from->topology()->dim() == mesh_to->topology()->dim());
+
+  auto to_v_to_f = mesh_to->topology()->connectivity(0, 1);
+  auto to_f_to_v = mesh_to->topology()->connectivity(1, 0);
+  assert(to_v_to_f);
+  assert(to_f_to_v);
+
+  auto dofmap_from = V_from.dofmap();
+  auto dofmap_to = V_to.dofmap();
+  assert(dofmap_from);
+  assert(dofmap_to);
+
+  assert(mesh_to->topology()->index_map(0));
+  assert(mesh_from->topology()->index_map(0));
+  const common::IndexMap& im_to = *mesh_to->topology()->index_map(0);
+  const common::IndexMap& im_from = *mesh_from->topology()->index_map(0);
+
+  dolfinx::la::SparsityPattern sp(
+      comm, {dofmap_from->index_map, dofmap_to->index_map},
+      {dofmap_from->index_map_bs(), dofmap_to->index_map_bs()});
+
+  assert(inclusion_map.size() == im_from.size_global());
+  for (int dof_from_global = 0; dof_from_global < im_from.size_global();
+       dof_from_global++)
+  {
+    std::int64_t dof_to_global = inclusion_map[dof_from_global];
+
+    std::vector<std::int32_t> local_dof_to_v{0};
+    im_to.global_to_local(std::vector<std::int64_t>{dof_to_global},
+                          local_dof_to_v);
+
+    auto local_dof_to = local_dof_to_v[0];
+
+    bool is_remote = (local_dof_to == -1);
+    bool is_ghost = local_dof_to >= im_to.size_local();
+    if (is_remote || is_ghost)
+      continue;
+
+    std::vector<std::int32_t> dof_from_v{0};
+    im_from.global_to_local(std::vector<std::int64_t>{dof_from_global},
+                            dof_from_v);
+
+    std::ranges::for_each(
+        to_v_to_f->links(local_dof_to),
+        [&](auto e)
+        {
+          sp.insert(
+              std::vector<int32_t>(to_f_to_v->links(e).size(), dof_from_v[0]),
+              to_f_to_v->links(e));
+        });
+  }
+  sp.finalize();
+  return sp;
+}
+
+template <std::floating_point T>
+void assemble_transfer_matrix(la::MatSet<T> auto mat_set,
+                              const dolfinx::fem::FunctionSpace<T>& V_from,
+                              const dolfinx::fem::FunctionSpace<T>& V_to,
+                              const std::vector<std::int64_t>& inclusion_map,
+                              std::function<T(std::int32_t)> weight)
+{
+  if (*V_from.element() != *V_to.element())
+    throw std::runtime_error(
+        "Transfer between different element types not supported");
+
+  if (V_from.element()->basix_element().family() != basix::element::family::P)
+    throw std::runtime_error("Only Lagrange elements supported");
+
+  if (V_from.element()->basix_element().degree() != 1)
+    throw std::runtime_error("Only piecewise linear elements supported");
+
+  // TODO: mixed elements? value shapes? DG?
+
+  auto mesh_from = V_from.mesh();
+  auto mesh_to = V_to.mesh();
+  assert(mesh_from);
+
+  MPI_Comm comm = mesh_from->comm();
+  {
+    // Check comms equal
+    int result;
+    MPI_Comm_compare(comm, mesh_to->comm(), &result);
+    assert(result == MPI_CONGRUENT);
+  }
+  assert(mesh_from->topology()->dim() == mesh_to->topology()->dim());
+
+  auto to_v_to_f = mesh_to->topology()->connectivity(0, 1);
+  auto to_f_to_v = mesh_to->topology()->connectivity(1, 0);
+  assert(to_v_to_f);
+  assert(to_f_to_v);
+
+  auto dofmap_from = V_from.dofmap();
+  auto dofmap_to = V_to.dofmap();
+  assert(dofmap_from);
+  assert(dofmap_to);
+
+  assert(mesh_to->topology()->index_map(0));
+  assert(mesh_from->topology()->index_map(0));
+  const common::IndexMap& im_to = *mesh_to->topology()->index_map(0);
+  const common::IndexMap& im_from = *mesh_from->topology()->index_map(0);
+
+  assert(inclusion_map.size() == im_from.size_global());
+
+  for (int dof_from_global = 0; dof_from_global < im_from.size_global();
+       dof_from_global++)
+  {
+    std::int64_t dof_to_global = inclusion_map[dof_from_global];
+
+    std::vector<std::int32_t> local_dof_to_v{0};
+    im_to.global_to_local(std::vector<std::int64_t>{dof_to_global},
+                          local_dof_to_v);
+
+    auto local_dof_to = local_dof_to_v[0];
+
+    bool is_remote = (local_dof_to == -1);
+    bool is_ghost = local_dof_to >= im_to.size_local();
+    if (is_remote || is_ghost)
+      continue;
+
+    std::vector<std::int32_t> dof_from_v{0};
+    im_from.global_to_local(std::vector<std::int64_t>{dof_from_global},
+                            dof_from_v);
+
+    for (auto e : to_v_to_f->links(local_dof_to))
+    {
+      for (auto n : to_f_to_v->links(e))
+      {
+        // For now we only support distance <= 1 -> this should be somehow
+        // asserted.
+        std::int32_t distance = (n == local_dof_to) ? 0 : 1;
+        mat_set(std::vector<int32_t>{dof_from_v[0]}, std::vector<int32_t>{n},
+                std::vector<T>{weight(distance)});
+      }
+    }
+  }
+}
+
+} // namespace dolfinx::multigrid
\ No newline at end of file
diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h
index 430bcc495bc..8675f16876b 100644
--- a/cpp/dolfinx/refinement/refine.h
+++ b/cpp/dolfinx/refinement/refine.h
@@ -6,21 +6,89 @@
 
 #pragma once
 
+#include "dolfinx/common/MPI.h"
 #include "dolfinx/graph/AdjacencyList.h"
+#include "dolfinx/graph/partitioners.h"
 #include "dolfinx/mesh/Mesh.h"
 #include "dolfinx/mesh/Topology.h"
 #include "dolfinx/mesh/cell_types.h"
+#include "dolfinx/mesh/graphbuild.h"
 #include "dolfinx/mesh/utils.h"
 #include "interval.h"
 #include "plaza.h"
 #include <algorithm>
 #include <concepts>
+#include <mpi.h>
 #include <optional>
 #include <spdlog/spdlog.h>
+#include <stdexcept>
 #include <utility>
 
 namespace dolfinx::refinement
 {
+
+/**
+ * @brief Create a cell partitioner which maintains the partition of a coarse
+ * mesh.
+ *
+ * @tparam T floating point type of mesh geometry.
+ * @param parent_mesh mesh indicating the partition scheme to use, i.e. the
+ * coarse mesh.
+ * @param parent_cell list of cell indices mapping cells of the new refined mesh
+ * into the coarse mesh, usually output of `refinement::refine`.
+ * @return The created cell partition function.
+ */
+template <std::floating_point T>
+mesh::CellPartitionFunction
+create_identity_partitioner(const mesh::Mesh<T>& parent_mesh,
+                            std::span<std::int32_t> parent_cell)
+{
+  // TODO: optimize for non ghosted mesh?
+
+  return [&parent_mesh,
+          parent_cell](MPI_Comm comm, int /*nparts*/,
+                       std::vector<mesh::CellType> cell_types,
+                       std::vector<std::span<const std::int64_t>> cells)
+             -> graph::AdjacencyList<std::int32_t>
+  {
+    auto cell_im
+        = parent_mesh.topology()->index_map(parent_mesh.topology()->dim());
+
+    std::int32_t num_cells
+        = cells.front().size() / mesh::num_cell_vertices(cell_types.front());
+    std::vector<std::int32_t> destinations(num_cells);
+
+    int rank = dolfinx::MPI::rank(comm);
+    for (std::int32_t i = 0; i < destinations.size(); i++)
+    {
+      bool ghost_parent_cell = parent_cell[i] > cell_im->size_local();
+      if (ghost_parent_cell)
+      {
+        destinations[i]
+            = cell_im->owners()[parent_cell[i] - cell_im->size_local()];
+      }
+      else
+      {
+        destinations[i] = rank;
+      }
+    }
+
+    if (comm == MPI_COMM_NULL)
+    {
+      return graph::regular_adjacency_list(std::move(destinations), 1);
+    }
+
+    auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells);
+    std::vector<std::int32_t> node_disp(MPI::size(comm) + 1, 0);
+    std::int32_t local_size = dual_graph.num_nodes();
+    MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t<std::int32_t>,
+                  node_disp.data() + 1, 1, dolfinx::MPI::mpi_t<std::int32_t>,
+                  comm);
+    std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin());
+    return compute_destination_ranks(comm, dual_graph, node_disp, destinations);
+  };
+}
+
 /// @brief Refine a mesh with markers.
 ///
 /// The refined mesh can be optionally re-partitioned across processes.
@@ -35,20 +103,17 @@ namespace dolfinx::refinement
 /// refined mesh will **not** include ghosts cells (cells connected by
 /// facet to an owned cell) even if the parent mesh is ghosted.
 ///
-/// @warning Passing `nullptr` for `partitioner`, the refined mesh will
-/// **not** have ghosts cells even if the parent mesh is ghosted. The
-/// possibility to not re-partition the refined mesh and include ghost
-/// cells in the refined mesh will be added in a future release.
 ///
 /// @param[in] mesh Input mesh to be refined.
 /// @param[in] edges Indices of the edges that should be split in the
 /// refinement. If not provided (`std::nullopt`), uniform refinement is
 /// performed.
-/// @param[in] partitioner Partitioner to be used to distribute the
-/// refined mesh. If not callable, refined cells will be on the same
+/// @param[in] partitioner (Optional) partitioner to be used to distribute the
+/// refined mesh. If not provided (`std::nullopt`) the partition of the coarse
+/// mesh will be maintained. 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 not selected, an empty list is returned.
+/// cells.
 /// @return New mesh, and optional parent cell indices and parent facet
 /// indices.
 template <std::floating_point T>
@@ -56,9 +121,8 @@ 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,
-       const mesh::CellPartitionFunction& partitioner
-       = mesh::create_cell_partitioner(mesh::GhostMode::none),
-       Option option = Option::none)
+       std::optional<mesh::CellPartitionFunction> partitioner = std::nullopt,
+       Option option = Option::parent_cell)
 {
   auto topology = mesh.topology();
   assert(topology);
@@ -70,9 +134,20 @@ refine(const mesh::Mesh<T>& mesh,
             ? interval::compute_refinement_data(mesh, edges, option)
             : plaza::compute_refinement_data(mesh, edges, option);
 
+  if (!partitioner.has_value())
+  {
+    if (!parent_cell)
+    {
+      throw std::runtime_error(
+          "Identity partitioner relies on parent cell computation");
+    }
+    assert(parent_cell);
+    partitioner = create_identity_partitioner(mesh, parent_cell.value());
+  }
+
   mesh::Mesh<T> mesh1 = mesh::create_mesh(
       mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(),
-      mesh.comm(), new_vertex_coords, xshape, partitioner);
+      mesh.comm(), new_vertex_coords, xshape, *partitioner);
 
   // Report the number of refined cells
   const int D = topology->dim();
@@ -84,4 +159,5 @@ refine(const mesh::Mesh<T>& mesh,
 
   return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)};
 }
+
 } // namespace dolfinx::refinement
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index 55b2f716e83..dab2396b8e9 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -57,6 +57,8 @@ add_executable(
   mesh/refinement/interval.cpp
   mesh/refinement/option.cpp
   mesh/refinement/rectangle.cpp
+  multigrid/transfer_matrix.cpp
+  multigrid/inclusion.cpp
   ${CMAKE_CURRENT_BINARY_DIR}/expr.c
   ${CMAKE_CURRENT_BINARY_DIR}/poisson.c
 )
diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp
new file mode 100644
index 00000000000..d77d14990e3
--- /dev/null
+++ b/cpp/test/multigrid/inclusion.cpp
@@ -0,0 +1,121 @@
+// 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 <concepts>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <vector>
+
+#include <mpi.h>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
+
+#include <basix/cell.h>
+
+#include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/mesh/generation.h>
+#include <dolfinx/multigrid/inclusion.h>
+#include <dolfinx/refinement/refine.h>
+
+using namespace dolfinx;
+using namespace Catch::Matchers;
+
+TEMPLATE_TEST_CASE("Gather global", "[multigrid]", double, float)
+{
+  using T = TestType;
+  MPI_Comm comm = MPI_COMM_WORLD;
+  auto comm_size = dolfinx::MPI::size(comm);
+  auto rank = dolfinx::MPI::rank(comm);
+  std::vector<T> local{static_cast<T>(rank), static_cast<T>(rank + 1)};
+
+  std::vector<T> global = multigrid::gather_global<T>(
+      std::span(local.data(), local.size()), comm_size * 2, comm);
+
+  CHECK(global.size() == static_cast<std::size_t>(2 * comm_size));
+  for (int i = 0; i < comm_size; i++)
+  {
+    CHECK(global[2 * i] == Catch::Approx(i));
+    CHECK(global[2 * i + 1] == Catch::Approx(i + 1));
+  }
+}
+
+template <std::floating_point T>
+void CHECK_inclusion_map(const dolfinx::mesh::Mesh<T>& from,
+                         const dolfinx::mesh::Mesh<T>& to,
+                         const std::vector<std::int64_t>& map)
+{
+  const common::IndexMap& im_from = *from.topology()->index_map(0);
+  const common::IndexMap& im_to = *to.topology()->index_map(0);
+
+  std::vector<T> global_x_to = multigrid::gather_global(
+      to.geometry().x().subspan(0, im_to.size_local() * 3),
+      im_to.size_global() * 3, to.comm());
+
+  REQUIRE(static_cast<std::int64_t>(map.size())
+          == im_from.size_local() + im_from.num_ghosts());
+  for (std::int64_t i = 0; i < static_cast<std::int64_t>(map.size()); i++)
+  {
+    CHECK(std::abs(from.geometry().x()[3 * i] - global_x_to[3 * map[i]])
+          < std::numeric_limits<T>::epsilon());
+    CHECK(std::abs(from.geometry().x()[3 * i + 1] - global_x_to[3 * map[i] + 1])
+          < std::numeric_limits<T>::epsilon());
+    CHECK(std::abs(from.geometry().x()[3 * i + 2] - global_x_to[3 * map[i] + 2])
+          < std::numeric_limits<T>::epsilon());
+  }
+}
+
+/// Performs one uniform refinement and checks the inclusion map between coarse
+/// and fine mesh.
+template <std::floating_point T>
+void TEST_inclusion(dolfinx::mesh::Mesh<T>&& mesh_coarse)
+{
+  mesh_coarse.topology()->create_entities(1);
+
+  auto [mesh_fine, parent_cell, parent_facet]
+      = refinement::refine(mesh_coarse, std::nullopt);
+  mesh_fine.topology()->create_connectivity(1, 0);
+  mesh_fine.topology()->create_connectivity(0, 1);
+  std::vector<std::int64_t> inclusion_map
+      = multigrid::inclusion_mapping(mesh_coarse, mesh_fine, true);
+
+  CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map);
+}
+
+TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double,
+                   float)
+{
+  for (auto n : {10})
+  {
+    TEST_inclusion(dolfinx::mesh::create_interval<TestType>(MPI_COMM_WORLD, n,
+                                                            {0.0, 1.0}));
+  }
+}
+
+TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double,
+                   float)
+{
+  for (auto n : {5})
+  {
+    TEST_inclusion(dolfinx::mesh::create_rectangle<TestType>(
+        MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {n, n}, mesh::CellType::triangle));
+  }
+}
+
+TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double,
+                   float)
+{
+  for (auto n : {5})
+  {
+    TEST_inclusion(dolfinx::mesh::create_box<TestType>(
+        MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {n, n, n},
+        mesh::CellType::tetrahedron));
+  }
+}
diff --git a/cpp/test/multigrid/transfer_matrix.cpp b/cpp/test/multigrid/transfer_matrix.cpp
new file mode 100644
index 00000000000..d5066d98fb8
--- /dev/null
+++ b/cpp/test/multigrid/transfer_matrix.cpp
@@ -0,0 +1,444 @@
+// 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 <basix/cell.h>
+#include <concepts>
+#include <cstdint>
+#include <memory>
+#include <optional>
+#include <vector>
+
+#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/common/IndexMap.h>
+#include <dolfinx/common/MPI.h>
+#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/multigrid/inclusion.h>
+#include <dolfinx/multigrid/transfer_matrix.h>
+#include <dolfinx/refinement/plaza.h>
+#include <dolfinx/refinement/refine.h>
+
+using namespace dolfinx;
+using namespace Catch::Matchers;
+
+template <std::floating_point T>
+constexpr auto weight
+    = [](std::int32_t d) -> T { return d == 0 ? T(1) : T(.5); };
+
+template <typename T = double>
+std::shared_ptr<const fem::FiniteElement<T>>
+P1_element(basix::cell::type cell_type)
+{
+  auto element
+      = basix::create_element<T>(basix::element::family::P, cell_type, 1,
+                                 basix::element::lagrange_variant::unset,
+                                 basix::element::dpc_variant::unset, false);
+  return std::make_shared<const fem::FiniteElement<T>>(element);
+}
+
+TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, 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);
+
+  auto element = P1_element<T>(basix::cell::type::interval);
+  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> inclusion_map
+      = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine);
+
+  CHECK_THAT(inclusion_map, RangeEquals(std::vector<std::int64_t>{0, 2, 4}));
+
+  la::SparsityPattern sp
+      = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  multigrid::assemble_transfer_matrix<T>(transfer_matrix.mat_set_values(),
+                                         *V_coarse, *V_fine, inclusion_map,
+                                         weight<T>);
+
+  // clang-format off
+  std::vector<T> expected{1.0, 0.5, 0.0, 0.0, 0.0,
+                          0.0, 0.5, 1.0, 0.5, 0.0,
+                          0.0, 0.0, 0.0, 0.5, 1.0};
+  // clang-format on
+  CHECK_THAT(transfer_matrix.to_dense(),
+             RangeEquals(expected,
+                         [](auto a, auto b) {
+                           return std::abs(a - b)
+                                  <= std::numeric_limits<T>::epsilon();
+                         }));
+}
+
+// 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::none, part));
+
+//   // TODO: see https://github.com/FEniCS/dolfinx/issues/3358
+//   //   auto part_refine
+//   //       = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
+//   mesh::CellPartitionFunction part_refine
+//       = [&](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 (dolfinx::MPI::size(comm) == 1)
+//       data = {{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}};
+//     else if (dolfinx::MPI::size(comm) == 2)
+//     {
+//       if (rank == 0)
+//         data = {{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0, 1}, {1, 0}};
+//       else
+//         data = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
+//     }
+//     else if (dolfinx::MPI::size(comm) == 3)
+//     {
+//       if (rank == 0)
+//         data = {{2, 0}, {0, 2}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}};
+//       else if (rank == 1)
+//         data = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1, 2}};
+//       else
+//         data = {{2, 1}, {2}, {2}, {2}, {2}, {2}, {2}, {2}};
+//     }
+//     else
+//       FAIL("Test only supports <= 3 processes");
+//     return graph::AdjacencyList<std::int32_t>(std::move(data));
+//   };
+
+//   auto [mesh_fine, parent_cell, parent_facet] = refinement::refine(
+//       *mesh_coarse, std::nullopt, part_refine, refinement::Option::none);
+
+//   auto element = P1_element(basix::cell::type::interval);
+
+//   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> inclusion_map;
+//   switch (comm_size)
+//   {
+//   case 1:
+//     inclusion_map = {0, 2, 4, 6, 8};
+//     break;
+//   case 2:
+//     inclusion_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18};
+//     break;
+//   case 3:
+//     inclusion_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 27, 25, 23, 20};
+//     break;
+//   }
+
+//   CHECK_THAT(multigrid::inclusion_mapping(*mesh_coarse, mesh_fine),
+//              RangeEquals(inclusion_map));
+
+//   la::SparsityPattern sp
+//       = multigrid::create_sparsity_pattern(*V_coarse, *V_fine,
+//       inclusion_map);
+
+//   la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+//   multigrid::assemble_transfer_matrix<T>(transfer_matrix.mat_set_values(),
+//                                         *V_coarse, *V_fine, inclusion_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.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+//         0.0, 0.5, 1.0, 0.5,
+//         /* 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.5, 1.0, 0.5, 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.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,
+//         /* 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, 1.0, 0.5, 0.5, 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.0, 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)
+//                                   <= std::numeric_limits<T>::epsilon();
+//                          }));
+// }
+
+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);
+
+  auto element = P1_element(basix::cell::type::triangle);
+
+  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> inclusion_map
+      = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine);
+  CHECK_THAT(inclusion_map, RangeEquals(std::vector<std::int64_t>{4, 1, 5, 8}));
+
+  la::SparsityPattern sp
+      = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  multigrid::assemble_transfer_matrix<T>(transfer_matrix.mat_set_values(),
+                                         *V_coarse, *V_fine, inclusion_map,
+                                         weight<T>);
+  transfer_matrix.scatter_rev();
+  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)
+                                  <= std::numeric_limits<T>::epsilon();
+                         }));
+}
+
+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);
+
+  auto element = P1_element(basix::cell::type::tetrahedron);
+
+  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> inclusion_map
+      = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine);
+  CHECK_THAT(inclusion_map, RangeEquals(std::vector<std::int64_t>{
+                                0, 6, 15, 25, 17, 9, 11, 22}));
+
+  la::SparsityPattern sp
+      = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  multigrid::assemble_transfer_matrix<T>(transfer_matrix.mat_set_values(),
+                                         *V_coarse, *V_fine, inclusion_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)
+                                  <= std::numeric_limits<T>::epsilon();
+                         }));
+}
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 4b92f1bb14e..0f8bbd6ad7b 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -45,6 +45,7 @@ nanobind_add_module(
   dolfinx/wrappers/la.cpp
   dolfinx/wrappers/log.cpp
   dolfinx/wrappers/mesh.cpp
+  dolfinx/wrappers/multigrid.cpp
   dolfinx/wrappers/petsc.cpp
   dolfinx/wrappers/refinement.cpp
 )
diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py
index 46385f0d345..cc11abd7db5 100644
--- a/python/dolfinx/mesh.py
+++ b/python/dolfinx/mesh.py
@@ -518,7 +518,7 @@ def transfer_meshtag(
 def refine(
     msh: Mesh,
     edges: typing.Optional[np.ndarray] = None,
-    partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none),
+    partitioner: typing.Optional[typing.Callable] = None,
     option: RefinementOption = RefinementOption.none,
 ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
     """Refine a mesh.
diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py
new file mode 100644
index 00000000000..b045466e427
--- /dev/null
+++ b/python/dolfinx/multigrid.py
@@ -0,0 +1,44 @@
+# 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
+
+from typing import Callable
+
+import numpy as np
+from numpy.typing import NDArray
+
+from dolfinx.cpp.la import SparsityPattern
+from dolfinx.cpp.multigrid import assemble_transfer_matrix as _assemble_transfer_matrix
+from dolfinx.cpp.multigrid import create_sparsity_pattern as _create_sparsity_pattern
+from dolfinx.cpp.multigrid import inclusion_mapping as _inclusion_mapping
+from dolfinx.fem import FunctionSpace
+from dolfinx.la import MatrixCSR
+from dolfinx.mesh import Mesh
+
+__all__ = ["assemble_transfer_matrix", "create_sparsity_pattern", "inclusion_mapping"]
+
+
+def inclusion_mapping(
+    mesh_from: Mesh, mesh_to: Mesh, allow_all_to_all: bool = False
+) -> NDArray[np.int64]:
+    return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object, allow_all_to_all)
+
+
+def assemble_transfer_matrix(
+    T: MatrixCSR,
+    V_from: FunctionSpace,
+    V_to: FunctionSpace,
+    inclusion_map: NDArray[np.int64],
+    weight: Callable[[int], float],
+):
+    _assemble_transfer_matrix(
+        T._cpp_object, V_from._cpp_object, V_to._cpp_object, inclusion_map, weight
+    )
+
+
+def create_sparsity_pattern(
+    V_from: FunctionSpace, V_to: FunctionSpace, inclusion_map: NDArray[np.int64]
+) -> SparsityPattern:
+    return _create_sparsity_pattern(V_from._cpp_object, V_to._cpp_object, inclusion_map)
diff --git a/python/dolfinx/wrappers/dolfinx.cpp b/python/dolfinx/wrappers/dolfinx.cpp
index 9c0864f9ef1..21bcec93f4f 100644
--- a/python/dolfinx/wrappers/dolfinx.cpp
+++ b/python/dolfinx/wrappers/dolfinx.cpp
@@ -25,6 +25,8 @@ void la(nb::module_& m);
 void mesh(nb::module_& m);
 void nls(nb::module_& m);
 void refinement(nb::module_& m);
+void multigrid(nb::module_& m);
+
 } // namespace dolfinx_wrappers
 
 NB_MODULE(cpp, m)
@@ -74,6 +76,10 @@ NB_MODULE(cpp, m)
   nb::module_ refinement = m.def_submodule("refinement", "Refinement module");
   dolfinx_wrappers::refinement(refinement);
 
+  // Create multigrid submodule
+  nb::module_ multigrid = m.def_submodule("multigrid", "Multigrid module");
+  dolfinx_wrappers::multigrid(multigrid);
+
 #if defined(HAS_PETSC) && defined(HAS_PETSC4PY)
   // PETSc-specific wrappers
   nb::module_ nls = m.def_submodule("nls", "Nonlinear solver module");
diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp
new file mode 100644
index 00000000000..e0b22fa4118
--- /dev/null
+++ b/python/dolfinx/wrappers/multigrid.cpp
@@ -0,0 +1,71 @@
+// 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 <cstdint>
+#include <vector>
+
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/function.h>
+#include <nanobind/stl/vector.h>
+
+#include <dolfinx/fem/FunctionSpace.h>
+#include <dolfinx/la/MatrixCSR.h>
+#include <dolfinx/multigrid/inclusion.h>
+#include <dolfinx/multigrid/transfer_matrix.h>
+
+#include "dolfinx_wrappers/array.h"
+
+namespace nb = nanobind;
+
+namespace dolfinx_wrappers
+{
+
+void multigrid(nb::module_& m)
+{
+  m.def(
+      "inclusion_mapping",
+      [](const dolfinx::mesh::Mesh<double>& mesh_from,
+         const dolfinx::mesh::Mesh<double>& mesh_to, bool allow_all_to_all)
+      {
+        std::vector<std::int64_t> map
+            = dolfinx::multigrid::inclusion_mapping<double>(mesh_from, mesh_to,
+                                                            allow_all_to_all);
+        return dolfinx_wrappers::as_nbarray(std::move(map));
+      },
+      nb::arg("mesh_from"), nb::arg("mesh_to"), nb::arg("allow_all_to_all"),
+      "Computes inclusion mapping between two meshes");
+
+  m.def(
+      "create_sparsity_pattern",
+      [](const dolfinx::fem::FunctionSpace<double>& V_from,
+         const dolfinx::fem::FunctionSpace<double>& V_to,
+         nb::ndarray<const std::int64_t, nb::numpy>& inclusion_map)
+      {
+        // TODO: cahnge to accepting range;
+        auto vec = std::vector(inclusion_map.data(),
+                               inclusion_map.data() + inclusion_map.size());
+        return dolfinx::multigrid::create_sparsity_pattern<double>(V_from, V_to,
+                                                                   vec);
+      },
+      nb::arg("V_from"), nb::arg("V_to"), nb::arg("inclusion_map"));
+
+  m.def(
+      "assemble_transfer_matrix",
+      [](dolfinx::la::MatrixCSR<double>& A,
+         const dolfinx::fem::FunctionSpace<double>& V_from,
+         const dolfinx::fem::FunctionSpace<double>& V_to,
+         const std::vector<std::int64_t>& inclusion_map,
+         std::function<double(std::int32_t)> weight)
+      {
+        dolfinx::multigrid::assemble_transfer_matrix(
+            A.mat_set_values(), V_from, V_to, inclusion_map, weight);
+      },
+      nb::arg("A"), nb::arg("V_from"), nb::arg("V_to"),
+      nb::arg("inclusion_map"), nb::arg("weight"),
+      "Assembles a transfer matrix between two function spaces.");
+}
+
+} // namespace dolfinx_wrappers
diff --git a/python/test/unit/multigrid/test_inclusion.py b/python/test/unit/multigrid/test_inclusion.py
new file mode 100644
index 00000000000..06fee848579
--- /dev/null
+++ b/python/test/unit/multigrid/test_inclusion.py
@@ -0,0 +1,48 @@
+# 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
+
+from mpi4py import MPI
+
+import pytest
+
+from dolfinx.mesh import GhostMode, create_interval, create_unit_cube, create_unit_square, refine
+from dolfinx.multigrid import inclusion_mapping
+
+
+@pytest.mark.parametrize(
+    "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet]
+)
+def test_1d(ghost_mode):
+    mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode)
+    mesh_fine, _, _ = refine(mesh)
+    inclusion_mapping(mesh, mesh_fine, True)
+    # TODO: extend with future operations on inclusion mappings
+
+
+@pytest.mark.parametrize(
+    "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet]
+)
+def test_2d(ghost_mode):
+    mesh = create_unit_square(MPI.COMM_WORLD, 5, 5, ghost_mode=ghost_mode)
+    mesh.topology.create_entities(1)
+    mesh_fine, _, _ = refine(mesh)
+    inclusion_mapping(mesh, mesh_fine, True)
+    # TODO: extend with future operations on inclusion mappings
+
+
+@pytest.mark.parametrize(
+    "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet]
+)
+def test_3d(ghost_mode):
+    mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5, ghost_mode=ghost_mode)
+    mesh.topology.create_entities(1)
+    mesh_fine, _, _ = refine(mesh)
+    inclusion_mapping(mesh, mesh_fine, True)
+    # TODO: extend with future operations on inclusion mappings
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])
diff --git a/python/test/unit/multigrid/test_transfer.py b/python/test/unit/multigrid/test_transfer.py
new file mode 100644
index 00000000000..3bfbfcc66be
--- /dev/null
+++ b/python/test/unit/multigrid/test_transfer.py
@@ -0,0 +1,58 @@
+# 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
+
+from mpi4py import MPI
+
+import pytest
+
+from dolfinx.fem import functionspace
+from dolfinx.la import matrix_csr
+from dolfinx.mesh import GhostMode, create_interval, refine
+from dolfinx.multigrid import assemble_transfer_matrix, create_sparsity_pattern, inclusion_mapping
+
+
+@pytest.mark.parametrize(
+    "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet]
+)
+def test_1d(ghost_mode):
+    mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode)
+    mesh_fine, _, _ = refine(mesh)
+    inclusion_map = inclusion_mapping(mesh, mesh_fine, True)
+    V = functionspace(mesh, ("Lagrange", 1, (1,)))
+    V_fine = functionspace(mesh_fine, ("Lagrange", 1, (1,)))
+    assert V.element == V_fine.element
+    V_fine.mesh.topology.create_connectivity(1, 0)
+    V_fine.mesh.topology.create_connectivity(0, 1)
+    sp = create_sparsity_pattern(V, V_fine, inclusion_map)
+    T = matrix_csr(sp)
+    assemble_transfer_matrix(T, V, V_fine, inclusion_map, lambda i: 1.0 if i == 0 else 0.5)
+    # continue with assembly of matrix
+
+
+@pytest.mark.parametrize("degree", [2, 3, 4])
+@pytest.mark.parametrize(
+    "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet]
+)
+def test_1d_higher_order(degree, ghost_mode):
+    mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode)
+    mesh_fine, _, _ = refine(mesh)
+
+    # this is a strictly geometric operation, and thus should pass
+    inclusion_map = inclusion_mapping(mesh, mesh_fine, True)
+
+    V = functionspace(mesh, ("Lagrange", degree, (1,)))
+    V_fine = functionspace(mesh_fine, ("Lagrange", degree, (1,)))
+
+    V_fine.mesh.topology.create_connectivity(1, 0)
+    V_fine.mesh.topology.create_connectivity(0, 1)
+
+    # Check not supported throws
+    with pytest.raises(Exception):
+        create_sparsity_pattern(V, V_fine, inclusion_map)
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])