Add transfer tags to submesh

cpp/dolfinx/fem/Function.h

@@ -164,7 +164,7 @@ class Function
           std::pair<std::vector<value_type>, std::vector<std::size_t>>(
               md::mdspan<const geometry_type,
                          md::extents<std::size_t, 3, md::dynamic_extent>>)>& f,
-      CellRange auto&& cells)
+      mesh::CellRange auto&& cells)
   {
     assert(_function_space);
     assert(_function_space->element());
@@ -249,7 +249,7 @@ class Function
   ///
   /// @pre `cells0` and `cells1` must have the same length.
   void interpolate(const Function<value_type, geometry_type>& u0,
-                   CellRange auto&& cells0, CellRange auto&& cells1)
+                   mesh::CellRange auto&& cells0, mesh::CellRange auto&& cells1)
   {
     fem::interpolate(*this, cells1, u0, cells0);
   }
@@ -261,7 +261,7 @@ class Function
   /// @param[in] u Function to be interpolated.
   /// @param[in] cells Cells to interpolate from.
   void interpolate(const Function<value_type, geometry_type>& u,
-                   CellRange auto&& cells)
+                   mesh::CellRange auto&& cells)
   {
     fem::interpolate(*this, u, cells);
   }
@@ -297,7 +297,7 @@ class Function
   ///
   /// @pre `cells0` `cells1` must have the same length.
   void interpolate(const Expression<value_type, geometry_type>& e0,
-                   CellRange auto&& cells0, CellRange auto&& cells1)
+                   mesh::CellRange auto&& cells0, mesh::CellRange auto&& cells1)
   {
     // Extract mesh
     const mesh::Mesh<geometry_type>* mesh0 = nullptr;
@@ -398,7 +398,7 @@ class Function
   /// interpolate from if `e0` has Function coefficients. If no mesh can
   /// be associated with `e0` then the mesh associated with `this` is used.
   void interpolate(const Expression<value_type, geometry_type>& e0,
-                   CellRange auto&& cells)
+                   mesh::CellRange auto&& cells)
   {
     interpolate(e0, cells, cells);
   }
@@ -432,7 +432,7 @@ class Function
   /// interpolation points of `this` with cells in `u`. Can be computed
   /// with `fem::create_interpolation_data`.
   void interpolate(const Function<value_type, geometry_type>& u,
-                   CellRange auto&& cells, double tol, int maxit,
+                   mesh::CellRange auto&& cells, double tol, int maxit,
                    const geometry::PointOwnershipData<U>& interpolation_data)
   {
     fem::interpolate(*this, u, cells, tol, maxit, interpolation_data);
@@ -455,7 +455,7 @@ class Function
   /// @param[in] maxit Maximum number of Newton iterations in non-affine
   /// pull-back. If the mesh geometry is affine this argument is ignored.
   void eval(std::span<const geometry_type> x, std::array<std::size_t, 2> xshape,
-            CellRange auto&& cells, std::span<value_type> u,
+            mesh::CellRange auto&& cells, std::span<value_type> u,
             std::array<std::size_t, 2> ushape, double tol, int maxit) const
   {
     if (cells.empty())

cpp/dolfinx/fem/interpolate.h

@@ -36,12 +36,6 @@
   { x.extent(1) } -> std::integral;
 };
 
-/// Requirement on range of cell indices.
-template <typename R>
-concept CellRange = std::ranges::input_range<R> and std::ranges::sized_range<R>
-                    and std::is_integral_v<
-                        std::remove_const_t<std::ranges::range_value_t<R>>>;
-
 /// @brief Compute the evaluation points in the physical space at which
 /// an expression should be computed to interpolate it in a finite
 /// element space.
@@ -55,7 +49,7 @@
 template <std::floating_point T>
 std::vector<T> interpolation_coords(const fem::FiniteElement<T>& element,
                                     const mesh::Geometry<T>& geometry,
-                                    CellRange auto&& cells)
+                                    mesh::CellRange auto&& cells)
 {
   // Find CoordinateElement appropriate to element
   auto cmap_index = [&geometry](mesh::CellType cell_type)
@@ -139,7 +133,8 @@
 /// calling `interpolation_coords`.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate(Function<T, U>& u, std::span<const T> f,
-                 std::array<std::size_t, 2> fshape, CellRange auto&& cells);
+                 std::array<std::size_t, 2> fshape,
+                 mesh::CellRange auto&& cells);
 
 namespace impl
 {
@@ -368,8 +363,9 @@
 /// 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, CellRange auto&& cells1,
-                          const Function<T, U>& u0, CellRange auto&& cells0)
+void interpolate_same_map(Function<T, U>& u1, mesh::CellRange auto&& cells1,
+                          const Function<T, U>& u0,
+                          mesh::CellRange auto&& cells0)
 {
   auto V0 = u0.function_space();
   assert(V0);
@@ -469,9 +465,10 @@
 /// @pre Function%s `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, CellRange auto&& cells1,
+void interpolate_nonmatching_maps(Function<T, U>& u1,
+                                  mesh::CellRange auto&& cells1,
                                   const Function<T, U>& u0,
-                                  CellRange auto&& cells0)
+                                  mesh::CellRange auto&& cells0)
 {
   // Get mesh
   auto V0 = u0.function_space();
@@ -719,7 +716,7 @@
 /// @param [out] coeffs Output Function coefficients.
 template <dolfinx::scalar T, std::floating_point U>
 void point_evaluation(const FiniteElement<U>& element, bool symmetric,
-                      const DofMap& dofmap, CellRange auto&& cells,
+                      const DofMap& dofmap, mesh::CellRange auto&& cells,
                       std::span<const std::uint32_t> cell_info,
                       std::span<const T> f, std::array<std::size_t, 2> fshape,
                       std::span<T> coeffs)
@@ -821,7 +818,8 @@
 /// @param [out] coeffs Output Function coefficients.
 template <dolfinx::scalar T, std::floating_point U>
 void identity_mapped_evaluation(const FiniteElement<U>& element, bool symmetric,
-                                const DofMap& dofmap, CellRange auto&& cells,
+                                const DofMap& dofmap,
+                                mesh::CellRange auto&& cells,
                                 std::span<const std::uint32_t> cell_info,
                                 std::span<const T> f,
                                 std::array<std::size_t, 2> fshape,
@@ -897,7 +895,7 @@
 /// @param [out] coeffs Output Function coefficients.
 template <dolfinx::scalar T, std::floating_point U>
 void piola_mapped_evaluation(const FiniteElement<U>& element, bool symmetric,
-                             const DofMap& dofmap, CellRange auto&& cells,
+                             const DofMap& dofmap, mesh::CellRange auto&& cells,
                              std::span<const std::uint32_t> cell_info,
                              std::span<const T> f,
                              std::array<std::size_t, 2> fshape,
@@ -1077,7 +1075,7 @@
 template <std::floating_point T>
 geometry::PointOwnershipData<T> create_interpolation_data(
     const mesh::Geometry<T>& geometry0, const FiniteElement<T>& element0,
-    const mesh::Mesh<T>& mesh1, CellRange auto&& cells, T padding)
+    const mesh::Mesh<T>& mesh1, mesh::CellRange auto&& cells, T padding)
 {
   // Collect all the points at which values are needed to define the
   // interpolating function
@@ -1097,7 +1095,8 @@
 
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate(Function<T, U>& u, std::span<const T> f,
-                 std::array<std::size_t, 2> fshape, CellRange auto&& cells)
+                 std::array<std::size_t, 2> fshape,
+                 mesh::CellRange auto&& cells)
 {
   // TODO: Index for mixed-topology, zero for now
   const int index = 0;
@@ -1180,7 +1179,7 @@
 /// fem::create_interpolation_data.
 template <dolfinx::scalar T, std::floating_point U>
 void interpolate(Function<T, U>& u1, const Function<T, U>& u0,
-                 CellRange auto&& cells, double tol, int maxit,
+                 mesh::CellRange auto&& cells, double tol, int maxit,
                  const geometry::PointOwnershipData<U>& interpolation_data)
 {
   auto mesh1 = u1.function_space()->mesh();
@@ -1257,8 +1256,8 @@
 ///
 /// @pre `cells0` and `cells1` have the same size.
 template <dolfinx::scalar T, std::floating_point U>
-void interpolate(Function<T, U>& u1, CellRange auto&& cells1,
-                 const Function<T, U>& u0, CellRange auto&& cells0)
+void interpolate(Function<T, U>& u1, mesh::CellRange auto&& cells1,
+                 const Function<T, U>& u0, mesh::CellRange auto&& cells0)
 {
   if (cells0.size() != cells1.size())
     throw std::runtime_error("Length of cell lists do not match.");

cpp/dolfinx/mesh/EntityMap.cpp

@@ -6,7 +6,6 @@
 
 #include "EntityMap.h"
 #include "Topology.h"
-#include <ranges>
 #include <span>
 #include <unordered_map>
 #include <vector>
@@ -25,55 +24,6 @@ std::shared_ptr<const Topology> EntityMap::sub_topology() const
 {
   return _sub_topology;
 }
-//-----------------------------------------------------------------------------
-std::vector<std::int32_t>
-EntityMap::sub_topology_to_topology(std::span<const std::int32_t> entities,
-                                    bool inverse) const
-{
-  if (!inverse)
-  {
-    // In this case, we want to map from entity indices in
-    // `_sub_topology` to corresponding entities in `_topology`. Hence,
-    // for each index in `entities`, we get the corresponding index in
-    // `_topology` using `_sub_topology_to_topology`
-    auto mapped
-        = entities
-          | std::views::transform([this](std::int32_t i)
-                                  { return _sub_topology_to_topology[i]; });
-    return std::vector<std::int32_t>(mapped.begin(), mapped.end());
-  }
-  else
-  {
-    // In this case, we are mapping from entity indices in `_topology`
-    // to entity indices in `_sub_topology`. Hence, we first need to
-    // construct the "inverse" of `_sub_topology_to_topology`
-    std::unordered_map<std::int32_t, std::int32_t> topology_to_sub_topology;
-    topology_to_sub_topology.reserve(_sub_topology_to_topology.size());
-    for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i)
-    {
-      topology_to_sub_topology.insert(
-          {_sub_topology_to_topology[i], static_cast<std::int32_t>(i)});
-    }
 
-    // For each entity index in `entities` (which are indices in
-    // `_topology`), get the corresponding entity in `_sub_topology`.
-    // Since `_sub_topology` consists of a subset of entities in
-    // `_topology`, there are entities in topology that may not exist in
-    // `_sub_topology`. If this is the case, mark those entities with
-    // -1.
-
-    auto mapped
-        = entities
-          | std::views::transform(
-              [&topology_to_sub_topology](std::int32_t i)
-              {
-                // Map the entity if it exists. If it doesn't, mark with
-                // -1.
-                auto it = topology_to_sub_topology.find(i);
-                return (it != topology_to_sub_topology.end()) ? it->second : -1;
-              });
-    return std::vector<std::int32_t>(mapped.begin(), mapped.end());
-  }
-}
 //-----------------------------------------------------------------------------
 } // namespace dolfinx::mesh

cpp/dolfinx/mesh/EntityMap.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2025 Jørgen S. Dokken and Joseph P. Dean
+// Copyright (C) 2025-2026 Jørgen S. Dokken and Joseph P. Dean
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -9,6 +9,7 @@
 #include "Topology.h"
 #include <concepts>
 #include <dolfinx/common/IndexMap.h>
+#include <ranges>
 #include <span>
 #include <vector>
 
@@ -99,9 +100,55 @@
   /// `this->sub_topology()`.
   /// @return A list of mapped entity indices. Entities that do not
   /// exist in the target topology are marked as -1.
-  std::vector<std::int32_t>
-  sub_topology_to_topology(std::span<const std::int32_t> entities,
-                           bool inverse) const;
+  std::vector<std::int32_t> sub_topology_to_topology(CellRange auto&& entities,
+                                                     bool inverse) const
+  {
+    if (!inverse)
+    {
+      // In this case, we want to map from entity indices in
+      // `_sub_topology` to corresponding entities in `_topology`. Hence,
+      // for each index in `entities`, we get the corresponding index in
+      // `_topology` using `_sub_topology_to_topology`
+      auto mapped
+          = std::forward<decltype(entities)>(entities)
+            | std::views::transform([this](std::int32_t i)
+                                    { return _sub_topology_to_topology[i]; });
+      return std::vector<std::int32_t>(mapped.begin(), mapped.end());
+    }
+    else
+    {
+      // In this case, we are mapping from entity indices in `_topology`
+      // to entity indices in `_sub_topology`. Hence, we first need to
+      // construct the "inverse" of `_sub_topology_to_topology`
+      std::unordered_map<std::int32_t, std::int32_t> topology_to_sub_topology;
+      topology_to_sub_topology.reserve(_sub_topology_to_topology.size());
+      for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i)
+      {
+        topology_to_sub_topology.insert(
+            {_sub_topology_to_topology[i], static_cast<std::int32_t>(i)});
+      }
+
+      // For each entity index in `entities` (which are indices in
+      // `_topology`), get the corresponding entity in `_sub_topology`.
+      // Since `_sub_topology` consists of a subset of entities in
+      // `_topology`, there are entities in topology that may not exist in
+      // `_sub_topology`. If this is the case, mark those entities with
+      // -1.
+
+      auto mapped = std::forward<decltype(entities)>(entities)
+                    | std::views::transform(
+                        [&topology_to_sub_topology](std::int32_t i)
+                        {
+                          // Map the entity if it exists. If it doesn't, mark
+                          // with -1.
+                          auto it = topology_to_sub_topology.find(i);
+                          return (it != topology_to_sub_topology.end())
+                                     ? it->second
+                                     : -1;
+                        });
+      return std::vector<std::int32_t>(mapped.begin(), mapped.end());
+    }
+  }
 
 private:
   // Dimension of the entities

cpp/dolfinx/mesh/Topology.h

@@ -7,6 +7,7 @@
 #pragma once
 
 #include <array>
+#include <concepts>
 #include <cstdint>
 #include <dolfinx/common/MPI.h>
 #include <dolfinx/graph/AdjacencyList.h>
@@ -26,6 +27,12 @@
 
 namespace dolfinx::mesh
 {
+/// Requirement on range of cell indices.
+template <typename R>
+concept CellRange = std::ranges::input_range<R> and std::ranges::sized_range<R>
+                    and std::is_integral_v<
+                        std::remove_const_t<std::ranges::range_value_t<R>>>;
+
 enum class CellType : std::int8_t;
 
 /// @brief Topology stores the topology of a mesh, consisting of mesh