Unification of one-sided sub-entity integrals + ridge integrals (#3900)

* Unify assemblers of vertex, exterior facet and ridge integrals (new)

* Assemble scalar added

* Use vertex_perms

* Fix docs

* Fix correct usage of perms

* add 3D vector test

* Remove code duplication

* Huge duplication removal

* Remove more duplication

* Even more duplication

* Use real dtype for basix element

* Use switch

* Use more switches

* Add default

* add vertex and ridge matrix integrals + lifting

* Adapt lifting test to check matrix and lifting of ridge and vertex integrals

* Merge main into branch

* Improve name of inner-most assembly function

* Apply suggestions from code review

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

* Rename "IntegralType::exterior_facet" to "IntegralType::facet" + address other comments by Garth

* Fix spacing

* Use renaming

* Simplify

* Rename lifting

* Remove outdated comment

* Fix typo

* Explicit lambda function declaration

* Revert IntegralType::facet backt to IntegralType::exterior_facet

* More reverting

* More reverting

* Revert spacing

* Revert spacing

* Further reverting

* Yet another revert

* Apply suggestions from code review

Co-authored-by: Michal Habera <michal.habera@gmail.com>

* Apply suggestions from code review

Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com>

---------

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>
Co-authored-by: Michal Habera <michal.habera@gmail.com>
Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com>

Author: 4 people

cpp/dolfinx/common/sort.h

@@ -67,7 +67,7 @@ struct __radix_sort
   /// @code
   /// std::array<std::int16_t, 3> a{2, 3, 1};
   /// std::array<std::int16_t, 3> i{0, 1, 2};
-  /// dolfixn::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0,
+  /// dolfinx::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0,
   /// 1} and a[i] = {1, 2, 3};
   /// @endcode
   /// @tparam R Type of range to be sorted.

cpp/dolfinx/fem/Form.h

@@ -38,9 +38,10 @@ class Function;
 enum class IntegralType : std::int8_t
 {
   cell = 0,           ///< Cell
-  exterior_facet = 1, ///< Exterior facet
+  exterior_facet = 1, ///< Facet
   interior_facet = 2, ///< Interior facet
-  vertex = 3          ///< Vertex
+  vertex = 3,         ///< Vertex
+  ridge = 4           ///< Ridge
 };
 
 /// @brief Represents integral data, containing the kernel, and a list

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -156,16 +156,24 @@ void assemble_cells_matrix(
   }
 }
 
-/// @brief Execute kernel over exterior facets and accumulate result in
-/// a matrix.
+/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result
+/// in a matrix.
+///
+/// Each entity is represented by (i) a cell that the entity is attached to
+/// and (ii) the local index of the entity  with respect to the cell. The
+/// kernel is executed for each entity. The kernel can access data
+/// (e.g., coefficients, basis functions) associated with the attached cell.
+/// However, entities may be attached to more than one cell. This function
+/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen
+/// from cell used to define the entity.
 ///
 /// @tparam T Matrix/form scalar type.
 /// @param[in] mat_set Function that accumulates computed entries into a
 /// matrix.
 /// @param[in] x_dofmap Dofmap for the mesh geometry.
 /// @param[in] x Mesh geometry (coordinates).
-/// @param[in] facets Facet indices (in the integration domain mesh) to
-/// execute the kernel over.
+/// @param[in] entities Integration entities (in the integration domain mesh) to
+/// execute the kernel over. These are pairs (cell, local entity index)
 /// @param[in] dofmap0 Test function (row) degree-of-freedom data
 /// holding the (0) dofmap, (1) dofmap block size and (2) dofmap cell
 /// indices.
@@ -188,17 +196,17 @@ void assemble_cells_matrix(
 /// function mesh.
 /// @param[in] cell_info1 Cell permutation information for the trial
 /// function mesh.
-/// @param[in] perms Facet permutation integer. Empty if facet
+/// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
 template <dolfinx::scalar T>
-void assemble_exterior_facets(
+void assemble_entities(
     la::MatSet<T> auto mat_set, mdspan2_t x_dofmap,
     md::mdspan<const scalar_value_t<T>,
                md::extents<std::size_t, md::dynamic_extent, 3>>
         x,
     md::mdspan<const std::int32_t,
                std::extents<std::size_t, md::dynamic_extent, 2>>
-        facets,
+        entities,
     std::tuple<mdspan2_t, int,
                md::mdspan<const std::int32_t,
                           std::extents<std::size_t, md::dynamic_extent, 2>>>
@@ -215,11 +223,11 @@ void assemble_exterior_facets(
     std::span<const std::uint32_t> cell_info1,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
 {
-  if (facets.empty())
+  if (entities.empty())
     return;
 
-  const auto [dmap0, bs0, facets0] = dofmap0;
-  const auto [dmap1, bs1, facets1] = dofmap1;
+  const auto [dmap0, bs0, entities0] = dofmap0;
+  const auto [dmap1, bs1, entities1] = dofmap1;
 
   // Data structures used in assembly
   std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
@@ -228,30 +236,30 @@ void assemble_exterior_facets(
   const int ndim0 = bs0 * num_dofs0;
   const int ndim1 = bs1 * num_dofs1;
   std::vector<T> Ae(ndim0 * ndim1);
-  assert(facets0.size() == facets.size());
-  assert(facets1.size() == facets.size());
-  for (std::size_t f = 0; f < facets.extent(0); ++f)
+  assert(entities0.size() == entities.size());
+  assert(entities1.size() == entities.size());
+  for (std::size_t f = 0; f < entities.extent(0); ++f)
   {
-    // Cell in the integration domain, local facet index relative to the
+    // Cell in the integration domain, local entity index relative to the
     // integration domain cell, and cells in the test and trial function
     // meshes
-    std::int32_t cell = facets(f, 0);
-    std::int32_t local_facet = facets(f, 1);
-    std::int32_t cell0 = facets0(f, 0);
-    std::int32_t cell1 = facets1(f, 0);
+    std::int32_t cell = entities(f, 0);
+    std::int32_t local_entity = entities(f, 1);
+    std::int32_t cell0 = entities0(f, 0);
+    std::int32_t cell1 = entities1(f, 0);
 
     // Get cell coordinates/geometry
     auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
     for (std::size_t i = 0; i < x_dofs.size(); ++i)
       std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
 
     // Permutations
-    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
+    std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity);
 
     // Tabulate tensor
     std::ranges::fill(Ae, 0);
     kernel(Ae.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
-           &local_facet, &perm, nullptr);
+           &local_entity, &perm, nullptr);
     P0(Ae, cell_info0, cell0, ndim1);
     P1T(Ae, cell_info1, cell1, ndim0);
 
@@ -605,7 +613,7 @@ void assemble_matrix(
           cell_info0, cell_info1);
     }
 
-    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
+    md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> facet_perms;
     if (a.needs_facet_permutations())
     {
       mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
@@ -614,40 +622,8 @@ void assemble_matrix(
       mesh->topology_mutable()->create_entity_permutations();
       const std::vector<std::uint8_t>& p
           = mesh->topology()->get_facet_permutations();
-      perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
-                         num_facets_per_cell);
-    }
-
-    for (int i = 0;
-         i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i)
-    {
-      if (num_cell_types > 1)
-      {
-        throw std::runtime_error("Exterior facet integrals with mixed "
-                                 "topology aren't supported yet");
-      }
-
-      using mdspanx2_t
-          = md::mdspan<const std::int32_t,
-                       md::extents<std::size_t, md::dynamic_extent, 2>>;
-
-      auto fn = a.kernel(IntegralType::exterior_facet, i, 0);
-      assert(fn);
-      auto& [coeffs, cstride]
-          = coefficients.at({IntegralType::exterior_facet, i});
-
-      std::span f = a.domain(IntegralType::exterior_facet, i, 0);
-      mdspanx2_t facets(f.data(), f.size() / 2, 2);
-      std::span f0 = a.domain_arg(IntegralType::exterior_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);
-      mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
-      assert((facets.size() / 2) * cstride == coeffs.size());
-      impl::assemble_exterior_facets(
-          mat_set, x_dofmap, x, facets, {dofs0, bs0, facets0}, P0,
-          {dofs1, bs1, facets1}, P1T, bc0, bc1, fn,
-          md::mdspan(coeffs.data(), facets.extent(0), cstride), constants,
-          cell_info0, cell_info1, perms);
+      facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
+                               num_facets_per_cell);
     }
 
     for (int i = 0;
@@ -685,7 +661,47 @@ void assemble_matrix(
            mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
           P1T, bc0, bc1, fn,
           mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), constants,
-          cell_info0, cell_info1, perms);
+          cell_info0, cell_info1, facet_perms);
+    }
+
+    for (auto itg_type : {fem::IntegralType::exterior_facet,
+                          fem::IntegralType::vertex, fem::IntegralType::ridge})
+    {
+      md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms
+          = itg_type == fem::IntegralType::exterior_facet
+                ? facet_perms
+                : md::mdspan<const std::uint8_t,
+                             md::dextents<std::size_t, 2>>{};
+
+      for (int i = 0; i < a.num_integrals(itg_type, cell_type_idx); ++i)
+      {
+        if (num_cell_types > 1)
+        {
+          throw std::runtime_error("Exterior facet integrals with mixed "
+                                   "topology aren't supported yet");
+        }
+
+        using mdspanx2_t
+            = md::mdspan<const std::int32_t,
+                         md::extents<std::size_t, md::dynamic_extent, 2>>;
+
+        auto fn = a.kernel(itg_type, i, 0);
+        assert(fn);
+        auto& [coeffs, cstride] = coefficients.at({itg_type, i});
+
+        std::span e = a.domain(itg_type, i, 0);
+        mdspanx2_t entities(e.data(), e.size() / 2, 2);
+        std::span e0 = a.domain_arg(itg_type, 0, i, 0);
+        mdspanx2_t entities0(e0.data(), e0.size() / 2, 2);
+        std::span e1 = a.domain_arg(itg_type, 1, i, 0);
+        mdspanx2_t entities1(e1.data(), e1.size() / 2, 2);
+        assert((entities.size() / 2) * cstride == coeffs.size());
+        impl::assemble_entities(
+            mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0,
+            {dofs1, bs1, entities1}, P1T, bc0, bc1, fn,
+            md::mdspan(coeffs.data(), entities.extent(0), cstride), constants,
+            cell_info0, cell_info1, perms);
+      }
     }
   }
 }