From a51a9fb32a25e2c74f4582245db9f337d1c9e432 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Fri, 22 Mar 2024 14:34:14 +0000
Subject: [PATCH 01/29] First attempt on submesh interpolation

---
 cpp/dolfinx/fem/Function.h      |  51 ++++++--
 cpp/dolfinx/fem/interpolate.h   | 212 ++++++++++++++++++--------------
 python/dolfinx/fem/function.py  |   4 +
 python/dolfinx/wrappers/fem.cpp |  13 +-
 4 files changed, 178 insertions(+), 102 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 411a4237e0..4d56f23827 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -153,28 +153,32 @@ class Function
   /// @brief Interpolate a provided Function.
   /// @param[in] v The function to be interpolated
   /// @param[in] cells The cells to interpolate on
+  /// @param[in] cell_map Map from indices in `cells` to cell indices in `v`
   /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
   /// nonmatching meshes. This data can be generated with
   /// generate_nonmatching_meshes_interpolation_data (optional).
   void interpolate(
       const Function<value_type, geometry_type>& v,
       std::span<const std::int32_t> cells,
+      std::span<const std::int32_t> cell_map,
       const std::tuple<std::span<const std::int32_t>,
                        std::span<const std::int32_t>,
                        std::span<const geometry_type>,
                        std::span<const std::int32_t>>& nmm_interpolation_data
       = {})
   {
-    fem::interpolate(*this, v, cells, nmm_interpolation_data);
+    fem::interpolate(*this, v, cells, cell_map, nmm_interpolation_data);
   }
 
   /// @brief Interpolate a provided Function.
   /// @param[in] v The function to be interpolated
+  /// @param[in] cell_map Map from cells in self to cell indices in `v`
   /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
   /// nonmatching meshes. This data can be generated with
   /// generate_nonmatching_meshes_interpolation_data (optional).
   void interpolate(
       const Function<value_type, geometry_type>& v,
+      std::span<const std::int32_t> cell_map,
       const std::tuple<std::span<const std::int32_t>,
                        std::span<const std::int32_t>,
                        std::span<const geometry_type>,
@@ -184,12 +188,12 @@ class Function
     assert(_function_space);
     assert(_function_space->mesh());
     int tdim = _function_space->mesh()->topology()->dim();
-    auto cell_map = _function_space->mesh()->topology()->index_map(tdim);
-    assert(cell_map);
-    std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
+    auto cell_imap = _function_space->mesh()->topology()->index_map(tdim);
+    assert(cell_imap);
+    std::int32_t num_cells = cell_imap->size_local() + cell_imap->num_ghosts();
     std::vector<std::int32_t> cells(num_cells, 0);
     std::iota(cells.begin(), cells.end(), 0);
-    interpolate(v, cells, nmm_interpolation_data);
+    interpolate(v, cells, cell_map, nmm_interpolation_data);
   }
 
   /// Interpolate an expression function on a list of cells
@@ -283,8 +287,12 @@ class Function
   /// `FiniteElement::interpolation_points()` for the element associated
   /// with `u`.
   /// @param[in] cells The cells to interpolate on
+  /// @param[in] expr_mesh Mesh expression is defined on
+  /// @param[in] cell_map Map from `cells` to cells in expression
   void interpolate(const Expression<value_type, geometry_type>& e,
-                   std::span<const std::int32_t> cells)
+                   std::span<const std::int32_t> cells,
+                   const dolfinx::mesh::Mesh& expr_mesh,
+                   std::span<const std::int32_t> cell_map)
   {
     // Check that spaces are compatible
     assert(_function_space);
@@ -331,7 +339,28 @@ class Function
 
     // Evaluate Expression at points
     assert(_function_space->mesh());
-    e.eval(*_function_space->mesh(), cells, fdata,
+
+    //
+    std::vector<std::int32_t> cells_expr;
+    cells_expr.reserve(num_cells);
+    // Get mesh and check if mesh is shared
+    if (auto mesh_v = _function_space()->mesh(); expr_mesh == mesh_v)
+    {
+      cells_epxr.insert(cells_epxr.end(), cells_u.begin(), cells_u.end());
+    }
+    // If meshes are different and input mapping is given
+    else if (cell_map.size() > 0)
+    {
+      std::transform(cells_u.begin(), cells_u.end(),
+                     std::back_inserter(cells_epxr),
+                     [&cell_map](std::int32_t c) { return cell_map[c]; });
+    }
+    else
+    {
+      std::runtime_error("Meshes are different and no cell map is provided")
+    }
+
+    e.eval(*_function_space->mesh(), cells_expr, fdata,
            {num_cells, num_points * value_size});
 
     // Reshape evaluated data to fit interpolate
@@ -356,7 +385,11 @@ class Function
 
   /// Interpolate an Expression (based on UFL) on all cells
   /// @param[in] e The function to be interpolated
-  void interpolate(const Expression<value_type, geometry_type>& e)
+  /// @param[in] expr_mesh Mesh expression is defined on
+  /// @param[in] cell_map Map from `cells` to cells in expression
+  void interpolate(const Expression<value_type, geometry_type>& e,
+                   const dolfinx::mesh::Mesh& expr_mesh,
+                   std::span<const std::int32_t> cell_map)
   {
     assert(_function_space);
     assert(_function_space->mesh());
@@ -366,7 +399,7 @@ class Function
     std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
     std::vector<std::int32_t> cells(num_cells, 0);
     std::iota(cells.begin(), cells.end(), 0);
-    interpolate(e, cells);
+    interpolate(e, cells, expr_mesh, cell_map);
   }
 
   /// @brief Evaluate the Function at points.
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 9507382876..0da8877c19 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -321,38 +321,46 @@ void interpolation_apply(U&& Pi, V&& data, std::span<T> coeffs, int bs)
 /// map.
 /// @param[out] u1 The function to interpolate to
 /// @param[in] u0 The function to interpolate from
-/// @param[in] cells The cells to interpolate on
+/// @param[in] cells_u1 The cells to interpolate on
+/// @param[in] cells_u0 Equivalent cell in u0 for each cell in u1
 /// @pre The functions `u1` and `u0` must share the same mesh and the
 /// elements must share the same basis function map. Neither is checked
 /// by the function.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
-                          std::span<const std::int32_t> cells)
+                          std::span<const std::int32_t> cells_u1,
+                          std::span<const std::int32_t> cells_u0)
 {
   auto V0 = u0.function_space();
   assert(V0);
   auto V1 = u1.function_space();
   assert(V1);
-  auto mesh = V0->mesh();
-  assert(mesh);
+  auto mesh0 = V0->mesh();
+  assert(mesh0);
+
+  auto mesh1 = V1->mesh();
+  assert(mesh1);
 
   auto element0 = V0->element();
   assert(element0);
   auto element1 = V1->element();
   assert(element1);
 
-  const int tdim = mesh->topology()->dim();
-  auto map = mesh->topology()->index_map(tdim);
+  const int tdim = mesh0->topology()->dim();
+  auto map = mesh0->topology()->index_map(tdim);
   assert(map);
   std::span<T> u1_array = u1.x()->mutable_array();
   std::span<const T> u0_array = u0.x()->array();
 
-  std::span<const std::uint32_t> cell_info;
+  std::span<const std::uint32_t> cell_info0;
+  std::span<const std::uint32_t> cell_info1;
   if (element1->needs_dof_transformations()
       or element0->needs_dof_transformations())
   {
-    mesh->topology_mutable()->create_entity_permutations();
-    cell_info = std::span(mesh->topology()->get_cell_permutation_info());
+    mesh0->topology_mutable()->create_entity_permutations();
+    cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
+    mesh1->topology_mutable()->create_entity_permutations();
+    cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
   }
 
   // Get dofmaps
@@ -379,14 +387,14 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
 
   // Iterate over mesh and interpolate on each cell
   using X = typename dolfinx::scalar_value_type_t<T>;
-  for (auto c : cells)
+  for (std::size_t c = 0; c < cells_u0.size(); c++)
   {
-    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
+    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_u0[c]);
     for (std::size_t i = 0; i < dofs0.size(); ++i)
       for (int k = 0; k < bs0; ++k)
         local0[bs0 * i + k] = u0_array[bs0 * dofs0[i] + k];
 
-    apply_dof_transformation(local0, cell_info, c, 1);
+    apply_dof_transformation(local0, cell_info0, cells_u0[c], 1);
 
     // FIXME: Get compile-time ranges from Basix
     // Apply interpolation operator
@@ -395,8 +403,8 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
       for (std::size_t j = 0; j < im_shape[1]; ++j)
         local1[i] += static_cast<X>(i_m[im_shape[1] * i + j]) * local0[j];
 
-    apply_inverse_dof_transform(local1, cell_info, c, 1);
-    std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
+    apply_inverse_dof_transform(local1, cell_info1, cells_u1[c], 1);
+    std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells_u1[c]);
     for (std::size_t i = 0; i < dofs1.size(); ++i)
       for (int k = 0; k < bs1; ++k)
         u1_array[bs1 * dofs1[i] + k] = local1[bs1 * i + k];
@@ -410,37 +418,45 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
 /// standard isoparametric mapping.
 /// @param[out] u1 The function to interpolate to
 /// @param[in] u0 The function to interpolate from
-/// @param[in] cells The cells to interpolate on
+/// @param[in] cells_u1 The cells to interpolate on
+/// @param[in] cells_u0 Equivalent cell in u0 for each cell in u1
 /// @pre The functions `u1` and `u0` must share the same mesh. This is
 /// not checked by the function.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
-                                  std::span<const std::int32_t> cells)
+                                  std::span<const std::int32_t> cells_u1,
+                                  std::span<const std::int32_t> cells_u0)
 {
   // Get mesh
   auto V0 = u0.function_space();
   assert(V0);
-  auto mesh = V0->mesh();
-  assert(mesh);
+  auto mesh0 = V0->mesh();
+  assert(mesh0);
 
   // Mesh dims
-  const int tdim = mesh->topology()->dim();
-  const int gdim = mesh->geometry().dim();
+  const int tdim = mesh0->topology()->dim();
+  const int gdim = mesh0->geometry().dim();
 
   // Get elements
   auto V1 = u1.function_space();
   assert(V1);
+  auto mesh1 = V1->mesh();
+  assert(mesh1);
   auto element0 = V0->element();
   assert(element0);
   auto element1 = V1->element();
   assert(element1);
 
-  std::span<const std::uint32_t> cell_info;
+  std::span<const std::uint32_t> cell_info0;
+  std::span<const std::uint32_t> cell_info1;
+
   if (element1->needs_dof_transformations()
       or element0->needs_dof_transformations())
   {
-    mesh->topology_mutable()->create_entity_permutations();
-    cell_info = std::span(mesh->topology()->get_cell_permutation_info());
+    mesh0->topology_mutable()->create_entity_permutations();
+    cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
+    mesh1->topology_mutable()->create_entity_permutations();
+    cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
   }
 
   // Get dofmaps
@@ -464,10 +480,10 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
   const std::size_t value_size_ref0 = element0->reference_value_size() / bs0;
   const std::size_t value_size0 = V0->value_size() / bs0;
 
-  const CoordinateElement<U>& cmap = mesh->geometry().cmap();
-  auto x_dofmap = mesh->geometry().dofmap();
+  const CoordinateElement<U>& cmap = mesh0->geometry().cmap();
+  auto x_dofmap = mesh0->geometry().dofmap();
   const std::size_t num_dofs_g = cmap.dim();
-  std::span<const U> x_g = mesh->geometry().x();
+  std::span<const U> x_g = mesh0->geometry().x();
 
   // Evaluate coordinate map basis at reference interpolation points
   const std::array<std::size_t, 4> phi_shape
@@ -536,11 +552,11 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
   // Iterate over mesh and interpolate on each cell
   std::span<const T> array0 = u0.x()->array();
   std::span<T> array1 = u1.x()->mutable_array();
-  for (auto c : cells)
+  for (std::size_t c = 0; c < cells_u0.size(); c++)
   {
     // Get cell geometry (coordinate dofs)
     auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
-        x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
+        x_dofmap, cells_u0[c], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
     for (std::size_t i = 0; i < num_dofs_g; ++i)
     {
       const int pos = 3 * x_dofs[i];
@@ -580,7 +596,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
       apply_dof_transformation0(
           std::span(basis_reference0_b.data() + p * dim0 * value_size_ref0,
                     dim0 * value_size_ref0),
-          cell_info, c, value_size_ref0);
+          cell_info0, cells_u0[c], value_size_ref0);
     }
 
     for (std::size_t i = 0; i < basis0.extent(0); ++i)
@@ -602,7 +618,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
 
     // Copy expansion coefficients for v into local array
     const int dof_bs0 = dofmap0->bs();
-    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
+    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_u0[c]);
     for (std::size_t i = 0; i < dofs0.size(); ++i)
       for (int k = 0; k < dof_bs0; ++k)
         coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k];
@@ -645,7 +661,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
         mapped_values0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0,
         MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
     interpolation_apply(Pi_1, values, std::span(local1), bs1);
-    apply_inverse_dof_transform1(local1, cell_info, c, 1);
+    apply_inverse_dof_transform1(local1, cell_info1, cells_u1[c], 1);
 
     // Copy local coefficients to the correct position in u dof array
     const int dof_bs1 = dofmap1->bs();
@@ -1097,14 +1113,16 @@ create_nonmatching_meshes_interpolation_data(const mesh::Mesh<T>& mesh0,
 /// @brief Interpolate from one finite element Function to another one.
 /// @param[out] u The function to interpolate into
 /// @param[in] v The function to be interpolated
-/// @param[in] cells List of cell indices to interpolate on
+/// @param[in] cells_u List of cell indices to interpolate on
+/// @param[in] cell_map Mapping from cells in `u` to cells in `v`
 /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
 /// nonmatching meshes. This data can be generated with
 /// create_nonmatching_meshes_interpolation_data (optional).
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate(
     Function<T, U>& u, const Function<T, U>& v,
-    std::span<const std::int32_t> cells,
+    std::span<const std::int32_t> cells_u,
+    std::span<const std::int32_t> cell_map,
     const std::tuple<std::span<const std::int32_t>,
                      std::span<const std::int32_t>, std::span<const U>,
                      std::span<const std::int32_t>>& nmm_interpolation_data
@@ -1118,7 +1136,7 @@ void interpolate(
   auto cell_map0 = mesh->topology()->index_map(mesh->topology()->dim());
   assert(cell_map0);
   std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
-  if (u.function_space() == v.function_space() and cells.size() == num_cells0)
+  if (u.function_space() == v.function_space() and cells_u.size() == num_cells0)
   {
     // Same function spaces and on whole mesh
     std::span<T> u1_array = u.x()->mutable_array();
@@ -1127,79 +1145,91 @@ void interpolate(
   }
   else
   {
+    std::vector<std::int32_t> cells_v;
+    cells_v.reserve(cells_u.size());
     // Get mesh and check that functions share the same mesh
-    if (auto mesh_v = v.function_space()->mesh(); mesh != mesh_v)
+    if (auto mesh_v = v.function_space()->mesh(); mesh == mesh_v)
     {
-      impl::interpolate_nonmatching_meshes(u, v, cells, nmm_interpolation_data);
-      return;
+      cells_v.insert(cells_v.end(), cells_u.begin(), cells_u.end());
     }
-    else
+    // If meshes are different and input mapping is given
+    else if (cell_map.size() > 0)
+    {
+      std::transform(cells_u.begin(), cells_u.end(),
+                     std::back_inserter(cells_v),
+                     [&cell_map](std::int32_t c) { return cell_map[c]; });
+    }
+    // Non-matching meshes
+    if (cells_v.size() == 0)
     {
+      impl::interpolate_nonmatching_meshes(u, v, cells_u,
+                                           nmm_interpolation_data);
+      return;
+    }
 
-      // Get elements and check value shape
-      auto fs0 = v.function_space();
-      auto element0 = fs0->element();
-      assert(element0);
-      auto fs1 = u.function_space();
-      auto element1 = fs1->element();
-      assert(element1);
-      if (fs0->value_shape().size() != fs1->value_shape().size()
-          or !std::equal(fs0->value_shape().begin(), fs0->value_shape().end(),
-                         fs1->value_shape().begin()))
-      {
-        throw std::runtime_error(
-            "Interpolation: elements have different value dimensions");
-      }
+    // Get elements and check value shape
+    auto fs0 = v.function_space();
+    auto element0 = fs0->element();
+    assert(element0);
+    auto fs1 = u.function_space();
+    auto element1 = fs1->element();
+    assert(element1);
+    if (fs0->value_shape().size() != fs1->value_shape().size()
+        or !std::equal(fs0->value_shape().begin(), fs0->value_shape().end(),
+                       fs1->value_shape().begin()))
+    {
+      throw std::runtime_error(
+          "Interpolation: elements have different value dimensions");
+    }
 
-      if (element1 == element0 or *element1 == *element0)
-      {
-        // Same element, different dofmaps (or just a subset of cells)
+    if (element1 == element0 or *element1 == *element0)
+    {
+      // Same element, different dofmaps (or just a subset of cells)
 
-        const int tdim = mesh->topology()->dim();
-        auto cell_map = mesh->topology()->index_map(tdim);
-        assert(cell_map);
+      const int tdim = mesh->topology()->dim();
+      auto cell_map1 = mesh->topology()->index_map(tdim);
+      assert(cell_map1);
 
-        assert(element1->block_size() == element0->block_size());
+      assert(element1->block_size() == element0->block_size());
 
-        // Get dofmaps
-        std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
-        assert(dofmap0);
-        std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
-        assert(dofmap1);
+      // Get dofmaps
+      std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
+      assert(dofmap0);
+      std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
+      assert(dofmap1);
 
-        std::span<T> u1_array = u.x()->mutable_array();
-        std::span<const T> u0_array = v.x()->array();
+      std::span<T> u1_array = u.x()->mutable_array();
+      std::span<const T> u0_array = v.x()->array();
 
-        // Iterate over mesh and interpolate on each cell
-        const int bs0 = dofmap0->bs();
-        const int bs1 = dofmap1->bs();
-        for (auto c : cells)
+      // Iterate over mesh and interpolate on each cell
+      const int bs0 = dofmap0->bs();
+      const int bs1 = dofmap1->bs();
+      for (std::size_t c = 0; c < cells_u.size(); ++c)
+      {
+        std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_v[c]);
+        std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells_u[c]);
+        assert(bs0 * dofs0.size() == bs1 * dofs1.size());
+        for (std::size_t i = 0; i < dofs0.size(); ++i)
         {
-          std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
-          std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
-          assert(bs0 * dofs0.size() == bs1 * dofs1.size());
-          for (std::size_t i = 0; i < dofs0.size(); ++i)
+          for (int k = 0; k < bs0; ++k)
           {
-            for (int k = 0; k < bs0; ++k)
-            {
-              int index = bs0 * i + k;
-              std::div_t dv1 = std::div(index, bs1);
-              u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
-                  = u0_array[bs0 * dofs0[i] + k];
-            }
+            int index = bs0 * i + k;
+            std::div_t dv1 = std::div(index, bs1);
+            u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
+                = u0_array[bs0 * dofs0[i] + k];
           }
         }
       }
-      else if (element1->map_type() == element0->map_type())
-      {
-        // Different elements, same basis function map type
-        impl::interpolate_same_map(u, v, cells);
-      }
-      else
-      {
-        //  Different elements with different maps for basis functions
-        impl::interpolate_nonmatching_maps(u, v, cells);
-      }
+    }
+    else if (element1->map_type() == element0->map_type())
+    {
+      // Different elements, same basis function map type
+      impl::interpolate_same_map(u, v, cells_u, cells_v);
+    }
+    else
+    {
+      //  Different elements with different maps for basis functions
+      impl::interpolate_nonmatching_maps(u, v, cells_u, cells_v);
     }
   }
 }
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index da209ab070..ae2676a461 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -420,6 +420,7 @@ def interpolate(
         self,
         u: typing.Union[typing.Callable, Expression, Function],
         cells: typing.Optional[np.ndarray] = None,
+        cell_map: typing.Optional[np.ndarray] = None,
         nmm_interpolation_data: typing.Optional[PointOwnershipData] = None,
     ) -> None:
         """Interpolate an expression
@@ -428,6 +429,7 @@ def interpolate(
             u: The function, Expression or Function to interpolate.
             cells: The cells to interpolate over. If `None` then all
                 cells are interpolated over.
+            cell_map: Mapping from cells in receiving mesh to cells in `u` mesh
             nmm_interpolation_data: Data needed to interpolate functions defined on other meshes
         """
         if nmm_interpolation_data is None:
@@ -438,6 +440,8 @@ def interpolate(
                 dest_points=np.empty(0, dtype=x_dtype),
                 dest_cells=np.empty(0, dtype=np.int32),
             )
+        if cell_map is None:
+            cell_map = np.empty(0, dtype=np.int32)
 
         @singledispatch
         def _interpolate(u, cells: typing.Optional[np.ndarray] = None):
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index aaa172ada0..0a2043e5dd 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -391,6 +391,8 @@ void declare_objects(nb::module_& m, const std::string& type)
           [](dolfinx::fem::Function<T, U>& self,
              dolfinx::fem::Function<T, U>& u,
              nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> cells,
+             nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>
+                 cell_map,
              const std::tuple<
                  nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>,
                  nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>,
@@ -414,6 +416,7 @@ void declare_objects(nb::module_& m, const std::string& type)
                         std::get<3>(interpolation_data).data(),
                         std::get<3>(interpolation_data).size()));
             self.interpolate(u, std::span(cells.data(), cells.size()),
+                             std::span(cell_map.data(), cell_map.size()),
                              _interpolation_data);
           },
           nb::arg("u"), nb::arg("cells"), nb::arg("nmm_interpolation_data"),
@@ -453,8 +456,14 @@ void declare_objects(nb::module_& m, const std::string& type)
           "interpolate",
           [](dolfinx::fem::Function<T, U>& self,
              const dolfinx::fem::Expression<T, U>& expr,
-             nb::ndarray<const std::int32_t, nb::c_contig> cells)
-          { self.interpolate(expr, std::span(cells.data(), cells.size())); },
+             nb::ndarray<const std::int32_t, nb::c_contig> cells,
+             const dolfinx::mesh::Mesh<U>& expr_mesh,
+             nb::ndarray<const std::int32_t, nb::c_contig> cell_map)
+          {
+            self.interpolate(expr, std::span(cells.data(), cells.size()),
+                             expr_mesh,
+                             std::span(cell_map.data(), cell_map.size()));
+          },
           nb::arg("expr"), nb::arg("cells"),
           "Interpolate an Expression on a set of cells")
       .def_prop_ro(

From 2368272c6ef9b49738c5681d0182107de6f9e51a Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Sun, 24 Mar 2024 17:58:36 +0000
Subject: [PATCH 02/29] ADd some more conversion

---
 cpp/dolfinx/fem/Function.h      | 23 ++++++++++++-----------
 python/dolfinx/fem/function.py  | 33 +++++++++++++++++++++++----------
 python/dolfinx/wrappers/fem.cpp |  7 ++++---
 3 files changed, 39 insertions(+), 24 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 4d56f23827..158978b35f 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -291,7 +291,7 @@ class Function
   /// @param[in] cell_map Map from `cells` to cells in expression
   void interpolate(const Expression<value_type, geometry_type>& e,
                    std::span<const std::int32_t> cells,
-                   const dolfinx::mesh::Mesh& expr_mesh,
+                   const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
                    std::span<const std::int32_t> cell_map)
   {
     // Check that spaces are compatible
@@ -344,20 +344,20 @@ class Function
     std::vector<std::int32_t> cells_expr;
     cells_expr.reserve(num_cells);
     // Get mesh and check if mesh is shared
-    if (auto mesh_v = _function_space()->mesh(); expr_mesh == mesh_v)
+    if (auto mesh_v = _function_space->mesh();
+        expr_mesh.topology() == mesh_v->topology())
     {
-      cells_epxr.insert(cells_epxr.end(), cells_u.begin(), cells_u.end());
+      cells_expr.insert(cells_expr.end(), cells.begin(), cells.end());
     }
     // If meshes are different and input mapping is given
     else if (cell_map.size() > 0)
     {
-      std::transform(cells_u.begin(), cells_u.end(),
-                     std::back_inserter(cells_epxr),
+      std::transform(cells.begin(), cells.end(), std::back_inserter(cells_expr),
                      [&cell_map](std::int32_t c) { return cell_map[c]; });
     }
     else
     {
-      std::runtime_error("Meshes are different and no cell map is provided")
+      std::runtime_error("Meshes are different and no cell map is provided");
     }
 
     e.eval(*_function_space->mesh(), cells_expr, fdata,
@@ -388,18 +388,19 @@ class Function
   /// @param[in] expr_mesh Mesh expression is defined on
   /// @param[in] cell_map Map from `cells` to cells in expression
   void interpolate(const Expression<value_type, geometry_type>& e,
-                   const dolfinx::mesh::Mesh& expr_mesh,
+                   const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
                    std::span<const std::int32_t> cell_map)
   {
     assert(_function_space);
     assert(_function_space->mesh());
     const int tdim = _function_space->mesh()->topology()->dim();
-    auto cell_map = _function_space->mesh()->topology()->index_map(tdim);
-    assert(cell_map);
-    std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
+    auto cell_imap = _function_space->mesh()->topology()->index_map(tdim);
+    assert(cell_imap);
+    std::int32_t num_cells = cell_imap->size_local() + cell_imap->num_ghosts();
     std::vector<std::int32_t> cells(num_cells, 0);
     std::iota(cells.begin(), cells.end(), 0);
-    interpolate(e, cells, expr_mesh, cell_map);
+    interpolate(e, cells, expr_mesh,
+                std::span(cell_map.data(), cell_map.size()));
   }
 
   /// @brief Evaluate the Function at points.
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index ae2676a461..ea5fbf9b08 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -421,6 +421,7 @@ def interpolate(
         u: typing.Union[typing.Callable, Expression, Function],
         cells: typing.Optional[np.ndarray] = None,
         cell_map: typing.Optional[np.ndarray] = None,
+        expr_mesh: typing.Optional[Mesh] = None,
         nmm_interpolation_data: typing.Optional[PointOwnershipData] = None,
     ) -> None:
         """Interpolate an expression
@@ -440,8 +441,18 @@ def interpolate(
                 dest_points=np.empty(0, dtype=x_dtype),
                 dest_cells=np.empty(0, dtype=np.int32),
             )
+
+        if cells is None:
+            mesh = self.function_space.mesh
+            map = mesh.topology.index_map(mesh.topology.dim)
+            cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
+
         if cell_map is None:
-            cell_map = np.empty(0, dtype=np.int32)
+            # If cell map is not provided create identity map
+            cell_map = np.arange(len(cells), dtype=np.int32)
+            assert expr_mesh is None
+            expr_mesh = self.function_space.mesh
+
 
         @singledispatch
         def _interpolate(u, cells: typing.Optional[np.ndarray] = None):
@@ -451,7 +462,7 @@ def _interpolate(u, cells: typing.Optional[np.ndarray] = None):
         @_interpolate.register(Function)
         def _(u: Function, cells: typing.Optional[np.ndarray] = None):
             """Interpolate a fem.Function"""
-            self._cpp_object.interpolate(u._cpp_object, cells, nmm_interpolation_data)  # type: ignore
+            self._cpp_object.interpolate(u._cpp_object, cells, cell_map, nmm_interpolation_data)  # type: ignore
 
         @_interpolate.register(int)
         def _(u_ptr: int, cells: typing.Optional[np.ndarray] = None):
@@ -459,14 +470,16 @@ def _(u_ptr: int, cells: typing.Optional[np.ndarray] = None):
             self._cpp_object.interpolate_ptr(u_ptr, cells)  # type: ignore
 
         @_interpolate.register(Expression)
-        def _(expr: Expression, cells: typing.Optional[np.ndarray] = None):
-            """Interpolate Expression for the set of cells"""
-            self._cpp_object.interpolate(expr._cpp_object, cells)  # type: ignore
-
-        if cells is None:
-            mesh = self.function_space.mesh
-            map = mesh.topology.index_map(mesh.topology.dim)
-            cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
+        def _(expr: Expression,
+              cells: typing.Optional[np.ndarray] = None
+              ):
+            """Interpolate Expression from a given mesh onto the set of cells
+            Args:
+                expr: Expression to interpolate
+                cells: The cells to interpolate over. If `None` then all
+                    cells are interpolated over.
+            """
+            self._cpp_object.interpolate(expr._cpp_object, cells, expr_mesh, cell_map)  # type: ignore
 
         try:
             # u is a Function or Expression (or pointer to one)
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 0a2043e5dd..03bcf64fc1 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -419,7 +419,8 @@ void declare_objects(nb::module_& m, const std::string& type)
                              std::span(cell_map.data(), cell_map.size()),
                              _interpolation_data);
           },
-          nb::arg("u"), nb::arg("cells"), nb::arg("nmm_interpolation_data"),
+          nb::arg("u"), nb::arg("cells"), nb::arg("cell_map"),
+          nb::arg("nmm_interpolation_data"),
           "Interpolate a finite element function")
       .def(
           "interpolate_ptr",
@@ -464,8 +465,8 @@ void declare_objects(nb::module_& m, const std::string& type)
                              expr_mesh,
                              std::span(cell_map.data(), cell_map.size()));
           },
-          nb::arg("expr"), nb::arg("cells"),
-          "Interpolate an Expression on a set of cells")
+          nb::arg("expr"), nb::arg("cells"), nb::arg("expr_mesh"),
+          nb::arg("cell_map"), "Interpolate an Expression on a set of cells")
       .def_prop_ro(
           "x", nb::overload_cast<>(&dolfinx::fem::Function<T, U>::x),
           "Return the vector associated with the finite element Function")

From dca89f446b6e224d417cdc04bcf48a3a3a56213a Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 19:04:54 +0000
Subject: [PATCH 03/29] UPdate more interfaces

---
 cpp/dolfinx/fem/Function.h     | 10 +++++++---
 python/dolfinx/fem/function.py | 20 ++++++++++++++------
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 158978b35f..d762b9bb08 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -386,7 +386,8 @@ class Function
   /// Interpolate an Expression (based on UFL) on all cells
   /// @param[in] e The function to be interpolated
   /// @param[in] expr_mesh Mesh expression is defined on
-  /// @param[in] cell_map Map from `cells` to cells in expression
+  /// @param[in] cell_map Map from `cells` to cells in expression if
+  /// receiving function is defined on a different mesh than the expression
   void interpolate(const Expression<value_type, geometry_type>& e,
                    const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
                    std::span<const std::int32_t> cell_map)
@@ -399,8 +400,11 @@ class Function
     std::int32_t num_cells = cell_imap->size_local() + cell_imap->num_ghosts();
     std::vector<std::int32_t> cells(num_cells, 0);
     std::iota(cells.begin(), cells.end(), 0);
-    interpolate(e, cells, expr_mesh,
-                std::span(cell_map.data(), cell_map.size()));
+
+    if (cell_map.size() == 0)
+      interpolate(e, cells, expr_mesh, cell_map);
+    else
+      interpolate(e, cells, expr_mesh, std::span(cells.data(), cells.size()));
   }
 
   /// @brief Evaluate the Function at points.
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index ea5fbf9b08..62202d2823 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -448,11 +448,7 @@ def interpolate(
             cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
 
         if cell_map is None:
-            # If cell map is not provided create identity map
-            cell_map = np.arange(len(cells), dtype=np.int32)
-            assert expr_mesh is None
-            expr_mesh = self.function_space.mesh
-
+            cell_map = np.empty(0, dtype=np.int32)
 
         @singledispatch
         def _interpolate(u, cells: typing.Optional[np.ndarray] = None):
@@ -479,7 +475,19 @@ def _(expr: Expression,
                 cells: The cells to interpolate over. If `None` then all
                     cells are interpolated over.
             """
-            self._cpp_object.interpolate(expr._cpp_object, cells, expr_mesh, cell_map)  # type: ignore
+            if len(cell_map) == 0:
+                # If cell map is not provided create identity map
+                expr_cell_map = np.arange(len(cells), dtype=np.int32)
+                assert expr_mesh is None
+                mapping_mesh = self.function_space.mesh._cpp_object
+            else:
+                # If cell map is provided check that there is a mesh
+                # associated with the expression
+                expr_cell_map = cell_map
+                assert expr_mesh is not None
+                mapping_mesh = expr_mesh._cpp_object
+            self._cpp_object.interpolate(expr._cpp_object, cells,
+                                         mapping_mesh, expr_cell_map)  # type: ignore
 
         try:
             # u is a Function or Expression (or pointer to one)

From 4e59abb6dfccbc84543e8f8af48a5536e4e5c764 Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 19:51:41 +0000
Subject: [PATCH 04/29] Expose python arrays of mappings to C++. Fixes for
 expression interpolation + tests

---
 cpp/dolfinx/fem/Function.h                 |   3 +-
 python/dolfinx/fem/function.py             |   7 +-
 python/dolfinx/wrappers/mesh.cpp           |   6 +-
 python/test/unit/fem/test_interpolation.py | 110 ++++++++++++++++++++-
 4 files changed, 117 insertions(+), 9 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index d762b9bb08..9032a003b5 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -360,8 +360,7 @@ class Function
       std::runtime_error("Meshes are different and no cell map is provided");
     }
 
-    e.eval(*_function_space->mesh(), cells_expr, fdata,
-           {num_cells, num_points * value_size});
+    e.eval(expr_mesh, cells_expr, fdata, {num_cells, num_points * value_size});
 
     // Reshape evaluated data to fit interpolate
     // Expression returns matrix of shape (num_cells, num_points *
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 62202d2823..324f2abf86 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -466,9 +466,7 @@ def _(u_ptr: int, cells: typing.Optional[np.ndarray] = None):
             self._cpp_object.interpolate_ptr(u_ptr, cells)  # type: ignore
 
         @_interpolate.register(Expression)
-        def _(expr: Expression,
-              cells: typing.Optional[np.ndarray] = None
-              ):
+        def _(expr: Expression, cells: typing.Optional[np.ndarray] = None):
             """Interpolate Expression from a given mesh onto the set of cells
             Args:
                 expr: Expression to interpolate
@@ -486,8 +484,7 @@ def _(expr: Expression,
                 expr_cell_map = cell_map
                 assert expr_mesh is not None
                 mapping_mesh = expr_mesh._cpp_object
-            self._cpp_object.interpolate(expr._cpp_object, cells,
-                                         mapping_mesh, expr_cell_map)  # type: ignore
+            self._cpp_object.interpolate(expr._cpp_object, cells, mapping_mesh, expr_cell_map)  # type: ignore
 
         try:
             # u is a Function or Expression (or pointer to one)
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index 641c98639a..f0de46dfd3 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -301,8 +301,12 @@ void declare_mesh(nb::module_& m, std::string type)
       [](const dolfinx::mesh::Mesh<T>& mesh, int dim,
          nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> entities)
       {
-        return dolfinx::mesh::create_submesh(
+        auto [submesh, e_map, v_map, g_map] = dolfinx::mesh::create_submesh(
             mesh, dim, std::span(entities.data(), entities.size()));
+        auto _e_map = as_nbarray(e_map);
+        auto _v_map = as_nbarray(v_map);
+        auto _g_map = as_nbarray(g_map);
+        return std::tuple(submesh, _e_map, _v_map, _g_map);
       },
       nb::arg("mesh"), nb::arg("dim"), nb::arg("entities"));
 
diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 1834e6ce96..86055fe477 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -15,7 +15,7 @@
 import basix
 import ufl
 from basix.ufl import blocked_element, custom_element, element, enriched_element, mixed_element
-from dolfinx import default_real_type
+from dolfinx import default_real_type, default_scalar_type
 from dolfinx.fem import (
     Expression,
     Function,
@@ -29,6 +29,7 @@
     CellType,
     create_mesh,
     create_rectangle,
+    create_submesh,
     create_unit_cube,
     create_unit_square,
     locate_entities,
@@ -1025,3 +1026,110 @@ def f_test2(x):
 
     l2_error = assemble_scalar(form((u1 - u1_exact) ** 2 * dx_cell(cell_label), dtype=xtype))
     assert np.isclose(l2_error, 0.0, rtol=np.finfo(xtype).eps, atol=np.finfo(xtype).eps)
+
+
+def test_submesh_interpolation():
+    mesh = create_unit_square(MPI.COMM_WORLD, 6, 7)
+
+    def left_locator(x):
+        return x[0] <= 0.5 + 1e-14
+
+    def ref_func(x):
+        return x[0] + x[1] ** 2
+
+    tdim = mesh.topology.dim
+    cells = locate_entities(mesh, tdim, left_locator)
+    submesh, sub_to_parent, _, _ = create_submesh(mesh, tdim, cells)
+
+    V = functionspace(mesh, ("Lagrange", 2))
+    u = Function(V)
+    u.interpolate(ref_func)
+
+    V_sub = functionspace(submesh, ("DG", 3))
+    u_sub = Function(V_sub)
+
+    # Map from parent to sub mesh
+    u_sub.interpolate(u, cell_map=sub_to_parent)
+
+    u_sub_exact = Function(V_sub)
+    u_sub_exact.interpolate(ref_func)
+
+    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=1e-14)
+
+    # Map from sub to parent
+    W = functionspace(mesh, ("DG", 4))
+    w = Function(W)
+
+    cell_imap = mesh.topology.index_map(tdim)
+    num_cells = cell_imap.size_local + cell_imap.num_ghosts
+    parent_to_sub = np.full(num_cells, -1, dtype=np.int32)
+    parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
+
+    # Mapping back needs to be restricted to the subset of cells in the submesh
+    w.interpolate(u_sub, cells=sub_to_parent, cell_map=parent_to_sub)
+
+    w_exact = Function(W)
+    w_exact.interpolate(ref_func, cells=cells)
+
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=1e-14)
+
+
+def test_submesh_expression_interpolation():
+    mesh = create_unit_square(MPI.COMM_WORLD, 10, 8, cell_type=CellType.quadrilateral)
+
+    def left_locator(x):
+        return x[0] <= 0.5 + 1e-14
+
+    def ref_func(x):
+        return -3 * x[0] ** 2 + x[1] ** 2
+
+    def grad_ref_func(x):
+        values = np.zeros((2, x.shape[1]), dtype=default_scalar_type)
+        values[0] = -6 * x[0]
+        values[1] = 2 * x[1]
+        return values
+
+    def grad_squared(x):
+        grad = grad_ref_func(x)
+        return grad[0] ** 2 + grad[1] ** 2
+
+    tdim = mesh.topology.dim
+    cells = locate_entities(mesh, tdim, left_locator)
+    submesh, sub_to_parent, _, _ = create_submesh(mesh, tdim, cells)
+
+    V = functionspace(mesh, ("Lagrange", 2))
+    u = Function(V)
+    u.interpolate(ref_func)
+
+    V_sub = functionspace(submesh, ("N2curl", 1))
+    u_sub = Function(V_sub)
+
+    parent_expr = Expression(ufl.grad(u), V_sub.element.interpolation_points())
+
+    # Map from parent to sub mesh
+
+    u_sub.interpolate(parent_expr, expr_mesh=mesh, cell_map=sub_to_parent)
+
+    u_sub_exact = Function(V_sub)
+    u_sub_exact.interpolate(grad_ref_func)
+
+    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=1e-14)
+
+    # Map from sub to parent
+    W = functionspace(mesh, ("DQ", 2))
+    w = Function(W)
+
+    cell_imap = mesh.topology.index_map(tdim)
+    num_cells = cell_imap.size_local + cell_imap.num_ghosts
+    parent_to_sub = np.full(num_cells, -1, dtype=np.int32)
+    parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
+
+    sub_expr = Expression(ufl.dot(u_sub, u_sub), W.element.interpolation_points())
+
+    # Mapping back needs to be restricted to the subset of cells in the submesh
+    w.interpolate(sub_expr, cells=sub_to_parent, expr_mesh=submesh, cell_map=parent_to_sub)
+
+    w_exact = Function(W)
+    w_exact.interpolate(grad_squared, cells=cells)
+
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=1e-14)

From 62fd08f76abd8e0e60d2c34e821e67932eea3c59 Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 19:57:29 +0000
Subject: [PATCH 05/29] Add more asserts to please mypy

---
 python/dolfinx/fem/function.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 324f2abf86..356cd54256 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -473,8 +473,10 @@ def _(expr: Expression, cells: typing.Optional[np.ndarray] = None):
                 cells: The cells to interpolate over. If `None` then all
                     cells are interpolated over.
             """
+            assert cell_map is not None
             if len(cell_map) == 0:
                 # If cell map is not provided create identity map
+                assert cells is not None
                 expr_cell_map = np.arange(len(cells), dtype=np.int32)
                 assert expr_mesh is None
                 mapping_mesh = self.function_space.mesh._cpp_object

From 8232486d73d5aa9abc16ba0c98e261a1f53c7882 Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 20:14:17 +0000
Subject: [PATCH 06/29] Fix C++ interface

---
 cpp/demo/interpolation_different_meshes/main.cpp | 3 ++-
 cpp/dolfinx/fem/Function.h                       | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index 5dbcb4d895..09651b6da4 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -74,7 +74,8 @@ int main(int argc, char* argv[])
             *u_hex->function_space()->mesh(),
             *u_hex->function_space()->element(),
             *u_tet->function_space()->mesh(), padding);
-    u_hex->interpolate(*u_tet, nmm_interpolation_data);
+    constexpr std::span<const std::int32_t> cell_map;
+    u_hex->interpolate(*u_tet, cell_map, nmm_interpolation_data);
 
 #ifdef HAS_ADIOS2
     io::VTXWriter<double> write_tet(mesh_tet->comm(), "u_tet.bp", {u_tet});
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 9032a003b5..3688b41f51 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -178,7 +178,7 @@ class Function
   /// generate_nonmatching_meshes_interpolation_data (optional).
   void interpolate(
       const Function<value_type, geometry_type>& v,
-      std::span<const std::int32_t> cell_map,
+      std::span<const std::int32_t> cell_map = std::span<const std::int32_t>(),
       const std::tuple<std::span<const std::int32_t>,
                        std::span<const std::int32_t>,
                        std::span<const geometry_type>,

From 62a49fa3798984f03548e9093b42b5d046e8f61e Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 20:26:23 +0000
Subject: [PATCH 07/29] Set some default zero-spans in C++ + fix
 hyperelasticity demo

---
 cpp/demo/hyperelasticity/main.cpp | 2 +-
 cpp/dolfinx/fem/Function.h        | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index e286f2bf39..70dee61254 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -258,7 +258,7 @@ int main(int argc, char* argv[])
 
     auto sigma = fem::Function<T>(S);
     sigma.name = "cauchy_stress";
-    sigma.interpolate(sigma_expression);
+    sigma.interpolate(sigma_expression, *mesh);
 
     // Save solution in VTK format
     io::VTKFile file_u(mesh->comm(), "u.pvd", "w");
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 3688b41f51..f0b2f6ee64 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -292,7 +292,8 @@ class Function
   void interpolate(const Expression<value_type, geometry_type>& e,
                    std::span<const std::int32_t> cells,
                    const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
-                   std::span<const std::int32_t> cell_map)
+                   std::span<const std::int32_t> cell_map
+                   = std::span<const std::int32_t>())
   {
     // Check that spaces are compatible
     assert(_function_space);
@@ -389,7 +390,8 @@ class Function
   /// receiving function is defined on a different mesh than the expression
   void interpolate(const Expression<value_type, geometry_type>& e,
                    const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
-                   std::span<const std::int32_t> cell_map)
+                   std::span<const std::int32_t> cell_map
+                   = std::span<const std::int32_t>())
   {
     assert(_function_space);
     assert(_function_space->mesh());

From 073902721f8d3f852411c69c74842db1a461eb3f Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 20:40:59 +0000
Subject: [PATCH 08/29] Set atol depending on default floating type

---
 python/test/unit/fem/test_interpolation.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 86055fe477..fc64934f3c 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1053,8 +1053,8 @@ def ref_func(x):
 
     u_sub_exact = Function(V_sub)
     u_sub_exact.interpolate(ref_func)
-
-    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=1e-14)
+    atol = 5 * np.finfo(default_scalar_type).resolution
+    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=atol)
 
     # Map from sub to parent
     W = functionspace(mesh, ("DG", 4))
@@ -1071,7 +1071,7 @@ def ref_func(x):
     w_exact = Function(W)
     w_exact.interpolate(ref_func, cells=cells)
 
-    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=1e-14)
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)
 
 
 def test_submesh_expression_interpolation():
@@ -1112,8 +1112,8 @@ def grad_squared(x):
 
     u_sub_exact = Function(V_sub)
     u_sub_exact.interpolate(grad_ref_func)
-
-    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=1e-14)
+    atol = 5 * np.finfo(default_scalar_type).resolution
+    np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=atol)
 
     # Map from sub to parent
     W = functionspace(mesh, ("DQ", 2))
@@ -1131,5 +1131,4 @@ def grad_squared(x):
 
     w_exact = Function(W)
     w_exact.interpolate(grad_squared, cells=cells)
-
-    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=1e-14)
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)

From 7111ed1f948999bc2cf966206783a8a5f8464f43 Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 20:59:02 +0000
Subject: [PATCH 09/29] Use exact solution on submesh when interpolating back
 to avoid accumulation of round-off errors

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index fc64934f3c..e89ac5f89d 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1066,7 +1066,7 @@ def ref_func(x):
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
     # Mapping back needs to be restricted to the subset of cells in the submesh
-    w.interpolate(u_sub, cells=sub_to_parent, cell_map=parent_to_sub)
+    w.interpolate(u_sub_exact, cells=sub_to_parent, cell_map=parent_to_sub)
 
     w_exact = Function(W)
     w_exact.interpolate(ref_func, cells=cells)

From 71f9adfd1933279fa945db4e19953145f32d5601 Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 21:12:49 +0000
Subject: [PATCH 10/29] Use exact solution in mapping back due to same issue

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index e89ac5f89d..6677062081 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1124,7 +1124,7 @@ def grad_squared(x):
     parent_to_sub = np.full(num_cells, -1, dtype=np.int32)
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
-    sub_expr = Expression(ufl.dot(u_sub, u_sub), W.element.interpolation_points())
+    sub_expr = Expression(ufl.dot(u_sub_exact, u_sub_exact), W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh
     w.interpolate(sub_expr, cells=sub_to_parent, expr_mesh=submesh, cell_map=parent_to_sub)

From af3855910384b51cefedee98c5a9ab505db4780f Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 24 Mar 2024 21:25:29 +0000
Subject: [PATCH 11/29] Less strict tolerance

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 6677062081..76a621d9d9 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1131,4 +1131,4 @@ def grad_squared(x):
 
     w_exact = Function(W)
     w_exact.interpolate(grad_squared, cells=cells)
-    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=5e-4)

From 2fb6ee891c60baa741735c857e7eeb02cd3540af Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Mon, 8 Apr 2024 16:35:51 +0200
Subject: [PATCH 12/29] Apply suggestions from code review

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>
---
 cpp/dolfinx/fem/Function.h    | 4 ++--
 cpp/dolfinx/fem/interpolate.h | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index f0b2f6ee64..0d922812be 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -351,7 +351,7 @@ class Function
       cells_expr.insert(cells_expr.end(), cells.begin(), cells.end());
     }
     // If meshes are different and input mapping is given
-    else if (cell_map.size() > 0)
+    else if (!cell_map.empty())
     {
       std::transform(cells.begin(), cells.end(), std::back_inserter(cells_expr),
                      [&cell_map](std::int32_t c) { return cell_map[c]; });
@@ -385,7 +385,7 @@ class Function
 
   /// Interpolate an Expression (based on UFL) on all cells
   /// @param[in] e The function to be interpolated
-  /// @param[in] expr_mesh Mesh expression is defined on
+  /// @param[in] mesh Mesh the Expression `e` is defined on.
   /// @param[in] cell_map Map from `cells` to cells in expression if
   /// receiving function is defined on a different mesh than the expression
   void interpolate(const Expression<value_type, geometry_type>& e,
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 0da8877c19..378cb1ac25 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -1160,7 +1160,7 @@ void interpolate(
                      [&cell_map](std::int32_t c) { return cell_map[c]; });
     }
     // Non-matching meshes
-    if (cells_v.size() == 0)
+    if (cells_v.empty())
     {
       impl::interpolate_nonmatching_meshes(u, v, cells_u,
                                            nmm_interpolation_data);

From eaa9e432a80d279dd8d123ead6bda003acd189c1 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 14:46:57 +0000
Subject: [PATCH 13/29] Add default arguments in C++ and add documentation to
 tests

---
 cpp/dolfinx/fem/Function.h                 |  7 ++--
 cpp/dolfinx/fem/interpolate.h              | 37 +++++++++++-----------
 python/test/unit/fem/test_interpolation.py |  2 ++
 3 files changed, 24 insertions(+), 22 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 0d922812be..80ac7e22c7 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -178,7 +178,8 @@ class Function
   /// generate_nonmatching_meshes_interpolation_data (optional).
   void interpolate(
       const Function<value_type, geometry_type>& v,
-      std::span<const std::int32_t> cell_map = std::span<const std::int32_t>(),
+      std::span<const std::int32_t> cell_map = std::span<const std::int32_t>()
+      = {},
       const std::tuple<std::span<const std::int32_t>,
                        std::span<const std::int32_t>,
                        std::span<const geometry_type>,
@@ -357,9 +358,7 @@ class Function
                      [&cell_map](std::int32_t c) { return cell_map[c]; });
     }
     else
-    {
       std::runtime_error("Meshes are different and no cell map is provided");
-    }
 
     e.eval(expr_mesh, cells_expr, fdata, {num_cells, num_points * value_size});
 
@@ -391,7 +390,7 @@ class Function
   void interpolate(const Expression<value_type, geometry_type>& e,
                    const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
                    std::span<const std::int32_t> cell_map
-                   = std::span<const std::int32_t>())
+                   = std::span<const std::int32_t>() = {})
   {
     assert(_function_space);
     assert(_function_space->mesh());
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 378cb1ac25..68a0516aba 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -321,15 +321,15 @@ void interpolation_apply(U&& Pi, V&& data, std::span<T> coeffs, int bs)
 /// map.
 /// @param[out] u1 The function to interpolate to
 /// @param[in] u0 The function to interpolate from
-/// @param[in] cells_u1 The cells to interpolate on
-/// @param[in] cells_u0 Equivalent cell in u0 for each cell in u1
+/// @param[in] cells1 The cells to interpolate on
+/// @param[in] cells0 Equivalent cell in u0 for each cell in u1
 /// @pre The functions `u1` and `u0` must share the same mesh and the
 /// elements must share the same basis function map. Neither is checked
 /// by the function.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
-                          std::span<const std::int32_t> cells_u1,
-                          std::span<const std::int32_t> cells_u0)
+                          std::span<const std::int32_t> cells1,
+                          std::span<const std::int32_t> cells0)
 {
   auto V0 = u0.function_space();
   assert(V0);
@@ -346,6 +346,7 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
   auto element1 = V1->element();
   assert(element1);
 
+  assert(mesh0->topology()->dim());
   const int tdim = mesh0->topology()->dim();
   auto map = mesh0->topology()->index_map(tdim);
   assert(map);
@@ -387,14 +388,14 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
 
   // Iterate over mesh and interpolate on each cell
   using X = typename dolfinx::scalar_value_type_t<T>;
-  for (std::size_t c = 0; c < cells_u0.size(); c++)
+  for (std::size_t c = 0; c < cells0.size(); c++)
   {
-    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_u0[c]);
+    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells0[c]);
     for (std::size_t i = 0; i < dofs0.size(); ++i)
       for (int k = 0; k < bs0; ++k)
         local0[bs0 * i + k] = u0_array[bs0 * dofs0[i] + k];
 
-    apply_dof_transformation(local0, cell_info0, cells_u0[c], 1);
+    apply_dof_transformation(local0, cell_info0, cells0[c], 1);
 
     // FIXME: Get compile-time ranges from Basix
     // Apply interpolation operator
@@ -403,8 +404,8 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
       for (std::size_t j = 0; j < im_shape[1]; ++j)
         local1[i] += static_cast<X>(i_m[im_shape[1] * i + j]) * local0[j];
 
-    apply_inverse_dof_transform(local1, cell_info1, cells_u1[c], 1);
-    std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells_u1[c]);
+    apply_inverse_dof_transform(local1, cell_info1, cells1[c], 1);
+    std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells1[c]);
     for (std::size_t i = 0; i < dofs1.size(); ++i)
       for (int k = 0; k < bs1; ++k)
         u1_array[bs1 * dofs1[i] + k] = local1[bs1 * i + k];
@@ -418,14 +419,14 @@ void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
 /// standard isoparametric mapping.
 /// @param[out] u1 The function to interpolate to
 /// @param[in] u0 The function to interpolate from
-/// @param[in] cells_u1 The cells to interpolate on
-/// @param[in] cells_u0 Equivalent cell in u0 for each cell in u1
+/// @param[in] cells1 The cells to interpolate on
+/// @param[in] cells0 Equivalent cell in u0 for each cell in u1
 /// @pre The functions `u1` and `u0` must share the same mesh. This is
 /// not checked by the function.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
-                                  std::span<const std::int32_t> cells_u1,
-                                  std::span<const std::int32_t> cells_u0)
+                                  std::span<const std::int32_t> cells1,
+                                  std::span<const std::int32_t> cells0)
 {
   // Get mesh
   auto V0 = u0.function_space();
@@ -552,11 +553,11 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
   // Iterate over mesh and interpolate on each cell
   std::span<const T> array0 = u0.x()->array();
   std::span<T> array1 = u1.x()->mutable_array();
-  for (std::size_t c = 0; c < cells_u0.size(); c++)
+  for (std::size_t c = 0; c < cells0.size(); c++)
   {
     // Get cell geometry (coordinate dofs)
     auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
-        x_dofmap, cells_u0[c], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
+        x_dofmap, cells0[c], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
     for (std::size_t i = 0; i < num_dofs_g; ++i)
     {
       const int pos = 3 * x_dofs[i];
@@ -596,7 +597,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
       apply_dof_transformation0(
           std::span(basis_reference0_b.data() + p * dim0 * value_size_ref0,
                     dim0 * value_size_ref0),
-          cell_info0, cells_u0[c], value_size_ref0);
+          cell_info0, cells0[c], value_size_ref0);
     }
 
     for (std::size_t i = 0; i < basis0.extent(0); ++i)
@@ -618,7 +619,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
 
     // Copy expansion coefficients for v into local array
     const int dof_bs0 = dofmap0->bs();
-    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_u0[c]);
+    std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells0[c]);
     for (std::size_t i = 0; i < dofs0.size(); ++i)
       for (int k = 0; k < dof_bs0; ++k)
         coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k];
@@ -661,7 +662,7 @@ void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
         mapped_values0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0,
         MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
     interpolation_apply(Pi_1, values, std::span(local1), bs1);
-    apply_inverse_dof_transform1(local1, cell_info1, cells_u1[c], 1);
+    apply_inverse_dof_transform1(local1, cell_info1, cells1[c], 1);
 
     // Copy local coefficients to the correct position in u dof array
     const int dof_bs1 = dofmap1->bs();
diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 76a621d9d9..0220429c03 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1029,6 +1029,7 @@ def f_test2(x):
 
 
 def test_submesh_interpolation():
+    """Test interpolation of a function between a submesh and its parent"""
     mesh = create_unit_square(MPI.COMM_WORLD, 6, 7)
 
     def left_locator(x):
@@ -1075,6 +1076,7 @@ def ref_func(x):
 
 
 def test_submesh_expression_interpolation():
+    """Test interpolation of an expression between a submesh and its parent"""
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 8, cell_type=CellType.quadrilateral)
 
     def left_locator(x):

From 46f0118a05c0edacdda1ee852e4b6f64dd33fed5 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 14:54:32 +0000
Subject: [PATCH 14/29] Fix submesh C++ wrapping

---
 python/dolfinx/wrappers/mesh.cpp | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index f0de46dfd3..6b80317196 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -301,12 +301,13 @@ void declare_mesh(nb::module_& m, std::string type)
       [](const dolfinx::mesh::Mesh<T>& mesh, int dim,
          nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> entities)
       {
-        auto [submesh, e_map, v_map, g_map] = dolfinx::mesh::create_submesh(
+        auto submesh = dolfinx::mesh::create_submesh(
             mesh, dim, std::span(entities.data(), entities.size()));
-        auto _e_map = as_nbarray(e_map);
-        auto _v_map = as_nbarray(v_map);
-        auto _g_map = as_nbarray(g_map);
-        return std::tuple(submesh, _e_map, _v_map, _g_map);
+        auto _e_map = as_nbarray(std::move(std::get<1>(submesh)));
+        auto _v_map = as_nbarray(std::move(std::get<2>(submesh)));
+        auto _g_map = as_nbarray(std::move(std::get<3>(submesh)));
+        return std::tuple(std::move(std::get<0>(submesh)), _e_map, _v_map,
+                          _g_map);
       },
       nb::arg("mesh"), nb::arg("dim"), nb::arg("entities"));
 

From 72579a92c312e0561747663e1b234b031d93436b Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 15:03:18 +0000
Subject: [PATCH 15/29] Renaming + updating docstring

---
 cpp/dolfinx/fem/interpolate.h | 58 ++++++++++++++++++-----------------
 1 file changed, 30 insertions(+), 28 deletions(-)

diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 68a0516aba..a03ead0419 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -1112,17 +1112,19 @@ create_nonmatching_meshes_interpolation_data(const mesh::Mesh<T>& mesh0,
 }
 
 /// @brief Interpolate from one finite element Function to another one.
-/// @param[out] u The function to interpolate into
-/// @param[in] v The function to be interpolated
-/// @param[in] cells_u List of cell indices to interpolate on
-/// @param[in] cell_map Mapping from cells in `u` to cells in `v`
+/// @param[out] u1 The function to interpolate into
+/// @param[in] u0 The function to be interpolated
+/// @param[in] cells1 List of cell indices associated with the mesh of `u1` that
+/// will be interpolated onto
+/// @param[in] cell_map Mapping of cells in mesh associated with `u1` to cells
+/// associated with `u0`
 /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
 /// nonmatching meshes. This data can be generated with
 /// create_nonmatching_meshes_interpolation_data (optional).
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate(
-    Function<T, U>& u, const Function<T, U>& v,
-    std::span<const std::int32_t> cells_u,
+    Function<T, U>& u1, const Function<T, U>& u0,
+    std::span<const std::int32_t> cells1,
     std::span<const std::int32_t> cell_map,
     const std::tuple<std::span<const std::int32_t>,
                      std::span<const std::int32_t>, std::span<const U>,
@@ -1130,49 +1132,49 @@ void interpolate(
     = {})
 {
   assert(u.function_space());
-  assert(v.function_space());
-  auto mesh = u.function_space()->mesh();
+  assert(u0.function_space());
+  auto mesh = u1.function_space()->mesh();
   assert(mesh);
 
   auto cell_map0 = mesh->topology()->index_map(mesh->topology()->dim());
   assert(cell_map0);
   std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
-  if (u.function_space() == v.function_space() and cells_u.size() == num_cells0)
+  if (u1.function_space() == u0.function_space()
+      and cells1.size() == num_cells0)
   {
     // Same function spaces and on whole mesh
-    std::span<T> u1_array = u.x()->mutable_array();
-    std::span<const T> u0_array = v.x()->array();
+    std::span<T> u1_array = u1.x()->mutable_array();
+    std::span<const T> u0_array = u0.x()->array();
     std::copy(u0_array.begin(), u0_array.end(), u1_array.begin());
   }
   else
   {
-    std::vector<std::int32_t> cells_v;
-    cells_v.reserve(cells_u.size());
+    std::vector<std::int32_t> cells0;
+    cells0.reserve(cells1.size());
     // Get mesh and check that functions share the same mesh
-    if (auto mesh_v = v.function_space()->mesh(); mesh == mesh_v)
+    if (auto mesh_v = u0.function_space()->mesh(); mesh == mesh_v)
     {
-      cells_v.insert(cells_v.end(), cells_u.begin(), cells_u.end());
+      cells0.insert(cells0.end(), cells1.begin(), cells1.end());
     }
     // If meshes are different and input mapping is given
     else if (cell_map.size() > 0)
     {
-      std::transform(cells_u.begin(), cells_u.end(),
-                     std::back_inserter(cells_v),
+      std::transform(cells1.begin(), cells1.end(), std::back_inserter(cells0),
                      [&cell_map](std::int32_t c) { return cell_map[c]; });
     }
     // Non-matching meshes
-    if (cells_v.empty())
+    if (cells0.empty())
     {
-      impl::interpolate_nonmatching_meshes(u, v, cells_u,
+      impl::interpolate_nonmatching_meshes(u, u0, cells1,
                                            nmm_interpolation_data);
       return;
     }
 
     // Get elements and check value shape
-    auto fs0 = v.function_space();
+    auto fs0 = u0.function_space();
     auto element0 = fs0->element();
     assert(element0);
-    auto fs1 = u.function_space();
+    auto fs1 = u1.function_space();
     auto element1 = fs1->element();
     assert(element1);
     if (fs0->value_shape().size() != fs1->value_shape().size()
@@ -1194,9 +1196,9 @@ void interpolate(
       assert(element1->block_size() == element0->block_size());
 
       // Get dofmaps
-      std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
+      std::shared_ptr<const DofMap> dofmap0 = u0.function_space()->dofmap();
       assert(dofmap0);
-      std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
+      std::shared_ptr<const DofMap> dofmap1 = u1.function_space()->dofmap();
       assert(dofmap1);
 
       std::span<T> u1_array = u.x()->mutable_array();
@@ -1205,10 +1207,10 @@ void interpolate(
       // Iterate over mesh and interpolate on each cell
       const int bs0 = dofmap0->bs();
       const int bs1 = dofmap1->bs();
-      for (std::size_t c = 0; c < cells_u.size(); ++c)
+      for (std::size_t c = 0; c < cells1.size(); ++c)
       {
-        std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells_v[c]);
-        std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells_u[c]);
+        std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cells0[c]);
+        std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cells1[c]);
         assert(bs0 * dofs0.size() == bs1 * dofs1.size());
         for (std::size_t i = 0; i < dofs0.size(); ++i)
         {
@@ -1225,12 +1227,12 @@ void interpolate(
     else if (element1->map_type() == element0->map_type())
     {
       // Different elements, same basis function map type
-      impl::interpolate_same_map(u, v, cells_u, cells_v);
+      impl::interpolate_same_map(u1, u0, cells1, cells0);
     }
     else
     {
       //  Different elements with different maps for basis functions
-      impl::interpolate_nonmatching_maps(u, v, cells_u, cells_v);
+      impl::interpolate_nonmatching_maps(u1, u0, cells1, cells0);
     }
   }
 }

From 45cf1ed3449d3702ab3eca2d2ef1df446f8949c1 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 15:06:07 +0000
Subject: [PATCH 16/29] Fix conversion of u

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

diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index a03ead0419..0061818f9d 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -1131,7 +1131,7 @@ void interpolate(
                      std::span<const std::int32_t>>& nmm_interpolation_data
     = {})
 {
-  assert(u.function_space());
+  assert(u1.function_space());
   assert(u0.function_space());
   auto mesh = u1.function_space()->mesh();
   assert(mesh);

From d176828b7846f358027af25ce25d7f7d9618d73d Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 15:11:51 +0000
Subject: [PATCH 17/29] More renaming

---
 cpp/dolfinx/fem/interpolate.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 0061818f9d..6ddafaa4cc 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -1165,7 +1165,7 @@ void interpolate(
     // Non-matching meshes
     if (cells0.empty())
     {
-      impl::interpolate_nonmatching_meshes(u, u0, cells1,
+      impl::interpolate_nonmatching_meshes(u1, u0, cells1,
                                            nmm_interpolation_data);
       return;
     }
@@ -1201,8 +1201,8 @@ void interpolate(
       std::shared_ptr<const DofMap> dofmap1 = u1.function_space()->dofmap();
       assert(dofmap1);
 
-      std::span<T> u1_array = u.x()->mutable_array();
-      std::span<const T> u0_array = v.x()->array();
+      std::span<T> u1_array = u1.x()->mutable_array();
+      std::span<const T> u0_array = u0.x()->array();
 
       // Iterate over mesh and interpolate on each cell
       const int bs0 = dofmap0->bs();

From 460acb0b25e83da44ab3c86c4110a843543f0990 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 8 Apr 2024 15:16:56 +0000
Subject: [PATCH 18/29] Fix documentation

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

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 80ac7e22c7..b66cb75f31 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -384,7 +384,7 @@ class Function
 
   /// Interpolate an Expression (based on UFL) on all cells
   /// @param[in] e The function to be interpolated
-  /// @param[in] mesh Mesh the Expression `e` is defined on.
+  /// @param[in] expr_mesh Mesh the Expression `e` is defined on.
   /// @param[in] cell_map Map from `cells` to cells in expression if
   /// receiving function is defined on a different mesh than the expression
   void interpolate(const Expression<value_type, geometry_type>& e,

From 2f0679598f32c6ed5fbd4591ff38b95fae7b6656 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Wed, 10 Apr 2024 14:44:25 +0200
Subject: [PATCH 19/29] Apply suggestions from code review

Remove empty comment
---
 cpp/dolfinx/fem/Function.h | 1 -
 1 file changed, 1 deletion(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index b66cb75f31..9dcc181d8c 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -342,7 +342,6 @@ class Function
     // Evaluate Expression at points
     assert(_function_space->mesh());
 
-    //
     std::vector<std::int32_t> cells_expr;
     cells_expr.reserve(num_cells);
     // Get mesh and check if mesh is shared

From f8052ed97db0550eb61915ac470a4f0c0d51f82f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Sun, 14 Apr 2024 09:57:17 +0200
Subject: [PATCH 20/29] Improve documentation

Co-authored-by: Joe Dean <jpd62@cam.ac.uk>
---
 cpp/dolfinx/fem/Function.h     | 4 ++--
 python/dolfinx/fem/function.py | 5 ++++-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 9dcc181d8c..4faf83c589 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -153,7 +153,7 @@ class Function
   /// @brief Interpolate a provided Function.
   /// @param[in] v The function to be interpolated
   /// @param[in] cells The cells to interpolate on
-  /// @param[in] cell_map Map from indices in `cells` to cell indices in `v`
+  /// @param[in] cell_map A map from cells in the mesh associated with `this` function to cells in mesh associated with `v`
   /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
   /// nonmatching meshes. This data can be generated with
   /// generate_nonmatching_meshes_interpolation_data (optional).
@@ -288,7 +288,7 @@ class Function
   /// `FiniteElement::interpolation_points()` for the element associated
   /// with `u`.
   /// @param[in] cells The cells to interpolate on
-  /// @param[in] expr_mesh Mesh expression is defined on
+  /// @param[in] expr_mesh The mesh to evaluate the expression on
   /// @param[in] cell_map Map from `cells` to cells in expression
   void interpolate(const Expression<value_type, geometry_type>& e,
                    std::span<const std::int32_t> cells,
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index d729cdfb43..6acb45b171 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -426,7 +426,10 @@ def interpolate(
             u: The function, Expression or Function to interpolate.
             cells: The cells to interpolate over. If `None` then all
                 cells are interpolated over.
-            cell_map: Mapping from cells in receiving mesh to cells in `u` mesh
+            cell_map: Mapping from `cells` to to cells in the mesh that `u` is defined over.
+            expr_mesh: If an Expression with coefficients or constants from another mesh
+                than the function is supplied, the mesh associated with this expression has 
+                to be provided, along with `cell_map.`
             nmm_interpolation_data: Data needed to interpolate functions defined on other meshes
         """
         if nmm_interpolation_data is None:

From c5df0d6475a6f743da4e771d11c7eb03e9dbf5bc Mon Sep 17 00:00:00 2001
From: jorgensd <dokken92@gmail.com>
Date: Sun, 14 Apr 2024 10:00:09 +0200
Subject: [PATCH 21/29] Ruff format docstring

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

diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 6acb45b171..25eecbbe7f 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -428,7 +428,7 @@ def interpolate(
                 cells are interpolated over.
             cell_map: Mapping from `cells` to to cells in the mesh that `u` is defined over.
             expr_mesh: If an Expression with coefficients or constants from another mesh
-                than the function is supplied, the mesh associated with this expression has 
+                than the function is supplied, the mesh associated with this expression has
                 to be provided, along with `cell_map.`
             nmm_interpolation_data: Data needed to interpolate functions defined on other meshes
         """

From cce45bd36b0a0762ddee6d2930f459a17574bca1 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 16:53:25 +0000
Subject: [PATCH 22/29] Try less strict atol

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 0220429c03..272bf72679 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1133,4 +1133,5 @@ def grad_squared(x):
 
     w_exact = Function(W)
     w_exact.interpolate(grad_squared, cells=cells)
-    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=5e-4)
+    atol = 5 * np.finfo(default_scalar_type).resolution
+    np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)

From a1a57aca2e6de83f1efcf479859328d1ce49f1e7 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 17:26:14 +0000
Subject: [PATCH 23/29] Modify test to avoid squaring functions

---
 python/test/unit/fem/test_interpolation.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 272bf72679..b14dbc8ca3 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1091,9 +1091,10 @@ def grad_ref_func(x):
         values[1] = 2 * x[1]
         return values
 
-    def grad_squared(x):
+
+    def modified_grad(x):
         grad = grad_ref_func(x)
-        return grad[0] ** 2 + grad[1] ** 2
+        return -3 * grad[1], 2 * grad[0]
 
     tdim = mesh.topology.dim
     cells = locate_entities(mesh, tdim, left_locator)
@@ -1114,11 +1115,11 @@ def grad_squared(x):
 
     u_sub_exact = Function(V_sub)
     u_sub_exact.interpolate(grad_ref_func)
-    atol = 5 * np.finfo(default_scalar_type).resolution
+    atol = 6 * np.finfo(default_scalar_type).resolution
     np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=atol)
 
     # Map from sub to parent
-    W = functionspace(mesh, ("DQ", 2))
+    W = functionspace(mesh, ("DQ", 2, (mesh.geometry.dim,)))
     w = Function(W)
 
     cell_imap = mesh.topology.index_map(tdim)
@@ -1126,12 +1127,13 @@ def grad_squared(x):
     parent_to_sub = np.full(num_cells, -1, dtype=np.int32)
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
-    sub_expr = Expression(ufl.dot(u_sub_exact, u_sub_exact), W.element.interpolation_points())
+    # Map exact solution (based on quadrature points) back to parent mesh
+    sub_vec = ufl.as_vector((-3* u_sub_exact[1], 2 * u_sub_exact[0]))
+    sub_expr = Expression(sub_vec, W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh
     w.interpolate(sub_expr, cells=sub_to_parent, expr_mesh=submesh, cell_map=parent_to_sub)
 
     w_exact = Function(W)
-    w_exact.interpolate(grad_squared, cells=cells)
-    atol = 5 * np.finfo(default_scalar_type).resolution
+    w_exact.interpolate(modified_grad, cells=cells)
     np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)

From 0d390f25a1c7f455966e8b0af4400aa4e3d8b9e4 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 17:30:24 +0000
Subject: [PATCH 24/29] Ruff formatting

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index b14dbc8ca3..1bb2648859 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1091,7 +1091,6 @@ def grad_ref_func(x):
         values[1] = 2 * x[1]
         return values
 
-
     def modified_grad(x):
         grad = grad_ref_func(x)
         return -3 * grad[1], 2 * grad[0]
@@ -1128,7 +1127,7 @@ def modified_grad(x):
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
     # Map exact solution (based on quadrature points) back to parent mesh
-    sub_vec = ufl.as_vector((-3* u_sub_exact[1], 2 * u_sub_exact[0]))
+    sub_vec = ufl.as_vector((-3 * u_sub_exact[1], 2 * u_sub_exact[0]))
     sub_expr = Expression(sub_vec, W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh

From 2a04e4df4219ec3185d45a6e6ef4c49407568bdb Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 17:51:11 +0000
Subject: [PATCH 25/29] Scatter forward data

---
 python/test/unit/fem/test_interpolation.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 1bb2648859..cca7692652 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1135,4 +1135,5 @@ def modified_grad(x):
 
     w_exact = Function(W)
     w_exact.interpolate(modified_grad, cells=cells)
+    w_exact.x.scatter_forward()
     np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)

From a908db1d3747833d1887f9b6aa76ba644bf9a996 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 17:57:17 +0000
Subject: [PATCH 26/29] Modify test again

---
 python/test/unit/fem/test_interpolation.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index cca7692652..038c4d3352 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1093,7 +1093,7 @@ def grad_ref_func(x):
 
     def modified_grad(x):
         grad = grad_ref_func(x)
-        return -3 * grad[1], 2 * grad[0]
+        return -grad[1], 2 * grad[0]
 
     tdim = mesh.topology.dim
     cells = locate_entities(mesh, tdim, left_locator)
@@ -1118,7 +1118,7 @@ def modified_grad(x):
     np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=atol)
 
     # Map from sub to parent
-    W = functionspace(mesh, ("DQ", 2, (mesh.geometry.dim,)))
+    W = functionspace(mesh, ("DQ", 1, (mesh.geometry.dim,)))
     w = Function(W)
 
     cell_imap = mesh.topology.index_map(tdim)
@@ -1127,7 +1127,7 @@ def modified_grad(x):
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
     # Map exact solution (based on quadrature points) back to parent mesh
-    sub_vec = ufl.as_vector((-3 * u_sub_exact[1], 2 * u_sub_exact[0]))
+    sub_vec = ufl.as_vector((-u_sub_exact[1], 2 * u_sub_exact[0]))
     sub_expr = Expression(sub_vec, W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh

From 4cd78fc58a285b4413fa2f62616467efdb598b5f Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 18:28:48 +0000
Subject: [PATCH 27/29] Try scaling down gradient

---
 python/test/unit/fem/test_interpolation.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 038c4d3352..75f9b7f22c 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1093,7 +1093,7 @@ def grad_ref_func(x):
 
     def modified_grad(x):
         grad = grad_ref_func(x)
-        return -grad[1], 2 * grad[0]
+        return -0.2 * grad[1], 0.1 * grad[0]
 
     tdim = mesh.topology.dim
     cells = locate_entities(mesh, tdim, left_locator)
@@ -1114,7 +1114,7 @@ def modified_grad(x):
 
     u_sub_exact = Function(V_sub)
     u_sub_exact.interpolate(grad_ref_func)
-    atol = 6 * np.finfo(default_scalar_type).resolution
+    atol = 10 * np.finfo(default_scalar_type).resolution
     np.testing.assert_allclose(u_sub_exact.x.array, u_sub.x.array, atol=atol)
 
     # Map from sub to parent
@@ -1127,7 +1127,7 @@ def modified_grad(x):
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
     # Map exact solution (based on quadrature points) back to parent mesh
-    sub_vec = ufl.as_vector((-u_sub_exact[1], 2 * u_sub_exact[0]))
+    sub_vec = ufl.as_vector((-0.2 * u_sub_exact[1], -0.5 * u_sub_exact[0]))
     sub_expr = Expression(sub_vec, W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh

From 0d8193349802cbc00827cba34b1ebf2b2915a2e6 Mon Sep 17 00:00:00 2001
From: "Jorgen S. Dokken" <dokken@simula.no>
Date: Mon, 22 Apr 2024 18:44:30 +0000
Subject: [PATCH 28/29] Fix test

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

diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 75f9b7f22c..5a9d93205c 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -1127,7 +1127,7 @@ def modified_grad(x):
     parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
 
     # Map exact solution (based on quadrature points) back to parent mesh
-    sub_vec = ufl.as_vector((-0.2 * u_sub_exact[1], -0.5 * u_sub_exact[0]))
+    sub_vec = ufl.as_vector((-0.2 * u_sub_exact[1], 0.1 * u_sub_exact[0]))
     sub_expr = Expression(sub_vec, W.element.interpolation_points())
 
     # Mapping back needs to be restricted to the subset of cells in the submesh

From 4953137fef932e64a1984347c986b576f9130d9b Mon Sep 17 00:00:00 2001
From: "Jack S. Hale" <mail@jackhale.co.uk>
Date: Tue, 23 Apr 2024 09:05:17 +0200
Subject: [PATCH 29/29] clang-format

---
 cpp/dolfinx/fem/Function.h    | 3 ++-
 cpp/dolfinx/fem/interpolate.h | 8 ++------
 2 files changed, 4 insertions(+), 7 deletions(-)

diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 1d749e2fcb..050b4892db 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -153,7 +153,8 @@ class Function
   /// @brief Interpolate a provided Function.
   /// @param[in] v The function to be interpolated
   /// @param[in] cells The cells to interpolate on
-  /// @param[in] cell_map A map from cells in the mesh associated with `this` function to cells in mesh associated with `v`
+  /// @param[in] cell_map A map from cells in the mesh associated with `this`
+  /// function to cells in mesh associated with `v`
   /// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
   /// nonmatching meshes. This data can be generated with
   /// generate_nonmatching_meshes_interpolation_data (optional).
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index 53f4abe9b5..239c3384c6 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -31,13 +31,9 @@ class Function;
 template <typename T>
 concept MDSpan = requires(T x, std::size_t idx) {
   x(idx, idx);
-  {
-    x.extent(0)
-  } -> std::integral;
+  { x.extent(0) } -> std::integral;
 
-  {
-    x.extent(1)
-  } -> std::integral;
+  { x.extent(1) } -> std::integral;
 };
 
 /// @brief Compute the evaluation points in the physical space at which