Replace IntegralType with structure storing the necessary information

Author: jorgensd

cpp/demo/codim_0_assembly/main.cpp

@@ -100,12 +100,12 @@ 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, 1), *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}}}};
+        = {{fem::IntegralType(0, 1), {{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, 1), 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, 1), 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,17 +271,17 @@ 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::facet, *mesh->topology(), dfacets);
+        IntegralType(1, 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
     // in the map), and exterior facet domain marked as '1' in the UFL
     // file, and `domains` holds the necessary data to perform
     // integration of selected facets.
     std::map<
-        fem::IntegralType,
+        std::int32_t,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
-        subdomain_data{{fem::IntegralType::facet, {{1, domains}}}};
+        subdomain_data{{IntegralType(1, 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/CMakeLists.txt

@@ -34,6 +34,7 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/DofMap.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/ElementDofLayout.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElement.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/Form.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/dofmapbuilder.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/petsc.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.cpp

cpp/dolfinx/fem/Form.cpp

@@ -0,0 +1,24 @@
+// Copyright (C) 2025 Jørgen S. Dokken
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include "Form.h"
+
+bool dolfinx::fem::operator<(const dolfinx::fem::IntegralType& self,
+                             const dolfinx::fem::IntegralType& other)
+{
+  // First, compare codim
+  if (self.codim < other.codim)
+  {
+    return true;
+  }
+  // If codims are equal, compare num_cells
+  if (self.codim == other.codim)
+  {
+    return self.num_cells < other.num_cells;
+  }
+  // Otherwise, this object is not less than 'other'
+  return false;
+}

cpp/dolfinx/fem/Form.h

@@ -35,14 +35,24 @@ 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
-  facet = 1,          ///< 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 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;
+  }
 };
 
+// Overload the less-than operator
+bool operator<(const IntegralType& self, const IntegralType& other);
+
 /// @brief Represents integral data, containing the kernel, and a list
 /// of entities to integrate over and the indicies of the coefficient
 /// functions (relative to the Form) active for this integral.
@@ -274,10 +284,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::facet
-                   or type == IntegralType::interior_facet)
+          else if (type.codim = 1)
           {
             const mesh::Topology topology = *_mesh->topology();
             int tdim = topology.dim();
@@ -317,10 +326,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::facet
-                   or type == IntegralType::interior_facet)
+          else if (type.codim == 1)
           {
             const mesh::Topology topology = *_mesh->topology();
             int tdim = topology.dim();
@@ -454,16 +462,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::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.

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, 1);
+    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::facet, cell_type_idx);
-         ++i)
+    IntegralType facet = IntegralType(1, 1);
+    for (int i = 0; i < a.num_integrals(facet, cell_type_idx); ++i)
     {
       if (num_cell_types > 1)
       {
@@ -631,15 +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::facet, i, 0);
+      auto fn = a.kernel(facet, i, 0);
       assert(fn);
-      auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i});
+      auto& [coeffs, cstride] = coefficients.at({facet, i});
 
-      std::span f = a.domain(IntegralType::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::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::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(
@@ -649,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)
       {
@@ -665,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, 1);
+  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,13 +227,14 @@ T assemble_scalar(
                        num_facets_per_cell);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::facet, 0); ++i)
+  IntegralType facet_type = IntegralType(1, 1);
+  for (int i = 0; i < M.num_integrals(facet_type, 0); ++i)
   {
-    auto fn = M.kernel(IntegralType::facet, i, 0);
+    auto fn = M.kernel(facet_type, i, 0);
     assert(fn);
-    auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i});
+    auto& [coeffs, cstride] = coefficients.at({facet_type, i});
 
-    std::span facets = M.domain(IntegralType::facet, i, 0);
+    std::span facets = M.domain(facet_type, i, 0);
 
     // Two values per each adjacent cell (cell index and local facet
     // index)
@@ -250,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).
@@ -275,15 +277,15 @@ T assemble_scalar(
         perms, cdofs_b);
   }
 
-  for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i)
+  IntegralType vertex_type = IntegralType(1, 2);
+  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);
+    std::span<const std::int32_t> vertices = M.domain(vertex_type, i, 0);
     assert(vertices.size() * cstride == coeffs.size());
 
     constexpr std::size_t num_adjacent_cells = 1;