Correct and consolidate Form.h's entity-domain mapping

.github/workflows/fenicsx-refs.env

@@ -3,4 +3,4 @@ basix_ref=main
 ufl_repository=FEniCS/ufl
 ufl_ref=main
 ffcx_repository=FEniCS/ffcx
-ffcx_ref=main
+ffcx_ref=dokken/manifold-fixes

cpp/dolfinx/fem/Form.h

@@ -18,6 +18,7 @@
 #include <dolfinx/mesh/EntityMap.h>
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/cell_types.h>
+#include <format>
 #include <functional>
 #include <map>
 #include <memory>
@@ -47,6 +48,30 @@
   ridge = 4           ///< Ridge
 };
 
+/// @brief Topological dimension of the mesh entities an integral of the
+/// given type is over.
+///
+/// @param[in] type Integral type.
+/// @param[in] tdim Topological dimension of the integration domain.
+/// @return Entity dimension, equal to `tdim` for a cell integral.
+constexpr int integral_entity_dim(IntegralType type, int tdim)
+{
+  switch (type)
+  {
+  case IntegralType::exterior_facet:
+  case IntegralType::interior_facet:
+    return tdim - 1;
+  case IntegralType::ridge:
+    return tdim - 2;
+  case IntegralType::vertex:
+    return 0;
+  case IntegralType::cell:
+    return tdim;
+  }
+
+  throw std::invalid_argument("Unknown integral type.");
+}
+
 namespace impl
 {
 /// @brief Permutations of the cell-local entities that an integral of
@@ -66,23 +91,11 @@
 entity_permutations(mesh::Topology& topology, IntegralType type,
                     mesh::CellType cell_type)
 {
-  const int tdim = topology.dim();
-  int dim = tdim;
-  switch (type)
-  {
-  case IntegralType::exterior_facet:
-  case IntegralType::interior_facet:
-    dim = tdim - 1;
-    break;
-  case IntegralType::ridge:
-    dim = tdim - 2;
-    break;
-  case IntegralType::vertex:
-    dim = 0;
-    break;
-  case IntegralType::cell:
+  if (type == IntegralType::cell)
     return {};
-  }
+
+  const int tdim = topology.dim();
+  const int dim = integral_entity_dim(type, tdim);
 
   topology.create_entity_permutations(dim);
   const std::vector<std::uint8_t>& p = topology.get_entity_permutations(dim);
@@ -245,11 +258,11 @@
       return *it;
     };
 
-    // A helper function to compute the (cell, local_facet) pairs in the
-    // argument/coefficient domain from the (cell, local_facet) pairs in
+    // A helper function to compute the (cell, local_entity) pairs in the
+    // argument/coefficient domain from the (cell, local_entity) pairs in
     // `this->mesh()`.
-    auto compute_facet_domains
-        = [](const auto& int_ents_mesh, int codim, const auto& c_to_f,
+    auto compute_entity_domains
+        = [](const auto& int_ents_mesh, int codim, const auto& c_to_e,
              const auto& emap, bool inverse)
     {
       // TODO: This function would be much neater using
@@ -267,36 +280,77 @@
         for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
           entities.push_back(int_ents_mesh[i]);
       }
-      else if (codim == 1)
+      else
       {
-        // In the codim 1 case, we need to map facets in `this->mesh()`
-        // to cells in the argument/coefficient mesh, so here we extract
-        // the facet index using the cell-to-facet connectivity.
+        // Otherwise the integration entities are sub-entities of the
+        // cells of `this->mesh()` and are themselves cells of the
+        // argument/coefficient mesh, so here we extract the entity
+        // index using the cell-to-entity connectivity.
         for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
         {
           entities.push_back(
-              c_to_f->links(int_ents_mesh[i])[int_ents_mesh[i + 1]]);
+              c_to_e->links(int_ents_mesh[i])[int_ents_mesh[i + 1]]);
         }
       }
-      else
-        throw std::invalid_argument("Codimension > 1 not supported.");
 
       // Map from entity indices in `this->mesh()` to the corresponding
       // cell indices in the argument/coefficient mesh
       std::vector<std::int32_t> cells_mesh0
           = emap.sub_topology_to_topology(entities, inverse);
 
-      // Create a list of (cell, local_facet_index) pairs in the
-      // argument/coefficient domain. Since `create_submesh`preserves
-      // the local facet index (with respect to the cell), we can use
-      // the local facet indices from the input integration entities
+      // Create a list of (cell, local_entity_index) pairs in the
+      // argument/coefficient domain. Only the cell column is meaningful.
+      // For codim > 0 the entity is itself the cell, so it has
+      // no local index. The colummn is never used later on, but written
+      // like this for consistency for packing/assembly.
       std::vector<std::int32_t> e = int_ents_mesh;
       for (std::size_t i = 0; i < cells_mesh0.size(); ++i)
         e[2 * i] = cells_mesh0[i];
 
       return e;
     };
 
+    // Map the integration entities of one integral to the
+    // argument/coefficient domain, checking first that the mapping is
+    // expressible: `compute_entity_domains` maps an integration entity
+    // to a *cell* of that mesh, so unless the meshes have equal
+    // dimension the integral's entity dimension must be that mesh's.
+    auto map_entities = [tdim, &topology, &compute_entity_domains](
+                            IntegralType type, const auto& entities,
+                            const mesh::Mesh<geometry_type>& mesh0,
+                            const mesh::EntityMap& emap,
+                            bool inverse) -> std::vector<std::int32_t>
+    {
+      if (type == IntegralType::cell)
+        return emap.sub_topology_to_topology(entities, inverse);
+
+      if (type == IntegralType::vertex)
+      {
+        throw std::invalid_argument(
+            "Vertex integrals are not supported for a form with an argument "
+            "or coefficient on another mesh. Supported types are cell, "
+            "exterior facet, interior facet and ridge.");
+      }
+
+      const int dim0 = mesh0.topology()->dim();
+      const int codim = tdim - dim0;
+      const int edim = integral_entity_dim(type, tdim);
+      assert(codim >= 0);
+      if (codim > 0 and edim != dim0)
+      {
+        throw std::invalid_argument(std::format(
+            "Cannot map integration entities of dimension {} to cells of a "
+            "mesh of dimension {}. An argument or coefficient on another mesh "
+            "must live on the entities being integrated over.",
+            edim, dim0));
+      }
+
+      std::shared_ptr<const graph::AdjacencyList<std::int32_t>> c_to_e
+          = topology.connectivity(tdim, edim);
+      assert(codim == 0 or c_to_e);
+      return compute_entity_domains(entities, codim, c_to_e, emap, inverse);
+    };
+
     _edata.reserve(_function_spaces.size());
     for (auto& space : _function_spaces)
     {
@@ -324,21 +378,8 @@
         {
           auto [type, idx, kernel_idx] = key;
           std::vector<std::int32_t> e;
-          if (type == IntegralType::cell)
-            e = emap.sub_topology_to_topology(itg.entities, inverse);
-          else if (type == IntegralType::exterior_facet
-                   or type == IntegralType::interior_facet)
-          {
-            assert(mesh0);
-            int codim = tdim - mesh0->topology()->dim();
-            assert(codim >= 0);
-            auto c_to_f = topology.connectivity(tdim, tdim - 1);
-            assert(c_to_f);
-            e = compute_facet_domains(itg.entities, codim, c_to_f, emap,
-                                      inverse);
-          }
-          else
-            throw std::invalid_argument("Integral type not supported.");
+          assert(mesh0);
+          e = map_entities(type, itg.entities, *mesh0, emap, inverse);
 
           vdata.insert({key, std::move(e)});
         }
@@ -364,20 +405,8 @@
           bool inverse = emap.sub_topology() == mesh0->topology();
 
           std::vector<std::int32_t> e;
-          if (type == IntegralType::cell)
-            e = emap.sub_topology_to_topology(integral.entities, inverse);
-          else if (type == IntegralType::exterior_facet
-                   or type == IntegralType::interior_facet)
-          {
-            assert(mesh0);
-            int codim = tdim - mesh0->topology()->dim();
-            auto c_to_f = topology.connectivity(tdim, tdim - 1);
-            assert(c_to_f);
-            e = compute_facet_domains(integral.entities, codim, c_to_f, emap,
-                                      inverse);
-          }
-          else
-            throw std::invalid_argument("Integral type not supported.");
+          assert(mesh0);
+          e = map_entities(type, integral.entities, *mesh0, emap, inverse);
           _cdata.insert({{type, idx, c}, std::move(e)});
         }
       }

cpp/dolfinx/fem/utils.cpp

@@ -157,28 +157,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type,
 {
   const int tdim = topology.dim();
 
-  int dim = -1;
-  switch (integral_type)
-  {
-  case IntegralType::cell:
-    dim = tdim;
-    break;
-  case IntegralType::exterior_facet:
-    dim = tdim - 1;
-    break;
-  case IntegralType::interior_facet:
-    dim = tdim - 1;
-    break;
-  case IntegralType::vertex:
-    dim = 0;
-    break;
-  case IntegralType::ridge:
-    dim = tdim - 2;
-    break;
-  default:
-    throw std::invalid_argument(
-        "Cannot compute integration domains. Integral type not supported.");
-  }
+  const int dim = integral_entity_dim(integral_type, tdim);
 
   {
     // Create span of the owned entities (leaves off any ghosts)

cpp/dolfinx/fem/utils.h

@@ -707,27 +707,7 @@ Form<T, U> create_form_factory(
     for (IntegralType itg_type : {IntegralType::exterior_facet,
                                   IntegralType::vertex, IntegralType::ridge})
     {
-      std::size_t dim;
-      switch (itg_type)
-      {
-      case IntegralType::exterior_facet:
-      {
-        dim = tdim - 1;
-        break;
-      }
-      case IntegralType::ridge:
-      {
-        dim = tdim - 2;
-        break;
-      }
-      case IntegralType::vertex:
-      {
-        dim = 0;
-        break;
-      }
-      default:
-        throw std::invalid_argument("Unsupported integral type");
-      }
+      const std::size_t dim = integral_entity_dim(itg_type, tdim);
 
       const std::function<std::vector<std::int32_t>(const mesh::Topology&,
                                                     IntegralType)>

cpp/test/fem/form.cpp

@@ -1,4 +1,4 @@
-// Copyright (C) 2025 Garth N. Wells
+// Copyright (C) 2025-2026 Garth N. Wells and Jørgen S. Dokken
 //
 // This file is part of DOLFINX (https://www.fenicsproject.org)
 //
@@ -7,13 +7,15 @@
 #include "expr.h"
 #include <basix/finite-element.h>
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
 #include <dolfinx/fem/Expression.h>
 #include <dolfinx/fem/Function.h>
 #include <dolfinx/fem/FunctionSpace.h>
 #include <dolfinx/fem/assembler.h>
 #include <dolfinx/fem/utils.h>
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/generation.h>
+#include <dolfinx/mesh/utils.h>
 
 using namespace dolfinx;
 
@@ -70,3 +72,52 @@ TEST_CASE("Create Expression/Form (mismatch of mesh geometry)",
   test_form_cmap_compat(V);
   test_expression_cmap_compat(V);
 }
+
+TEST_CASE("Form with data on a mesh of the wrong dimension",
+          "[form_entity_domain]")
+{
+  // `Form` maps an integration entity to a *cell* of the argument's
+  // mesh, so a ridge integral cannot take data on a facet submesh. FFCx
+  // rejects the combination when compiling a kernel, so this is reached
+  // only by building a Form directly, as here.
+  auto mesh = std::make_shared<mesh::Mesh<double>>(mesh::create_box<double>(
+      MPI_COMM_WORLD, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {2, 2, 2},
+      mesh::CellType::tetrahedron, graph::partition_graph));
+  const int tdim = mesh->topology()->dim();
+  mesh->topology_mutable()->create_entities(tdim - 1);
+  mesh->topology_mutable()->create_entities(tdim - 2);
+  mesh->topology_mutable()->create_connectivity(tdim, tdim - 2);
+
+  // A submesh of the facets, i.e. one dimension too high for a ridge
+  auto facets = mesh::locate_entities(
+      *mesh, tdim - 1,
+      [](auto x) { return std::vector<std::int8_t>(x.extent(1), 1); });
+  auto [submesh, e_map, v_map, g_map]
+      = mesh::create_submesh(*mesh, tdim - 1, facets);
+  auto smesh = std::make_shared<mesh::Mesh<double>>(std::move(submesh));
+
+  auto element = basix::create_element<double>(
+      basix::element::family::P, basix::cell::type::triangle, 1,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+  auto V = std::make_shared<fem::FunctionSpace<double>>(
+      fem::create_functionspace<double>(
+          smesh, std::make_shared<fem::FiniteElement<double>>(
+                     element, std::vector<std::size_t>{})));
+
+  auto kernel = [](double*, const double*, const double*, const double*,
+                   const int*, const uint8_t*, void*) {};
+  std::map integrals{
+      std::pair{std::tuple{fem::IntegralType::ridge, -1, 0},
+                fem::integral_data<double>(kernel, std::vector<std::int32_t>{},
+                                           std::vector<int>{})}};
+
+  auto build_form = [&V, &integrals, &mesh, &e_map]()
+  {
+    return fem::Form<double, double>({V}, integrals, mesh, {}, {}, false,
+                                     {std::cref(e_map)});
+  };
+
+  CHECK_THROWS_WITH(build_form(), Catch::Matchers::ContainsSubstring(
+                                      "integration entities of dimension"));
+}

python/doc/source/release_notes.md

@@ -2,6 +2,18 @@
 
 ## v0.12.0 (draft)
 
+### Assembly over ridges with data on a codimension-2 submesh
+
+**Authors**: [Jørgen S. Dokken](https://github.com/jorgensd)
+
+A form integrated over the ridges of a mesh (`ufl.dR`) may now have an
+argument or coefficient on a codimension-2 submesh, in the same way that a
+facet integral may have one on a codimension-1 submesh.
+
+Vertex integrals with data on another mesh remain unsupported and now raise a
+message naming the types that are supported.
+
+
 ### Entity permutations are computed per dimension
 
 **Authors**: [Jørgen S. Dokken](https://github.com/jorgensd)