Replace IntegralType enum with something that is more expressive

.github/workflows/fenicsx-refs.env

@@ -3,4 +3,5 @@ basix_ref=main
 ufl_repository=FEniCS/ufl
 ufl_ref=main
 ffcx_repository=FEniCS/ffcx
-ffcx_ref=main
+ffcx_ref=dokken/integral-type
+

.github/workflows/windows.yml

@@ -50,6 +50,7 @@ jobs:
       - name: Load environment variables
         run: |
           cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $env:GITHUB_ENV
+
       - name: Set up Python
         uses: actions/setup-python@v6
         with:
@@ -95,6 +96,7 @@ jobs:
         run: |
           python -m pip -v install --check-build-dependencies --no-build-isolation --no-cache-dir .[ci] --config-settings=cmake.args=-DBasix_DIR="D:/a/dolfinx/basix-install/lib/cmake/basix" --config-settings=cmake.build-type="Release"
 
+
       - name: Checkout FFCx
         uses: actions/checkout@v5
         with:

cpp/demo/codim_0_assembly/main.cpp

@@ -100,12 +100,11 @@ int main(int argc, char* argv[])
     // domain `mesh`
     std::vector<std::int32_t> integration_entities
         = fem::compute_integration_domains(
-            fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2));
+            fem::IntegralType(0), *mesh->topology(), cell_marker.find(2));
     std::map<
         fem::IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
-        subdomain_data
-        = {{fem::IntegralType::cell, {{3, integration_entities}}}};
+        subdomain_data = {{fem::IntegralType(0), {{3, integration_entities}}}};
 
     // A mixed-domain form involves functions defined over multiple
     // meshes. The mesh passed to `create_form` is called the

cpp/demo/custom_kernel/main.cpp

@@ -75,7 +75,7 @@ double assemble_matrix0(std::shared_ptr<const fem::FunctionSpace<T>> V,
 {
   // Kernel data (ID, kernel function, cell indices to execute over)
   std::map integrals{
-      std::pair{std::tuple{fem::IntegralType::cell, 0, 0},
+      std::pair{std::tuple{fem::IntegralType(0), 0, 0},
                 fem::integral_data<T>(kernel, cells, std::vector<int>{})}};
 
   fem::Form<T, T> a({V, V}, integrals, V->mesh(), {}, {}, false, {});
@@ -105,7 +105,7 @@ double assemble_vector0(std::shared_ptr<const fem::FunctionSpace<T>> V,
 {
   auto mesh = V->mesh();
   std::map integrals{
-      std::pair{std::tuple{fem::IntegralType::cell, 0, 0},
+      std::pair{std::tuple{fem::IntegralType(0), 0, 0},
                 fem::integral_data<T>(kernel, cells, std::vector<int>{})}};
   fem::Form<T> L({V}, integrals, mesh, {}, {}, false, {});
   auto dofmap = V->dofmap();

cpp/demo/mixed_poisson/main.cpp

@@ -271,7 +271,7 @@ int main(int argc, char* argv[])
     // (applied on the `ds(1)` in the UFL file). First we get facet data
     // integration data for facets in dfacets.
     std::vector<std::int32_t> domains = fem::compute_integration_domains(
-        fem::IntegralType::exterior_facet, *mesh->topology(), dfacets);
+        fem::IntegralType(1), *mesh->topology(), dfacets);
 
     // Create data structure for the `ds(1)` integration domain in form
     // (see the UFL file). It is for en exterior facet integral (the key
@@ -281,7 +281,7 @@ int main(int argc, char* argv[])
     std::map<
         fem::IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
-        subdomain_data{{fem::IntegralType::exterior_facet, {{1, domains}}}};
+        subdomain_data{{fem::IntegralType(1), {{1, domains}}}};
 
     // Since we are doing a `ds(1)` integral on mesh and `u0` is defined
     // on the `submesh`, our form involves more than one mesh. The mesh

cpp/dolfinx/fem/Form.h

@@ -35,12 +35,44 @@ template <dolfinx::scalar T, std::floating_point U>
 class Function;
 
 /// @brief Type of integral
-enum class IntegralType : std::int8_t
+
+struct IntegralType
 {
-  cell = 0,           ///< Cell
-  exterior_facet = 1, ///< Exterior facet
-  interior_facet = 2, ///< Interior facet
-  vertex = 3          ///< Vertex
+  std::int32_t codim;     ///< Codimension of the integral
+  std::int32_t num_cells; ///< Number of cells in the integral
+
+  /// @brief Create an integral type.
+  /// @param codim Codimension of the integral.
+  /// @param num_cells Number of cells in the integral. Default is 1.
+  IntegralType(std::int32_t codim, std::int32_t num_cells = 1)
+      : codim(codim), num_cells(num_cells)
+  {
+  }
+
+  /// @brief Equality operator for integral type
+  /// @param other The other IntegralType to compare with
+  /// @return True if the integral types are equal, false otherwise
+  bool operator==(const IntegralType& other) const
+  {
+    return codim == other.codim && num_cells == other.num_cells;
+  }
+
+  /// @brief Less than operator.
+  /// Required for set operations.
+  /// @param other The other IntegralType to compare with
+  bool operator<(const IntegralType& other) const
+  {
+    // First, compare codim
+    if (codim < other.codim)
+      return true;
+
+    // If codims are equal, compare num_cells
+    if (codim == other.codim)
+      return num_cells < other.num_cells;
+
+    // Otherwise, this object is not less than 'other'
+    return false;
+  }
 };
 
 /// @brief Represents integral data, containing the kernel, and a list
@@ -130,9 +162,9 @@ class Form
   /// trial function spaces.
   /// @param[in] integrals Integrals in the form, where
   /// `integrals[IntegralType, i, kernel index]` returns the `i`th integral
-  /// (`integral_data`) of type `IntegralType` with kernel index `kernel index`.
-  /// The `i`-index refers to the position of a kernel when flattened by
-  /// sorted subdomain ids, sorted by subdomain ids. The subdomain ids can
+  /// (`integral_data`) of type `IntegralType` with kernel index `kernel
+  /// index`. The `i`-index refers to the position of a kernel when flattened
+  /// by sorted subdomain ids, sorted by subdomain ids. The subdomain ids can
   /// contain duplicate entries referring to different kernels over the same
   /// subdomain.
   /// @param[in] coefficients Coefficients in the form.
@@ -274,10 +306,9 @@ class Form
         {
           auto [type, idx, kernel_idx] = key;
           std::vector<std::int32_t> e;
-          if (type == IntegralType::cell)
+          if (type.codim == 0)
             e = emap.sub_topology_to_topology(itg.entities, inverse);
-          else if (type == IntegralType::exterior_facet
-                   or type == IntegralType::interior_facet)
+          else if (type.codim == 1)
           {
             const mesh::Topology topology = *_mesh->topology();
             int tdim = topology.dim();
@@ -317,10 +348,9 @@ class Form
           bool inverse = emap.sub_topology() == mesh0->topology();
 
           std::vector<std::int32_t> e;
-          if (type == IntegralType::cell)
+          if (type.codim == 0)
             e = emap.sub_topology_to_topology(integral.entities, inverse);
-          else if (type == IntegralType::exterior_facet
-                   or type == IntegralType::interior_facet)
+          else if (type.codim == 1)
           {
             const mesh::Topology topology = *_mesh->topology();
             int tdim = topology.dim();
@@ -454,16 +484,12 @@ class Form
   /// These are the entities in the mesh returned by ::mesh that are
   /// integrated over by a given integral (kernel).
   ///
-  /// - For IntegralType::cell, returns a list of cell indices.
-  /// - For IntegralType::exterior_facet, returns a list with shape
-  /// `(num_facets, 2)`, where `[cell_index, 0]` is the cell index and
-  /// `[cell_index, 1]` is the local facet index relative to the cell.
-  /// - For IntegralType::interior_facet the shape is `(num_facets, 4)`,
-  /// where `[cell_index, 0]` is one attached cell and `[cell_index, 1]`
-  /// is the is the local facet index relative to the cell, and
-  /// `[cell_index, 2]` is the other one attached cell and `[cell_index, 1]`
-  /// is the is the local facet index relative to this cell. Storage
-  /// is row-major.
+  /// - For IntegralType of codim 0 returns a list of cell indices.
+  /// - For any other codim return a set of tuples `[cell_index,
+  /// local_entity_index]`.
+  /// - If the number of cells in the integral type, return as many tuples per
+  /// entity as there are number of cells in the integral type. Storage is
+  /// row-major.
   ///
   /// @param[in] type Integral type.
   /// @param[in] idx Integral index in flattened list of integral kernels.
@@ -512,8 +538,8 @@ class Form
   /// is marked with -1.
   ///
   /// @param[in] type Integral type.
-  /// @param[in] rank Argument index, e.g. `0` for the test function space, `1`
-  /// for the trial function space.
+  /// @param[in] rank Argument index, e.g. `0` for the test function space,
+  /// `1` for the trial function space.
   /// @param[in] idx Integral identifier.
   /// @param[in] kernel_idx Kernel index (cell type).
   /// @return Entity indices in the argument function space mesh that is

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -588,15 +588,15 @@ void assemble_matrix(
       cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
       cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
     }
-
-    for (int i = 0; i < a.num_integrals(IntegralType::cell, cell_type_idx); ++i)
+    IntegralType cell_integral = IntegralType(0);
+    for (int i = 0; i < a.num_integrals(cell_integral, cell_type_idx); ++i)
     {
-      auto fn = a.kernel(IntegralType::cell, i, cell_type_idx);
+      auto fn = a.kernel(cell_integral, i, cell_type_idx);
       assert(fn);
-      std::span cells = a.domain(IntegralType::cell, i, cell_type_idx);
-      std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
-      std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx);
-      auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
+      std::span cells = a.domain(cell_integral, i, cell_type_idx);
+      std::span cells0 = a.domain_arg(cell_integral, 0, i, cell_type_idx);
+      std::span cells1 = a.domain_arg(cell_integral, 1, i, cell_type_idx);
+      auto& [coeffs, cstride] = coefficients.at({cell_integral, i});
       assert(cells.size() * cstride == coeffs.size());
       impl::assemble_cells_matrix(
           mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0}, P0,
@@ -618,8 +618,8 @@ void assemble_matrix(
                          num_facets_per_cell);
     }
 
-    for (int i = 0;
-         i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i)
+    IntegralType facet = IntegralType(1);
+    for (int i = 0; i < a.num_integrals(facet, cell_type_idx); ++i)
     {
       if (num_cell_types > 1)
       {
@@ -631,16 +631,15 @@ void assemble_matrix(
           = md::mdspan<const std::int32_t,
                        md::extents<std::size_t, md::dynamic_extent, 2>>;
 
-      auto fn = a.kernel(IntegralType::exterior_facet, i, 0);
+      auto fn = a.kernel(facet, i, 0);
       assert(fn);
-      auto& [coeffs, cstride]
-          = coefficients.at({IntegralType::exterior_facet, i});
+      auto& [coeffs, cstride] = coefficients.at({facet, i});
 
-      std::span f = a.domain(IntegralType::exterior_facet, i, 0);
+      std::span f = a.domain(facet, i, 0);
       mdspanx2_t facets(f.data(), f.size() / 2, 2);
-      std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0);
+      std::span f0 = a.domain_arg(facet, 0, i, 0);
       mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
-      std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0);
+      std::span f1 = a.domain_arg(facet, 1, i, 0);
       mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
       assert((facets.size() / 2) * cstride == coeffs.size());
       impl::assemble_exterior_facets(
@@ -650,8 +649,8 @@ void assemble_matrix(
           cell_info0, cell_info1, perms);
     }
 
-    for (int i = 0;
-         i < a.num_integrals(IntegralType::interior_facet, cell_type_idx); ++i)
+    IntegralType interior_facet = IntegralType(1, 2);
+    for (int i = 0; i < a.num_integrals(interior_facet, cell_type_idx); ++i)
     {
       if (num_cell_types > 1)
       {
@@ -666,14 +665,13 @@ void assemble_matrix(
           = md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                             md::dynamic_extent>>;
 
-      auto fn = a.kernel(IntegralType::interior_facet, i, 0);
+      auto fn = a.kernel(interior_facet, i, 0);
       assert(fn);
-      auto& [coeffs, cstride]
-          = coefficients.at({IntegralType::interior_facet, i});
+      auto& [coeffs, cstride] = coefficients.at({interior_facet, i});
 
-      std::span facets = a.domain(IntegralType::interior_facet, i, 0);
-      std::span facets0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0);
-      std::span facets1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0);
+      std::span facets = a.domain(interior_facet, i, 0);
+      std::span facets0 = a.domain_arg(interior_facet, 0, i, 0);
+      std::span facets1 = a.domain_arg(interior_facet, 1, i, 0);
       assert((facets.size() / 4) * 2 * cstride == coeffs.size());
       impl::assemble_interior_facets(
           mat_set, x_dofmap, x,

cpp/dolfinx/fem/assemble_scalar_impl.h

@@ -201,12 +201,13 @@ T assemble_scalar(
   std::vector<scalar_value_t<T>> cdofs_b(2 * 3 * x_dofmap.extent(1));
 
   T value = 0;
-  for (int i = 0; i < M.num_integrals(IntegralType::cell, 0); ++i)
+  IntegralType cell_integral = IntegralType(0);
+  for (int i = 0; i < M.num_integrals(cell_integral, 0); ++i)
   {
-    auto fn = M.kernel(IntegralType::cell, i, 0);
+    auto fn = M.kernel(cell_integral, i, 0);
     assert(fn);
-    auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
-    std::span<const std::int32_t> cells = M.domain(IntegralType::cell, i, 0);
+    auto& [coeffs, cstride] = coefficients.at({cell_integral, i});
+    std::span<const std::int32_t> cells = M.domain(cell_integral, i, 0);
     assert(cells.size() * cstride == coeffs.size());
     value += impl::assemble_cells(
         x_dofmap, x, cells, fn, constants,
@@ -226,14 +227,14 @@ T assemble_scalar(
                        num_facets_per_cell);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i)
+  IntegralType facet_type = IntegralType(1);
+  for (int i = 0; i < M.num_integrals(facet_type, 0); ++i)
   {
-    auto fn = M.kernel(IntegralType::exterior_facet, i, 0);
+    auto fn = M.kernel(facet_type, i, 0);
     assert(fn);
-    auto& [coeffs, cstride]
-        = coefficients.at({IntegralType::exterior_facet, i});
+    auto& [coeffs, cstride] = coefficients.at({facet_type, i});
 
-    std::span facets = M.domain(IntegralType::exterior_facet, i, 0);
+    std::span facets = M.domain(facet_type, i, 0);
 
     // Two values per each adjacent cell (cell index and local facet
     // index)
@@ -251,13 +252,13 @@ T assemble_scalar(
         cdofs_b);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i)
+  IntegralType interior_facet_type = IntegralType(1, 2);
+  for (int i = 0; i < M.num_integrals(interior_facet_type, 0); ++i)
   {
-    auto fn = M.kernel(IntegralType::interior_facet, i, 0);
+    auto fn = M.kernel(interior_facet_type, i, 0);
     assert(fn);
-    auto& [coeffs, cstride]
-        = coefficients.at({IntegralType::interior_facet, i});
-    std::span facets = M.domain(IntegralType::interior_facet, i, 0);
+    auto& [coeffs, cstride] = coefficients.at({interior_facet_type, i});
+    std::span facets = M.domain(interior_facet_type, i, 0);
 
     constexpr std::size_t num_adjacent_cells = 2;
     // Two values per each adj. cell (cell index and local facet index).
@@ -276,16 +277,16 @@ T assemble_scalar(
         perms, cdofs_b);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i)
+  IntegralType vertex_type = IntegralType(-1);
+  for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i)
   {
-    auto fn = M.kernel(IntegralType::vertex, i, 0);
+    auto fn = M.kernel(vertex_type, i, 0);
     assert(fn);
 
-    auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
+    auto& [coeffs, cstride] = coefficients.at({vertex_type, i});
 
-    std::span<const std::int32_t> vertices
-        = M.domain(IntegralType::vertex, i, 0);
-    assert(vertices.size() * cstride == coeffs.size());
+    std::span<const std::int32_t> vertices = M.domain(vertex_type, i, 0);
+    assert(vertices.size() / 2 * cstride == coeffs.size());
 
     constexpr std::size_t num_adjacent_cells = 1;
     // Two values per adj. cell (cell index and local vertex index).