Allow expressions with submeshes

cpp/demo/hyperelasticity/main.cpp

@@ -334,7 +334,7 @@ int main(int argc, char* argv[])
     // Compute Cauchy stress. Construct appropriate Basix element for
     // stress.
     fem::Expression sigma_expression = fem::create_expression<T, U>(
-        *expression_hyperelasticity_sigma, {{"u", u}}, {});
+        *expression_hyperelasticity_sigma, {{"u", u}}, {}, {});
 
     constexpr auto family = basix::element::family::P;
     auto cell_type

cpp/dolfinx/fem/Expression.h

@@ -1,4 +1,5 @@
-// Copyright (C) 2020-2025 Jack S. Hale, Michal Habera and Garth N. Wells
+// Copyright (C) 2020-2026 Jack S. Hale, Michal Habera, Garth N. Wells and
+// Jørgen S. Dokken
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -12,6 +13,7 @@
 #include <array>
 #include <concepts>
 #include <dolfinx/common/types.h>
+#include <dolfinx/mesh/EntityMap.h>
 #include <dolfinx/mesh/Mesh.h>
 #include <functional>
 #include <memory>
@@ -62,38 +64,35 @@
   /// @param[in] fn Function for tabulating the Expression.
   /// @param[in] value_shape Shape of Expression evaluated at single
   /// point.
+  /// @param[in] entity_maps Maps between entities of different meshes.
+  /// @param[in] coordinate_element_hash Hash for coordinate element used
+  /// to create the expression.
   /// @param[in] argument_space Function space for Argument. Used to
   /// computed a 1-form expression, e.g. can be used to create a matrix
   /// that when applied to a degree-of-freedom vector gives the
   /// expression values at the evaluation points.
-  Expression(const std::vector<std::shared_ptr<
-                 const Function<scalar_type, geometry_type>>>& coefficients,
-             const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
-                 constants,
-             std::span<const geometry_type> X,
-             std::array<std::size_t, 2> Xshape,
-             std::function<void(scalar_type*, const scalar_type*,
-                                const scalar_type*, const geometry_type*,
-                                const int*, const uint8_t*, void*)>
-                 fn,
-             const std::vector<std::size_t>& value_shape,
-             std::shared_ptr<const FunctionSpace<geometry_type>> argument_space
-             = nullptr)
+  Expression(
+      const std::vector<std::shared_ptr<
+          const Function<scalar_type, geometry_type>>>& coefficients,
+      const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
+          constants,
+      std::span<const geometry_type> X, std::array<std::size_t, 2> Xshape,
+      std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
+                         const geometry_type*, const int*, const uint8_t*,
+                         void*)>
+          fn,
+      const std::vector<std::size_t>& value_shape,
+      const std::vector<std::reference_wrapper<const dolfinx::mesh::EntityMap>>&
+          entity_maps,
+
+      std::uint64_t coordinate_element_hash,
+      std::shared_ptr<const FunctionSpace<geometry_type>> argument_space
+      = nullptr)
       : _argument_space(argument_space), _coefficients(coefficients),
         _constants(constants), _fn(fn), _value_shape(value_shape),
-        _x_ref(std::vector<geometry_type>(X.begin(), X.end()), Xshape)
-
-  {
-    for (auto& c : _coefficients)
-    {
-      assert(c);
-      if (c->function_space()->mesh()
-          != _coefficients.front()->function_space()->mesh())
-      {
-        throw std::runtime_error("Coefficients not all defined on same mesh.");
-      }
-    }
-  }
+        _x_ref(std::vector<geometry_type>(X.begin(), X.end()), Xshape),
+        _entity_maps(entity_maps),
+        _coordinate_element_hash(coordinate_element_hash) {};
 
   /// Move constructor
   Expression(Expression&& e) = default;
@@ -168,6 +167,19 @@
     return _x_ref;
   }
 
+  /// @brief Maps between entities of different meshes.
+  const std::vector<std::reference_wrapper<const dolfinx::mesh::EntityMap>>&
+  entity_maps() const
+  {
+    return _entity_maps;
+  }
+
+  /// @brief Hash for coordinate element used to create the expression.
+  std::uint64_t coordinate_element_hash() const
+  {
+    return _coordinate_element_hash;
+  }
+
 private:
   // Function space for Argument
   std::shared_ptr<const FunctionSpace<geometry_type>> _argument_space;
@@ -189,5 +201,12 @@
 
   // Evaluation points on reference cell
   std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _x_ref;
+
+  // Map between different mesh topologies
+  std::vector<std::reference_wrapper<const dolfinx::mesh::EntityMap>>
+      _entity_maps;
+
+  // Hash for coordinate element used to create the expression
+  std::uint64_t _coordinate_element_hash;
 };
 } // namespace dolfinx::fem

cpp/dolfinx/fem/assembler.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2018-2025 Garth N. Wells
+// Copyright (C) 2018-2026 Garth N. Wells and Jørgen S. Dokken
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -19,6 +19,7 @@
 #include <basix/mdspan.hpp>
 #include <cstdint>
 #include <dolfinx/common/types.h>
+#include <dolfinx/mesh/EntityMap.h>
 #include <memory>
 #include <optional>
 #include <span>
@@ -70,6 +71,12 @@
         std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
         element)
 {
+  // Check that domain is the same as mesh of the expression
+  if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash())
+  {
+    throw std::runtime_error(
+        "Expression was created on a different mesh. Cannot tabulate.");
+  }
   auto [X, Xshape] = e.X();
   impl::tabulate_expression(values, e.kernel(), Xshape, e.value_size(), coeffs,
                             constants, mesh, entities, element);
@@ -95,6 +102,13 @@
 void tabulate_expression(std::span<T> values, const fem::Expression<T, U>& e,
                          const mesh::Mesh<U>& mesh, fem::MDSpan2 auto entities)
 {
+  // Check that domain is the same as mesh of the expression
+  if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash())
+  {
+    throw std::runtime_error(
+        "Expression was created on a different mesh. Cannot tabulate.");
+  }
+
   std::optional<
       std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
       element = std::nullopt;
@@ -115,7 +129,8 @@
     std::vector<std::reference_wrapper<const Function<T, U>>> c;
     std::ranges::transform(coefficients, std::back_inserter(c),
                            [](auto c) -> const Function<T, U>& { return *c; });
-    fem::pack_coefficients(c, coffsets, entities, std::span(coeffs));
+    fem::pack_coefficients(c, mesh, entities, e.entity_maps(), coffsets,
+                           std::span(coeffs));
   }
   std::vector<T> constants = fem::pack_constants(e);
 

cpp/dolfinx/fem/pack.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2013-2025 Garth N. Wells
+// Copyright (C) 2013-2026 Garth N. Wells and Jørgen S. Dokken
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -17,6 +17,7 @@
 #include <basix/mdspan.hpp>
 #include <concepts>
 #include <dolfinx/mesh/Topology.h>
+#include <ranges>
 #include <span>
 #include <stdexcept>
 #include <type_traits>
@@ -332,17 +333,23 @@
 /// @brief Pack coefficient data over a list of cells or facets.
 ///
 /// Typically used to prepare coefficient data for an ::Expression.
-/// @tparam T
-/// @tparam U
+/// @tparam T data type of coefficients
+/// @tparam U floating point type of mesh geometry
 /// @param coeffs Coefficients to pack
-/// @param offsets Offsets
+/// @param mesh Mesh which the entities belong to
 /// @param entities Entities to pack over
-/// @param[out] c Packed coefficients.
+/// @param entity_maps Entity maps for the coefficient meshes
+/// @param offsets Offsets
+/// @param[in,out] c Packed coefficients.
 template <dolfinx::scalar T, std::floating_point U>
 void pack_coefficients(
     std::vector<std::reference_wrapper<const Function<T, U>>> coeffs,
-    std::span<const int> offsets, fem::MDSpan2 auto entities, std::span<T> c)
+    const mesh::Mesh<U>& mesh, fem::MDSpan2 auto entities,
+    const std::vector<std::reference_wrapper<const dolfinx::mesh::EntityMap>>&
+        entity_maps,
+    std::span<const int> offsets, std::span<T> c)
 {
+
   assert(!offsets.empty());
   const int cstride = offsets.back();
 
@@ -354,17 +361,105 @@
   {
     std::span<const std::uint32_t> cell_info
         = impl::get_cell_orientation_info(coeffs[coeff].get());
-
-    if constexpr (entities.rank() == 1)
+    // Get mesh of coefficient and check if entity map is required
+    auto mesh_c = coeffs[coeff].get().function_space()->mesh();
+    if (mesh_c->topology() == mesh.topology())
     {
-      impl::pack_coefficient_entity(std::span(c), cstride, coeffs[coeff].get(),
-                                    cell_info, entities, offsets[coeff]);
+      // If same mesh no mapping is needed
+      if constexpr (entities.rank() == 1)
+      {
+        impl::pack_coefficient_entity(std::span(c), cstride,
+                                      coeffs[coeff].get(), cell_info, entities,
+                                      offsets[coeff]);
+      }
+      else
+      {
+        auto cells = md::submdspan(entities, md::full_extent, 0);
+        impl::pack_coefficient_entity(std::span(c), cstride,
+                                      coeffs[coeff].get(), cell_info, cells,
+                                      offsets[coeff]);
+      }
     }
     else
     {
-      auto cells = md::submdspan(entities, md::full_extent, 0);
-      impl::pack_coefficient_entity(std::span(c), cstride, coeffs[coeff].get(),
-                                    cell_info, cells, offsets[coeff]);
+      // Helper function to get correct entity map
+      auto get_entity_map
+          = [mesh, &entity_maps](auto& mesh0) -> const mesh::EntityMap&
+      {
+        auto it = std::ranges::find_if(
+            entity_maps,
+            [mesh, mesh0](const mesh::EntityMap& em)
+            {
+              return ((em.topology() == mesh0->topology()
+                       and em.sub_topology() == mesh.topology()))
+                     or ((em.sub_topology() == mesh0->topology()
+                          and em.topology() == mesh.topology()));
+            });
+
+        if (it == entity_maps.end())
+        {
+          throw std::runtime_error(
+              "Incompatible mesh. argument entity_maps must be provided.");
+        }
+        return *it;
+      };
+      // Find correct entity map and determine direction of the map
+      const mesh::EntityMap& emap = get_entity_map(mesh_c);
+      bool inverse = emap.sub_topology() == mesh_c->topology();
+
+      std::vector<std::int32_t> e_b;
+      if constexpr (entities.rank() == 1)
+      {
+        e_b = emap.sub_topology_to_topology(
+            std::span(entities.data_handle(), entities.size()), inverse);
+        md::mdspan e(e_b.data(), e_b.size());
+        impl::pack_coefficient_entity(std::span(c), cstride,
+                                      coeffs[coeff].get(), cell_info, e,
+                                      offsets[coeff]);
+      }
+      else if (entities.rank() == 2)
+      {
+        const mesh::Topology topology = *mesh.topology();
+        int tdim = topology.dim();
+        int codim = tdim - mesh_c->topology()->dim();
+        if (codim == 0)
+        {
+          // If codim is zero we extract the cells and map them
+          auto cells = md::submdspan(entities, md::full_extent, 0);
+          std::vector<std::int32_t> contiguous_cells(cells.extent(0));
+          for (std::size_t i = 0; i < cells.extent(0); ++i)
+            contiguous_cells[i] = cells(i);
+          e_b = emap.sub_topology_to_topology(std::span(contiguous_cells),
+                                              inverse);
+          md::mdspan e(e_b.data(), e_b.size());
+          impl::pack_coefficient_entity(std::span(c), cstride,
+                                        coeffs[coeff].get(), cell_info, e,
+                                        offsets[coeff]);
+        }
+        else if (codim == 1)
+        {
+          // Codim 1 mesh, need to map (cell, local index) to facets and then
+          // to cells of the submesh
+          auto c_to_f = topology.connectivity(tdim, tdim - 1);
+          assert(c_to_f);
+          std::vector<std::int32_t> parent_entities;
+          parent_entities.reserve(entities.extent(0));
+          for (std::size_t e = 0; e < entities.extent(0); ++e)
+          {
+            auto pair = md::submdspan(entities, e, md::full_extent);
+            parent_entities.push_back(c_to_f->links(pair[0])[pair[1]]);
+          }
+          e_b = emap.sub_topology_to_topology(parent_entities, inverse);
+          md::mdspan e(e_b.data(), e_b.size());
+          impl::pack_coefficient_entity(std::span(c), cstride,
+                                        coeffs[coeff].get(), cell_info, e,
+                                        offsets[coeff]);
+        }
+      }
+      else
+      {
+        throw std::runtime_error("Entities span has unsupported rank.");
+      }
     }
   }
 }

cpp/dolfinx/fem/utils.h

@@ -953,21 +953,10 @@ Expression<T, U> create_expression(
     const ufcx_expression& e,
     const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
     const std::vector<std::shared_ptr<const Constant<T>>>& constants,
+    const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
+        entity_maps,
     std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
 {
-  if (!coefficients.empty())
-  {
-    assert(coefficients.front());
-    assert(coefficients.front()->function_space());
-    std::shared_ptr<const mesh::Mesh<U>> mesh
-        = coefficients.front()->function_space()->mesh();
-    if (mesh->geometry().cmap().hash() != e.coordinate_element_hash)
-    {
-      throw std::runtime_error(
-          "Expression and mesh geometric maps do not match.");
-    }
-  }
-
   if (e.rank > 0 and !argument_space)
   {
     throw std::runtime_error("Expression has Argument but no Argument "
@@ -1007,8 +996,10 @@ Expression<T, U> create_expression(
     throw std::runtime_error("Type not supported.");
 
   assert(tabulate_tensor);
+  std::uint64_t e_hash = e.coordinate_element_hash;
   return Expression(coefficients, constants, std::span<const U>(X), Xshape,
-                    tabulate_tensor, value_shape, argument_space);
+                    tabulate_tensor, value_shape, entity_maps, e_hash,
+                    argument_space);
 }
 
 /// @brief Create Expression from UFC input (with named coefficients and
@@ -1019,6 +1010,8 @@ Expression<T, U> create_expression(
     const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
         coefficients,
     const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
+    const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
+        entity_maps,
     std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
 {
   // Place coefficients in appropriate order
@@ -1057,7 +1050,8 @@ Expression<T, U> create_expression(
     }
   }
 
-  return create_expression(e, coeff_map, const_map, argument_space);
+  return create_expression(e, coeff_map, const_map, entity_maps,
+                           argument_space);
 }
 
 } // namespace dolfinx::fem