From f20421f633e344d0963dfa3ba82579cb28925c01 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Sun, 20 Oct 2024 22:39:01 +0200
Subject: [PATCH 1/3] Introduce transfer matrix operation

---
 cpp/dolfinx/CMakeLists.txt                 |   1 +
 cpp/dolfinx/transfer/CMakeLists.txt        |   4 +
 cpp/dolfinx/transfer/transfer_matrix.h     | 236 +++++++++++++
 cpp/test/CMakeLists.txt                    |   1 +
 cpp/test/transfer/transfer_matrix.cpp      | 377 +++++++++++++++++++++
 python/CMakeLists.txt                      |   1 +
 python/dolfinx/transfer.py                 |  23 ++
 python/dolfinx/wrappers/dolfinx.cpp        |   6 +
 python/dolfinx/wrappers/transfer.cpp       |  67 ++++
 python/test/unit/transfer/test_transfer.py |  22 ++
 10 files changed, 738 insertions(+)
 create mode 100644 cpp/dolfinx/transfer/CMakeLists.txt
 create mode 100644 cpp/dolfinx/transfer/transfer_matrix.h
 create mode 100644 cpp/test/transfer/transfer_matrix.cpp
 create mode 100644 python/dolfinx/transfer.py
 create mode 100644 python/dolfinx/wrappers/transfer.cpp
 create mode 100644 python/test/unit/transfer/test_transfer.py

diff --git a/cpp/dolfinx/CMakeLists.txt b/cpp/dolfinx/CMakeLists.txt
index 7adaaeb73cf..30c60f78576 100644
--- a/cpp/dolfinx/CMakeLists.txt
+++ b/cpp/dolfinx/CMakeLists.txt
@@ -19,6 +19,7 @@ set(DOLFINX_DIRS
     mesh
     nls
     refinement
+    transfer
 )
 
 # Add source to dolfinx target, and get sets of header files
diff --git a/cpp/dolfinx/transfer/CMakeLists.txt b/cpp/dolfinx/transfer/CMakeLists.txt
new file mode 100644
index 00000000000..68bcddd738e
--- /dev/null
+++ b/cpp/dolfinx/transfer/CMakeLists.txt
@@ -0,0 +1,4 @@
+set(HEADERS_transfer
+    ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h
+    PARENT_SCOPE
+)
diff --git a/cpp/dolfinx/transfer/transfer_matrix.h b/cpp/dolfinx/transfer/transfer_matrix.h
new file mode 100644
index 00000000000..add92cdc9d1
--- /dev/null
+++ b/cpp/dolfinx/transfer/transfer_matrix.h
@@ -0,0 +1,236 @@
+// 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 <vector>
+
+#include "dolfinx/common/IndexMap.h"
+#include "dolfinx/fem/FunctionSpace.h"
+#include "dolfinx/la/SparsityPattern.h"
+#include "dolfinx/la/utils.h"
+#include "dolfinx/mesh/Mesh.h"
+
+namespace dolfinx::transfer
+{
+
+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)
+{
+
+  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_global(), -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(); 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(); }))
+      {
+        map[to_global_from(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;
+  }
+
+  // map holds at this point for every original local index the corresponding
+  // mapped global index. All other entries are still -1, but are available on
+  // other processes.
+  std::vector<std::int64_t> result(map.size(), -1);
+  MPI_Allreduce(map.data(), result.data(), map.size(),
+                dolfinx::MPI::mpi_type<std::int64_t>(), MPI_MAX,
+                mesh_from.comm());
+
+  assert(std::ranges::all_of(result, [](auto e) { return e >= 0; }));
+  return result;
+}
+
+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)
+{
+  // TODO: P1 and which elements supported?
+  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)
+{
+  // TODO: P1 and which elements supported?
+  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);
+
+  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<double>{weight(distance)});
+      }
+    }
+  }
+}
+
+} // namespace dolfinx::transfer
\ No newline at end of file
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index 4a3dd8f3817..64315eada0e 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -47,6 +47,7 @@ add_executable(
   mesh/refinement/interval.cpp
   mesh/refinement/option.cpp
   mesh/refinement/rectangle.cpp
+  transfer/transfer_matrix.cpp
   ${CMAKE_CURRENT_BINARY_DIR}/poisson.c
 )
 target_link_libraries(unittests PRIVATE Catch2::Catch2WithMain dolfinx)
diff --git a/cpp/test/transfer/transfer_matrix.cpp b/cpp/test/transfer/transfer_matrix.cpp
new file mode 100644
index 00000000000..3487090d01b
--- /dev/null
+++ b/cpp/test/transfer/transfer_matrix.cpp
@@ -0,0 +1,377 @@
+// 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 <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/refinement/plaza.h>
+#include <dolfinx/refinement/refine.h>
+#include <dolfinx/transfer/transfer_matrix.h>
+
+using namespace dolfinx;
+using namespace Catch::Matchers;
+
+namespace
+{
+template <typename T>
+constexpr auto EPS = std::numeric_limits<T>::epsilon();
+} // namespace
+
+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);
+
+  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> inclusion_map
+      = transfer::inclusion_mapping(*mesh_coarse, mesh_fine);
+
+  CHECK_THAT(inclusion_map, RangeEquals(std::vector<std::int64_t>{0, 2, 4}));
+
+  la::SparsityPattern sp
+      = transfer::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  transfer::assemble_transfer_matrix<T>(transfer_matrix.mat_set_values(),
+                                        *V_coarse, *V_fine, inclusion_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));
+
+  // 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 = 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> 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(transfer::inclusion_mapping(*mesh_coarse, mesh_fine),
+             RangeEquals(inclusion_map));
+
+  la::SparsityPattern sp
+      = transfer::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  transfer::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) <= 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);
+
+  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> inclusion_map
+      = transfer::inclusion_mapping(*mesh_coarse, mesh_fine);
+  CHECK_THAT(inclusion_map, RangeEquals(std::vector<std::int64_t>{4, 1, 5, 8}));
+
+  la::SparsityPattern sp
+      = transfer::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  transfer::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) <= 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);
+
+  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> inclusion_map
+      = transfer::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
+      = transfer::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map);
+  la::MatrixCSR<T> transfer_matrix(std::move(sp), la::BlockMode::compact);
+  transfer::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) <= EPS<T>; }));
+}
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 0f440260214..dc13b4dc957 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -46,6 +46,7 @@ nanobind_add_module(
   dolfinx/wrappers/mesh.cpp
   dolfinx/wrappers/petsc.cpp
   dolfinx/wrappers/refinement.cpp
+  dolfinx/wrappers/transfer.cpp
 )
 
 set(CMAKE_CXX_EXTENSIONS OFF)
diff --git a/python/dolfinx/transfer.py b/python/dolfinx/transfer.py
new file mode 100644
index 00000000000..4e2a7670f3b
--- /dev/null
+++ b/python/dolfinx/transfer.py
@@ -0,0 +1,23 @@
+import numpy as np
+from numpy.typing import NDArray
+
+from dolfinx.cpp.la import SparsityPattern
+from dolfinx.cpp.transfer import (
+    assemble_transfer_matrix,
+)
+from dolfinx.cpp.transfer import create_sparsity_pattern as _create_sparsity_pattern
+from dolfinx.cpp.transfer import inclusion_mapping as _inclusion_mapping
+from dolfinx.fem import FunctionSpace
+from dolfinx.mesh import Mesh
+
+__all__ = ["inclusion_mapping", "create_sparsity_pattern", "assemble_transfer_matrix"]
+
+
+def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]:
+    return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object)
+
+
+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..9d2ca3e716d 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 transfer(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 transfer submodule
+  nb::module_ transfer = m.def_submodule("transfer", "Transfer module");
+  dolfinx_wrappers::transfer(transfer);
+
 #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/transfer.cpp b/python/dolfinx/wrappers/transfer.cpp
new file mode 100644
index 00000000000..16ffc4cbe5a
--- /dev/null
+++ b/python/dolfinx/wrappers/transfer.cpp
@@ -0,0 +1,67 @@
+// 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/vector.h>
+
+#include <dolfinx/la/MatrixCSR.h>
+#include <dolfinx/transfer/transfer_matrix.h>
+
+#include "array.h"
+
+namespace nb = nanobind;
+
+namespace dolfinx_wrappers
+{
+
+void transfer(nb::module_& m)
+{
+  m.def(
+      "inclusion_mapping",
+      [](const dolfinx::mesh::Mesh<double>& mesh_from,
+         const dolfinx::mesh::Mesh<double>& mesh_to)
+      {
+        std::vector<std::int64_t> map
+            = dolfinx::transfer::inclusion_mapping<double>(mesh_from, mesh_to);
+        return dolfinx_wrappers::as_nbarray(map);
+      },
+      nb::arg("mesh_from"), nb::arg("mesh_to"),
+      "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::transfer::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::transfer::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/transfer/test_transfer.py b/python/test/unit/transfer/test_transfer.py
new file mode 100644
index 00000000000..93c260442da
--- /dev/null
+++ b/python/test/unit/transfer/test_transfer.py
@@ -0,0 +1,22 @@
+from mpi4py import MPI
+
+import pytest
+
+from dolfinx.fem import functionspace
+from dolfinx.mesh import GhostMode, create_interval, refine
+from dolfinx.transfer import 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)
+    V = functionspace(mesh, ("Lagrange", 1, (1,)))
+    V_fine = functionspace(mesh_fine, ("Lagrange", 1, (1,)))
+    V_fine.mesh.topology.create_connectivity(1, 0)
+    V_fine.mesh.topology.create_connectivity(0, 1)
+    create_sparsity_pattern(V, V_fine, inclusion_map)
+    # continue with assembly of matrix

From 3b4e5d08cc52877994582875ffe4fee96273fd06 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Sun, 20 Oct 2024 22:51:44 +0200
Subject: [PATCH 2/3] Add missing headers

---
 python/dolfinx/transfer.py                 | 6 ++++++
 python/test/unit/transfer/test_transfer.py | 6 ++++++
 2 files changed, 12 insertions(+)

diff --git a/python/dolfinx/transfer.py b/python/dolfinx/transfer.py
index 4e2a7670f3b..b98781ac5ed 100644
--- a/python/dolfinx/transfer.py
+++ b/python/dolfinx/transfer.py
@@ -1,3 +1,9 @@
+# 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
+
 import numpy as np
 from numpy.typing import NDArray
 
diff --git a/python/test/unit/transfer/test_transfer.py b/python/test/unit/transfer/test_transfer.py
index 93c260442da..618d3be0a5f 100644
--- a/python/test/unit/transfer/test_transfer.py
+++ b/python/test/unit/transfer/test_transfer.py
@@ -1,3 +1,9 @@
+# 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 9d7b749b6b29054e6ddd477ae29885e7ed6e9e84 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 25 Jun 2024 10:23:50 +0200
Subject: [PATCH 3/3] Add geometric twogrid demo

---
 cpp/demo/CMakeLists.txt                   |   1 +
 cpp/demo/geometric-twogrid/CMakeLists.txt |  67 ++++++
 cpp/demo/geometric-twogrid/main.cpp       | 242 ++++++++++++++++++++++
 cpp/demo/geometric-twogrid/poisson.py     |  44 ++++
 4 files changed, 354 insertions(+)
 create mode 100644 cpp/demo/geometric-twogrid/CMakeLists.txt
 create mode 100644 cpp/demo/geometric-twogrid/main.cpp
 create mode 100644 cpp/demo/geometric-twogrid/poisson.py

diff --git a/cpp/demo/CMakeLists.txt b/cpp/demo/CMakeLists.txt
index 9111ed7fd78..fea72b664d7 100644
--- a/cpp/demo/CMakeLists.txt
+++ b/cpp/demo/CMakeLists.txt
@@ -19,6 +19,7 @@ endmacro(add_demo_subdirectory)
 add_demo_subdirectory(biharmonic)
 add_demo_subdirectory(codim_0_assembly)
 add_demo_subdirectory(custom_kernel)
+add_demo_subdirectory(geometric-multigrid)
 add_demo_subdirectory(hyperelasticity)
 add_demo_subdirectory(interpolation-io)
 add_demo_subdirectory(interpolation_different_meshes)
diff --git a/cpp/demo/geometric-twogrid/CMakeLists.txt b/cpp/demo/geometric-twogrid/CMakeLists.txt
new file mode 100644
index 00000000000..714bef6a828
--- /dev/null
+++ b/cpp/demo/geometric-twogrid/CMakeLists.txt
@@ -0,0 +1,67 @@
+# This file was generated by running
+#
+# python cmake/scripts/generate-cmakefiles from dolfinx/cpp
+#
+cmake_minimum_required(VERSION 3.19)
+
+set(PROJECT_NAME demo_geometric-twogrid)
+project(${PROJECT_NAME} LANGUAGES C CXX)
+
+# Set C++20 standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+if(NOT TARGET dolfinx)
+  find_package(DOLFINX REQUIRED)
+endif()
+
+include(CheckSymbolExists)
+set(CMAKE_REQUIRED_INCLUDES ${PETSC_INCLUDE_DIRS})
+check_symbol_exists(PETSC_USE_COMPLEX petscsystypes.h PETSC_SCALAR_COMPLEX)
+check_symbol_exists(PETSC_USE_REAL_DOUBLE petscsystypes.h PETSC_REAL_DOUBLE)
+
+# Add target to compile UFL files
+if(PETSC_SCALAR_COMPLEX EQUAL 1)
+  if(PETSC_REAL_DOUBLE EQUAL 1)
+    set(SCALAR_TYPE "--scalar_type=complex128")
+  else()
+    set(SCALAR_TYPE "--scalar_type=complex64")
+  endif()
+else()
+  if(PETSC_REAL_DOUBLE EQUAL 1)
+    set(SCALAR_TYPE "--scalar_type=float64")
+  else()
+    set(SCALAR_TYPE "--scalar_type=float32")
+  endif()
+endif()
+add_custom_command(
+  OUTPUT poisson.c
+  COMMAND ffcx ${CMAKE_CURRENT_SOURCE_DIR}/poisson.py ${SCALAR_TYPE}
+  VERBATIM
+  DEPENDS poisson.py
+  COMMENT "Compile poisson.py using FFCx"
+)
+
+set(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/poisson.c)
+target_link_libraries(${PROJECT_NAME} dolfinx)
+
+# Do not throw error for 'multi-line comments' (these are typical in rst which
+# includes LaTeX)
+include(CheckCXXCompilerFlag)
+check_cxx_compiler_flag("-Wno-comment" HAVE_NO_MULTLINE)
+set_source_files_properties(
+  main.cpp
+  PROPERTIES
+    COMPILE_FLAGS
+    "$<$<BOOL:${HAVE_NO_MULTLINE}>:-Wno-comment -Wall -Wextra -pedantic -Werror>"
+)
+
+# Test targets (used by DOLFINx testing system)
+set(TEST_PARAMETERS2 -np 2 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}")
+set(TEST_PARAMETERS3 -np 3 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}")
+add_test(NAME ${PROJECT_NAME}_mpi_2 COMMAND "mpirun" ${TEST_PARAMETERS2})
+add_test(NAME ${PROJECT_NAME}_mpi_3 COMMAND "mpirun" ${TEST_PARAMETERS3})
+add_test(NAME ${PROJECT_NAME}_serial COMMAND ${PROJECT_NAME})
diff --git a/cpp/demo/geometric-twogrid/main.cpp b/cpp/demo/geometric-twogrid/main.cpp
new file mode 100644
index 00000000000..7bd26721f82
--- /dev/null
+++ b/cpp/demo/geometric-twogrid/main.cpp
@@ -0,0 +1,242 @@
+// Copyright (C) 2024 Paul Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+/// @cond TODO
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <numbers>
+
+#include <petscksp.h>
+#include <petsclog.h>
+#include <petscmat.h>
+#include <petscpc.h>
+#include <petscpctypes.h>
+#include <petscsys.h>
+#include <petscsystypes.h>
+#include <petscvec.h>
+#include <petscviewer.h>
+
+#include <mpi.h>
+
+#include <basix/finite-element.h>
+
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/fem/DirichletBC.h>
+#include <dolfinx/fem/FunctionSpace.h>
+#include <dolfinx/fem/petsc.h>
+#include <dolfinx/fem/utils.h>
+#include <dolfinx/io/ADIOS2Writers.h>
+#include <dolfinx/io/VTKFile.h>
+#include <dolfinx/io/XDMFFile.h>
+#include <dolfinx/la/SparsityPattern.h>
+#include <dolfinx/la/petsc.h>
+#include <dolfinx/mesh/cell_types.h>
+#include <dolfinx/mesh/generation.h>
+#include <dolfinx/mesh/utils.h>
+#include <dolfinx/refinement/refine.h>
+#include <dolfinx/transfer/transfer_matrix.h>
+
+#include "poisson.h"
+
+using namespace dolfinx;
+using T = PetscScalar;
+using U = typename dolfinx::scalar_value_type_t<T>;
+
+struct PetscEnv
+{
+  PetscEnv(int argc, char** argv) { PetscInitialize(&argc, &argv, NULL, NULL); }
+
+  ~PetscEnv() { PetscFinalize(); }
+};
+
+// recommended to run with
+// ./demo_geometric-multigrid -pc_mg_log -all_ksp_monitor -ksp_converged_reason
+// -ksp_view
+int main(int argc, char** argv)
+{
+  PetscEnv petscEnv(argc, argv);
+  // PetscLogDefaultBegin();
+
+  auto element = basix::create_element<U>(
+      basix::element::family::P, basix::cell::type::tetrahedron, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+  auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_vertex);
+
+  auto mesh_coarse = std::make_shared<mesh::Mesh<U>>(dolfinx::mesh::create_box(
+      MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {10, 10, 10},
+      mesh::CellType::tetrahedron, part));
+
+  auto V_coarse = std::make_shared<fem::FunctionSpace<U>>(
+      fem::create_functionspace<U>(mesh_coarse, element, {}));
+
+  // refinement routine requires edges to be initialized
+  mesh_coarse->topology()->create_entities(1);
+  auto [mesh, parent_cells, parent_facets]
+      = dolfinx::refinement::refine(*mesh_coarse, std::nullopt);
+
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      std::make_shared<mesh::Mesh<U>>(mesh), element, {}));
+
+  auto f_ana = [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
+  {
+    std::vector<T> f;
+    for (std::size_t p = 0; p < x.extent(1); ++p)
+    {
+      auto x0 = x(0, p);
+      auto x1 = x(1, p);
+      auto x2 = x(2, p);
+      auto pi = std::numbers::pi;
+      f.push_back(2 * pi * pi * std::sin(pi * x0) * std::sin(pi * x1)
+                  * std::sin(pi * x2));
+    }
+    return {f, {f.size()}};
+  };
+  auto f = std::make_shared<fem::Function<T>>(V);
+  f->interpolate(f_ana);
+
+  {
+    io::VTKFile file(MPI_COMM_WORLD, "f.pvd", "w");
+    file.write<T>({*f}, 0.0);
+  }
+
+  auto a = std::make_shared<fem::Form<T>>(
+      fem::create_form<T>(*form_poisson_a, {V, V}, {}, {}, {}));
+  auto a_coarse = std::make_shared<fem::Form<T>>(
+      fem::create_form<T>(*form_poisson_a, {V_coarse, V_coarse}, {}, {}, {}));
+  auto L = std::make_shared<fem::Form<T>>(
+      fem::create_form<T>(*form_poisson_L, {V}, {{"f", f}}, {}, {}));
+  la::petsc::Matrix A(fem::petsc::create_matrix(*a), true);
+  la::petsc::Matrix A_coarse(fem::petsc::create_matrix(*a_coarse), true);
+
+  la::Vector<T> b(L->function_spaces()[0]->dofmap()->index_map,
+                  L->function_spaces()[0]->dofmap()->index_map_bs());
+
+  // TOOD: this somehow breaks?!?
+  // V->mesh()->topology_mutable()->create_connectivity(2, 3);
+  // mesh::exterior_facet_indices(*V->mesh()->topology())
+  auto bc = std::make_shared<const fem::DirichletBC<T>>(
+      0.0,
+      mesh::locate_entities_boundary(
+          *V->mesh(), 0,
+          [](auto x) { return std::vector<std::int8_t>(x.extent(1), true); }),
+      V);
+
+  // V_coarse->mesh()->topology_mutable()->create_connectivity(2, 3);
+  auto bc_coarse = std::make_shared<const fem::DirichletBC<T>>(
+      0.0,
+      mesh::locate_entities_boundary(
+          *V_coarse->mesh(), 0,
+          [](auto x) { return std::vector<std::int8_t>(x.extent(1), true); }),
+      V_coarse);
+
+  // assemble A
+  MatZeroEntries(A.mat());
+  fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A.mat(), ADD_VALUES), *a,
+                       {bc});
+  MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY);
+  MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY);
+  fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A.mat(), INSERT_VALUES), *V,
+                       {bc});
+  MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);
+
+  // assemble b
+  b.set(0.0);
+  fem::assemble_vector(b.mutable_array(), *L);
+  fem::apply_lifting<T, U>(b.mutable_array(), {a}, {{bc}}, {}, T(1));
+  b.scatter_rev(std::plus<T>());
+  fem::set_bc<T, U>(b.mutable_array(), {bc});
+
+  // assemble A_coarse
+  MatZeroEntries(A_coarse.mat());
+  fem::assemble_matrix(
+      la::petsc::Matrix::set_block_fn(A_coarse.mat(), ADD_VALUES), *a_coarse,
+      {bc_coarse});
+  MatAssemblyBegin(A_coarse.mat(), MAT_FLUSH_ASSEMBLY);
+  MatAssemblyEnd(A_coarse.mat(), MAT_FLUSH_ASSEMBLY);
+  fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A_coarse.mat(), INSERT_VALUES),
+                       *V_coarse, {bc_coarse});
+  MatAssemblyBegin(A_coarse.mat(), MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(A_coarse.mat(), MAT_FINAL_ASSEMBLY);
+
+  KSP ksp;
+  KSPCreate(MPI_COMM_WORLD, &ksp);
+  KSPSetType(ksp, "cg");
+
+  PC pc;
+  KSPGetPC(ksp, &pc);
+  KSPSetFromOptions(ksp);
+  PCSetType(pc, "mg");
+
+  PCMGSetLevels(pc, 2, NULL);
+  PCMGSetType(pc, PC_MG_MULTIPLICATIVE);
+  PCMGSetCycleType(pc, PC_MG_CYCLE_V);
+
+  // do not rely on coarse grid operators to be generated by
+  // restriction/prolongation
+  PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
+  PCMGSetOperators(pc, 0, A_coarse.mat(), A_coarse.mat());
+
+  // PCMGSetNumberSmooth(pc, 1);
+  PCSetFromOptions(pc);
+
+  mesh.topology_mutable()->create_connectivity(0, 1);
+  mesh.topology_mutable()->create_connectivity(1, 0);
+  auto inclusion_map = transfer::inclusion_mapping(*mesh_coarse, mesh);
+  la::SparsityPattern sp
+      = transfer::create_sparsity_pattern(*V_coarse, *V, inclusion_map);
+  la::petsc::Matrix restriction(MPI_COMM_WORLD, sp);
+  transfer::assemble_transfer_matrix<double>(
+      la::petsc::Matrix::set_block_fn(restriction.mat(), INSERT_VALUES),
+      *V_coarse, *V, inclusion_map,
+      [](std::int32_t d) -> double { return d == 0 ? 1. : .5; });
+
+  MatAssemblyBegin(restriction.mat(), MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(restriction.mat(), MAT_FINAL_ASSEMBLY);
+
+  // PETSc figures out to use transpose by dimensions!
+  PCMGSetInterpolation(pc, 1, restriction.mat());
+  PCMGSetRestriction(pc, 1, restriction.mat());
+
+  // MatView(A.mat(), PETSC_VIEWER_STDOUT_SELF);
+  KSPSetOperators(ksp, A.mat(), A.mat());
+  KSPSetUp(ksp);
+
+  auto u = std::make_shared<fem::Function<T>>(V);
+
+  la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
+  la::petsc::Vector _b(la::petsc::create_vector_wrap(b), false);
+
+  KSPSolve(ksp, _b.vec(), _u.vec());
+  u->x()->scatter_fwd();
+
+  {
+    io::VTKFile file(MPI_COMM_WORLD, "u.pvd", "w");
+    file.write<T>({*u}, 0.0);
+
+    io::XDMFFile file_xdmf(mesh.comm(), "u_xdmf.xdmf", "w");
+    file_xdmf.write_mesh(mesh);
+    file_xdmf.write_function(*u, 0.0);
+    file_xdmf.close();
+
+#ifdef HAS_ADIOS2
+    io::VTXWriter<U> vtx_writer(mesh.comm(), std::filesystem::path("u_vtx.bp"),
+                                {u}, io::VTXMeshPolicy::reuse);
+    vtx_writer.write(0);
+#endif
+  }
+
+  KSPDestroy(&ksp);
+
+  // PetscLogView(PETSC_VIEWER_STDOUT_SELF);
+}
+
+/// @endcond
\ No newline at end of file
diff --git a/cpp/demo/geometric-twogrid/poisson.py b/cpp/demo/geometric-twogrid/poisson.py
new file mode 100644
index 00000000000..191ce33347e
--- /dev/null
+++ b/cpp/demo/geometric-twogrid/poisson.py
@@ -0,0 +1,44 @@
+# The first step is to define the variational problem at hand. We define
+# the variational problem in UFL terms in a separate form file
+# {download}`demo_poisson/poisson.py`.  We begin by defining the finite
+# element:
+
+from basix.ufl import element
+from ufl import (
+    Coefficient,
+    FunctionSpace,
+    Mesh,
+    TestFunction,
+    TrialFunction,
+    dx,
+    grad,
+    inner,
+)
+
+e = element("Lagrange", "tetrahedron", 1)
+
+# The first argument to :py:class:`FiniteElement` is the finite element
+# family, the second argument specifies the domain, while the third
+# argument specifies the polynomial degree. Thus, in this case, our
+# element `element` consists of first-order, continuous Lagrange basis
+# functions on triangles (or in order words, continuous piecewise linear
+# polynomials on triangles).
+#
+# Next, we use this element to initialize the trial and test functions
+# ($u$ and $v$) and the coefficient functions ($f$ and $g$):
+
+coord_element = element("Lagrange", "tetrahedron", 1, shape=(3,))
+mesh = Mesh(coord_element)
+
+V = FunctionSpace(mesh, e)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Coefficient(V)
+g = Coefficient(V)
+
+# Finally, we define the bilinear and linear forms according to the
+# variational formulation of the equations:
+
+a = inner(grad(u), grad(v)) * dx
+L = inner(f, v) * dx