diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 35349f1cfbc..e0d40e051b0 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/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
diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h
index eb524f0a234..c637f9c0869 100644
--- a/cpp/dolfinx/fem/Expression.h
+++ b/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 @@ class Expression
   /// @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 @@ class Expression
     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 @@ class Expression
 
   // 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
diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
index 352196d502f..732fd988fe5 100644
--- a/cpp/dolfinx/fem/assembler.h
+++ b/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 @@ void tabulate_expression(
         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 @@ template <dolfinx::scalar T, std::floating_point U>
 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 @@ void tabulate_expression(std::span<T> values, const fem::Expression<T, U>& e,
     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);
 
diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h
index 65777820ee8..3706f3fb2c8 100644
--- a/cpp/dolfinx/fem/pack.h
+++ b/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 @@ void pack_coefficients(const Form<T, U>& form,
 /// @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 @@ void pack_coefficients(
   {
     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.");
+      }
     }
   }
 }
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 81fbc2f953d..763ea8cd8ea 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/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
diff --git a/cpp/test/fem/form.cpp b/cpp/test/fem/form.cpp
index 348ae15c5eb..fe52fa1128e 100644
--- a/cpp/test/fem/form.cpp
+++ b/cpp/test/fem/form.cpp
@@ -35,7 +35,7 @@ void test_expression_cmap_compat(const auto& V)
   // Create Expression that expects P1 geometry
   dolfinx::fem::Expression<double> expr1
       = dolfinx::fem::create_expression<double>(*expression_expr_Q6_P1,
-                                                {{"u1", u}}, {});
+                                                {{"u1", u}}, {}, {});
   auto [Xc, Xshape] = expr1.X();
   std::vector<double> grad_e(3 * Xshape[0] * cells.size());
   fem::tabulate_expression(std::span(grad_e), expr1, *mesh,
@@ -43,8 +43,11 @@ void test_expression_cmap_compat(const auto& V)
 
   // Create Expression that expects P2 geometry. Should throw because
   // mesh is P1.
-  CHECK_THROWS(dolfinx::fem::create_expression<double>(*expression_expr_Q6_P2,
-                                                       {{"u2", u}}, {}));
+  dolfinx::fem::Expression<double> expr2
+      = dolfinx::fem::create_expression<double>(*expression_expr_Q6_P2,
+                                                {{"u2", u}}, {}, {});
+  CHECK_THROWS(fem::tabulate_expression(
+      std::span(grad_e), expr2, *mesh, md::mdspan(cells.data(), cells.size())));
 }
 } // namespace
 
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 69c9cce7f2e..c64cf90d8c9 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -26,6 +26,7 @@
 if typing.TYPE_CHECKING:
     from mpi4py import MPI as _MPI
 
+    from dolfinx.mesh import EntityMap as _EntityMap
     from dolfinx.mesh import Mesh
 
 
@@ -117,6 +118,7 @@ def __init__(
         form_compiler_options: dict | None = None,
         jit_options: dict | None = None,
         dtype: npt.DTypeLike | None = None,
+        entity_maps: list[_EntityMap] | None = None,
     ):
         """Create an Expression.
 
@@ -131,6 +133,7 @@ def __init__(
             jit_options: Options controlling JIT compilation of C code.
             dtype: Type of the Expression values. If ``None``,
                 the dtype is deduced from the UFL expression.
+            entity_maps: Maps between different meshes.
 
         Note:
             This wrapper is responsible for the FFCx compilation of the
@@ -144,7 +147,10 @@ def __init__(
         # Get MPI communicator
         if comm is None:
             try:
-                domain = ufl.domain.extract_unique_domain(e)
+                domains = ufl.domain.extract_domains(e)
+                if len(domains) == 0:
+                    raise RuntimeError("Could not extract MPI communicator for Expression.")
+                domain = domains[0]
                 assert isinstance(domain, ufl.Mesh)
                 mesh = domain.ufl_cargo()
                 comm = mesh.comm
@@ -199,11 +205,17 @@ def _create_expression(dtype):
             else:
                 raise NotImplementedError(f"Type {dtype} not supported.")
 
+        _entity_maps = (
+            [entity_map._cpp_object for entity_map in entity_maps]
+            if entity_maps is not None
+            else []
+        )
         ffi = module.ffi
         self._cpp_object = _create_expression(dtype)(
             ffi.cast("uintptr_t", ffi.addressof(self._ufcx_expression)),
             coeffs,
             constants,
+            _entity_maps,
             self.argument_space,
         )
 
@@ -262,7 +274,7 @@ def eval(
                 raise TypeError("Passed values array does not have correct dtype.")
 
         constants = _cpp.fem.pack_constants(self._cpp_object)
-        coeffs = _cpp.fem.pack_coefficients(self._cpp_object, _entities)
+        coeffs = _cpp.fem.pack_coefficients(self._cpp_object, mesh._cpp_object, _entities)
         _cpp.fem.tabulate_expression(
             values, self._cpp_object, constants, coeffs, mesh._cpp_object, _entities
         )
diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h b/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h
index dce4671a76e..082eef4c068 100644
--- a/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h
+++ b/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h
@@ -212,6 +212,7 @@ void declare_assembly_functions(nanobind::module_& m)
   m.def(
       "pack_coefficients",
       [](const dolfinx::fem::Expression<T, U>& e,
+         const dolfinx::mesh::Mesh<U>& mesh,
          nb::ndarray<const std::int32_t, nb::c_contig> entities)
       {
         std::vector<int> coffsets = e.coefficient_offsets();
@@ -228,22 +229,22 @@ void declare_assembly_functions(nanobind::module_& m)
         if (entities.ndim() == 1)
         {
           dolfinx::fem::pack_coefficients(
-              c, coffsets, md::mdspan(entities.data(), entities.shape(0)),
-              std::span(coeffs));
+              c, mesh, md::mdspan(entities.data(), entities.shape(0)),
+              e.entity_maps(), coffsets, std::span(coeffs));
         }
         else
         {
           dolfinx::fem::pack_coefficients(
-              c, coffsets,
+              c, mesh,
               md::mdspan<const std::int32_t,
                          md::extents<std::size_t, md::dynamic_extent, 2>>(
                   entities.data(), entities.shape(0), entities.shape(1)),
-              std::span(coeffs));
+              e.entity_maps(), coffsets, std::span(coeffs));
         }
         return dolfinx_wrappers::as_nbarray(std::move(coeffs),
                                             {entities.shape(0), cstride});
       },
-      nb::arg("expression"), nb::arg("entities"),
+      nb::arg("expression"), nb::arg("mesh"), nb::arg("entities"),
       "Pack coefficients for a Expression.");
   m.def(
       "pack_constants",
diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h
index d208ff30b1b..84799555cbc 100644
--- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h
+++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h
@@ -683,6 +683,8 @@ void declare_objects(nb::module_& m, std::string type)
              nb::ndarray<const U, nb::ndim<2>, nb::c_contig> X,
              std::uintptr_t fn_addr,
              const std::vector<std::size_t>& value_shape,
+             const std::vector<const dolfinx::mesh::EntityMap*>& entity_maps,
+             std::uint64_t coordinate_element_hash,
              std::shared_ptr<const dolfinx::fem::FunctionSpace<U>>
                  argument_space)
           {
@@ -693,10 +695,12 @@ void declare_objects(nb::module_& m, std::string type)
             new (ex) dolfinx::fem::Expression<T, U>(
                 coefficients, constants, std::span(X.data(), X.size()),
                 {X.shape(0), X.shape(1)}, tabulate_expression_ptr, value_shape,
+                ptr_to_ref_wrapper_vec(entity_maps), coordinate_element_hash,
                 argument_space);
           },
           nb::arg("coefficients"), nb::arg("constants"), nb::arg("X"),
-          nb::arg("fn"), nb::arg("value_shape"), nb::arg("argument_space"))
+          nb::arg("fn"), nb::arg("value_shape"), nb::arg("entity_maps"),
+          nb::arg("coordinate_element_hash"), nb::arg("argument_space"))
       .def("X",
            [](const dolfinx::fem::Expression<T, U>& self)
            {
@@ -725,15 +729,17 @@ void declare_objects(nb::module_& m, std::string type)
              coefficients,
          const std::vector<std::shared_ptr<const dolfinx::fem::Constant<T>>>&
              constants,
+         const std::vector<const dolfinx::mesh::EntityMap*>& entity_maps,
          std::shared_ptr<const dolfinx::fem::FunctionSpace<U>> argument_space)
       {
         const ufcx_expression* p
             = reinterpret_cast<const ufcx_expression*>(expression);
-        return dolfinx::fem::create_expression<T, U>(*p, coefficients,
-                                                     constants, argument_space);
+        return dolfinx::fem::create_expression<T, U>(
+            *p, coefficients, constants, ptr_to_ref_wrapper_vec(entity_maps),
+            argument_space);
       },
       nb::arg("expression"), nb::arg("coefficients"), nb::arg("constants"),
-      nb::arg("argument_space").none(),
+      nb::arg("entity_maps"), nb::arg("argument_space").none(),
       "Create Expression from a pointer to ufc_form.");
 }
 
diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py
index c781eb4f261..c5804480150 100644
--- a/python/test/unit/fem/test_expression.py
+++ b/python/test/unit/fem/test_expression.py
@@ -1,4 +1,4 @@
-# Copyright (C) 2019-2024 Michal Habera and Jørgen S. Dokken
+# Copyright (C) 2019-2026 Michal Habera and Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -15,7 +15,7 @@
 from basix.ufl import quadrature_element
 from dolfinx import fem, la
 from dolfinx.fem import Constant, Expression, Function, form, functionspace
-from dolfinx.mesh import create_unit_square
+from dolfinx.mesh import create_rectangle, create_unit_square
 
 
 @pytest.mark.parametrize(
@@ -526,3 +526,112 @@ def test_rank1_blocked():
                 mask = np.ones(point_values.shape[2], dtype=bool)
                 mask[offset::vs] = False
                 np.testing.assert_allclose(point_values[i, j, mask], 0)
+
+
+@pytest.mark.parametrize("qdegree", [1, 3, 5])
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(np.complex64, marks=pytest.mark.xfail_win32_complex),
+        pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex),
+    ],
+)
+def test_submesh_codim_zero(dtype, qdegree):
+    xtype = dtype(0).real.dtype
+    mesh = create_unit_square(MPI.COMM_WORLD, 4, 3, dtype=xtype)
+
+    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
+    u = dolfinx.fem.Function(V, dtype=dtype)
+    u.interpolate(lambda x: x[0] + 2.0 * x[1])
+
+    def mark_left_cells(x):
+        return x[0] <= 0.5 + 100 * np.finfo(xtype).resolution
+
+    left_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, mark_left_cells)
+    submesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, left_cells)
+    u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, ("Lagrange", 2)), dtype=dtype)
+    u_sub.interpolate(lambda x: x[1] ** 2 + x[0] ** 2)
+
+    quadrature_points, _ = basix.make_quadrature(basix.CellType.triangle, qdegree)
+    quadrature_points = quadrature_points.astype(xtype)
+
+    sub_cellmap = submesh.topology.index_map(submesh.topology.dim)
+    num_sub_cells = sub_cellmap.size_local + sub_cellmap.num_ghosts
+    parent_cells = entity_map.sub_topology_to_topology(
+        np.arange(num_sub_cells, dtype=np.int32), inverse=False
+    )
+
+    expr = dolfinx.fem.Expression(
+        u * u_sub, quadrature_points, dtype=dtype, entity_maps=[entity_map]
+    )
+    values = expr.eval(mesh, parent_cells)
+
+    values_sub = expr.eval(submesh, np.arange(num_sub_cells, dtype=np.int32))
+
+    tol = 50 * np.finfo(dtype).eps
+    np.testing.assert_allclose(values, values_sub, atol=tol)
+
+    x = ufl.SpatialCoordinate(mesh)
+    u_exact = (x[0] + 2.0 * x[1]) * (x[1] ** 2 + x[0] ** 2)
+    cell_map = mesh.topology.index_map(mesh.topology.dim)
+    num_cells = cell_map.size_local + cell_map.num_ghosts
+    expr_exact = dolfinx.fem.Expression(u_exact, quadrature_points, dtype=dtype)
+    values_exact = expr_exact.eval(mesh, np.arange(num_cells, dtype=np.int32))
+    np.testing.assert_allclose(values, values_exact[parent_cells], atol=tol)
+
+
+@pytest.mark.parametrize("qdegree", [1, 3, 5])
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        np.float32,
+        np.float64,
+        pytest.param(np.complex64, marks=pytest.mark.xfail_win32_complex),
+        pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex),
+    ],
+)
+def test_submesh_codim_one(dtype, qdegree):
+    xtype = dtype(0).real.dtype
+    mesh = create_rectangle(
+        MPI.COMM_WORLD, np.array([[0, 0], [2.1, 2.0]], dtype=xtype), [5, 3], dtype=xtype
+    )
+    el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(), dtype=xtype)
+    V = dolfinx.fem.functionspace(mesh, el)
+    u = dolfinx.fem.Function(V, dtype=dtype)
+    u.interpolate(lambda x: x[0] + 2.0 * x[1])
+
+    tol = 50 * np.finfo(xtype).resolution
+
+    def mark_left_facets(x):
+        return np.isclose(x[0], 1.0, atol=tol)
+
+    left_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim - 1, mark_left_facets)
+    submesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, left_facets)
+    sub_el = basix.ufl.element(
+        "Lagrange", submesh.basix_cell(), 2, shape=(submesh.geometry.dim,), dtype=xtype
+    )
+    u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, sub_el), dtype=dtype)
+    u_sub.interpolate(lambda x: (x[1] ** 2, -(x[0] ** 2)))
+
+    quadrature_points, _ = basix.make_quadrature(basix.CellType.interval, qdegree)
+    quadrature_points = quadrature_points.astype(xtype)
+
+    n_h = ufl.FacetNormal(mesh)
+    expr = u * (ufl.dot(u_sub, n_h) + n_h[0] * u_sub[1])
+
+    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
+    expr = dolfinx.fem.Expression(expr, quadrature_points, dtype=dtype, entity_maps=[entity_map])
+
+    entities = dolfinx.fem.compute_integration_domains(
+        dolfinx.fem.IntegralType.exterior_facet, mesh.topology, left_facets
+    )
+
+    values = expr.eval(mesh, entities.reshape(-1, 2))
+
+    x = ufl.SpatialCoordinate(mesh)
+    expr_exact = (x[0] + 2.0 * x[1]) * (x[1] ** 2 - x[0] ** 2)
+    expr_exact = dolfinx.fem.Expression(expr_exact, quadrature_points, dtype=dtype)
+    values_exact = expr_exact.eval(mesh, entities.reshape(-1, 2))
+    np.testing.assert_allclose(values, values_exact, atol=tol)