Interpolation to and from sub-meshes for co-dimension 0

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");

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});

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 = 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>,
@@ -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,13 @@ 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<geometry_type>& expr_mesh,
+                   std::span<const std::int32_t> cell_map
+                   = std::span<const std::int32_t>())
   {
     // Check that spaces are compatible
     assert(_function_space);
@@ -331,8 +340,28 @@ class Function
 
     // Evaluate Expression at points
     assert(_function_space->mesh());
-    e.eval(*_function_space->mesh(), cells, fdata,
-           {num_cells, num_points * value_size});
+
+    //
+    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.topology() == mesh_v->topology())
+    {
+      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.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");
+    }
+
+    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 *
@@ -356,17 +385,27 @@ 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 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
+                   = std::span<const std::int32_t>())
   {
     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);
+
+    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.