Template over assembler kernel and use C++20 concepts (#2438)

* Use concepts for assembler kernel

* Doc update

* Alt syntax

* Use concepts

* Undo a change

* Small type fixes

* Small updates

* Simplifications

* Tidy up

* More concepts in AdjacencyList

* Add doc

* Doc fix

* Attempted work-around for old doxygen

* Re-apply fix

* concepts fix

* Update an action

* Use concepts for matrix add/set

.github/workflows/build-wheels.yml

@@ -168,13 +168,13 @@ jobs:
       # Artifact can be unzipped into $(pwd) and tested with e.g.:
       # docker run -ti -v $(pwd):/shared --env PIP_EXTRA_INDEX_URL=file:///shared python:3.9 /bin/bash -l
       - name: Upload simple503-test artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: "simple503-test"
           path: simple/*
 
       - name: Upload wheelhouse artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: wheelhouse
           path: wheelhouse/*
@@ -185,7 +185,7 @@ jobs:
 
       # For manual upload.
       - name: Upload simple503 artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: "simple503"
           path: simple/*
@@ -198,7 +198,7 @@ jobs:
           docker rm ${CONTAINER_ID}
 
       - name: Upload mpiexec artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: "mpiexec.hydra"
           path: mpiexec.hydra

.github/workflows/ccpp.yml

@@ -156,7 +156,7 @@ jobs:
         run: mpirun -np 2 python3 -m pytest python/test/unit/
 
       - name: Upload C++ documentation artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: doc-cpp-${{ matrix.petsc_arch }}-${{ matrix.petsc_int_type }}
           path: |
@@ -166,7 +166,7 @@ jobs:
           if-no-files-found: error
 
       - name: Upload Python documentation artifact
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: doc-python-${{ matrix.petsc_arch }}-${{ matrix.petsc_int_type }}
           path: |

cpp/dolfinx/common/utils.h

@@ -16,7 +16,7 @@
 namespace dolfinx::common
 {
 
-/// Sort two arrays based on the values in array @p indices. Any
+/// Sort two arrays based on the values in array `indices`. Any
 /// duplicate indices and the corresponding value are removed. In the
 /// case of duplicates, the entry with the smallest value is retained.
 /// @param[in] indices Array of indices

cpp/dolfinx/fem/DirichletBC.h

@@ -11,6 +11,7 @@
 #include "Function.h"
 #include "FunctionSpace.h"
 #include <array>
+#include <concepts>
 #include <functional>
 #include <memory>
 #include <span>
@@ -162,9 +163,8 @@ class DirichletBC
   ///
   /// @param[in] g The boundary condition value (`T` or convertible to
   /// `std::span<const T>`)
-  /// @param[in] dofs Degree-of-freedom block indices (
-  /// `std::vector<std::int32_t>`) to be constrained. The indices must
-  /// be sorted.
+  /// @param[in] dofs Degree-of-freedom block indices to be constrained.
+  /// The indices must be sorted.
   /// @param[in] V The function space to be constrained
   /// @note Can be used only with point-evaluation elements.
   /// @note The indices in `dofs` are for *blocks*, e.g. a block index
@@ -173,10 +173,10 @@ class DirichletBC
   /// @note The size of of `g` must be equal to the block size if `V`.
   /// Use the Function version if this is not the case, e.g. for some
   /// mixed spaces.
-  template <typename S, typename U,
-            typename = std::enable_if_t<
-                std::is_convertible_v<
-                    S, T> or std::is_convertible_v<S, std::span<const T>>>>
+  template <typename S, std::convertible_to<std::vector<std::int32_t>> U,
+            typename
+            = std::enable_if_t<std::is_convertible_v<S, T>
+                               or std::is_convertible_v<S, std::span<const T>>>>
   DirichletBC(const S& g, U&& dofs, std::shared_ptr<const FunctionSpace> V)
       : DirichletBC(std::make_shared<Constant<T>>(g), dofs, V)
   {
@@ -188,8 +188,7 @@ class DirichletBC
   ///@pre `dofs` must be sorted.
   ///
   /// @param[in] g The boundary condition value.
-  /// @param[in] dofs Degree-of-freedom block indices
-  /// (`std::vector<std::int32_t>`) to be constrained.
+  /// @param[in] dofs Degree-of-freedom block indices to be constrained.
   /// @param[in] V The function space to be constrained
   /// @note Can be used only with point-evaluation elements.
   /// @note The indices in `dofs` are for *blocks*, e.g. a block index
@@ -198,7 +197,7 @@ class DirichletBC
   /// @note The size of of `g` must be equal to the block size if `V`.
   /// Use the Function version if this is not the case, e.g. for some
   /// mixed spaces.
-  template <typename U>
+  template <std::convertible_to<std::vector<std::int32_t>> U>
   DirichletBC(std::shared_ptr<const Constant<T>> g, U&& dofs,
               std::shared_ptr<const FunctionSpace> V)
       : _function_space(V), _g(g), _dofs0(std::forward<U>(dofs)),
@@ -242,12 +241,11 @@ class DirichletBC
   /// @pre `dofs` must be sorted.
   ///
   /// @param[in] g The boundary condition value.
-  /// @param[in] dofs Degree-of-freedom block indices
-  /// (`std::vector<std::int32_t>`) to be constrained.
+  /// @param[in] dofs Degree-of-freedom block indices to be constrained.
   /// @note The indices in `dofs` are for *blocks*, e.g. a block index
   /// corresponds to 3 degrees-of-freedom if the dofmap associated with
   /// `g` has block size 3.
-  template <typename U>
+  template <std::convertible_to<std::vector<std::int32_t>> U>
   DirichletBC(std::shared_ptr<const Function<T>> g, U&& dofs)
       : _function_space(g->function_space()), _g(g),
         _dofs0(std::forward<U>(dofs)),
@@ -411,7 +409,7 @@ class DirichletBC
       auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
       const std::vector<T>& value = g->value;
       std::int32_t bs = _function_space->dofmap()->bs();
-      std::for_each(_dofs0.cbegin(), _dofs0.cend(),
+      std::for_each(_dofs0.begin(), _dofs0.end(),
                     [&x, &x0, &value, scale, bs](auto dof)
                     {
                       if (dof < (std::int32_t)x.size())

cpp/dolfinx/fem/DofMap.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include "ElementDofLayout.h"
+#include <concepts>
 #include <cstdlib>
 #include <dolfinx/common/MPI.h>
 #include <dolfinx/graph/AdjacencyList.h>
@@ -62,7 +63,7 @@ transpose_dofmap(const graph::AdjacencyList<std::int32_t>& dofmap,
                  std::int32_t num_cells);
 
 /// @brief Degree-of-freedom map.
-//
+///
 /// This class handles the mapping of degrees of freedom. It builds a
 /// dof map based on an ElementDofLayout on a specific mesh topology. It
 /// will reorder the dofs when running in parallel. Sub-dofmaps, both
@@ -72,21 +73,19 @@ class DofMap
 public:
   /// @brief Create a DofMap from the layout of dofs on a reference
   /// element, an IndexMap defining the distribution of dofs across
-  /// processes and a vector of indices
-  //
+  /// processes and a vector of indices.
+  ///
   /// @param[in] element The layout of the degrees of freedom on an
-  /// element (fem::ElementDofLayout)
+  /// element
   /// @param[in] index_map The map describing the parallel distribution
-  /// of the degrees of freedom
-  /// @param[in] index_map_bs The block size associated with the @p
-  /// index_map
-  /// @param[in] dofmap Adjacency list
-  /// (graph::AdjacencyList<std::int32_t>) with the degrees-of-freedom
-  /// for each cell
-  /// @param[in] bs The block size of the @p dofmap
-  template <typename E, typename U,
-            typename = std::enable_if_t<std::is_same<
-                graph::AdjacencyList<std::int32_t>, std::decay_t<U>>::value>>
+  /// of the degrees of freedom.
+  /// @param[in] index_map_bs The block size associated with the
+  /// `index_map`.
+  /// @param[in] dofmap Adjacency list with the degrees-of-freedom for
+  /// each cell.
+  /// @param[in] bs The block size of the `dofmap`.
+  template <std::convertible_to<fem::ElementDofLayout> E,
+            std::convertible_to<graph::AdjacencyList<std::int32_t>> U>
   DofMap(E&& element, std::shared_ptr<const common::IndexMap> index_map,
          int index_map_bs, U&& dofmap, int bs)
       : index_map(index_map), _index_map_bs(index_map_bs),

cpp/dolfinx/fem/Form.h

@@ -286,8 +286,8 @@ class Form
     return it->second.second;
   }
 
-  /// Get the list of (cell_index, local_facet_index) pairs for the ith
-  /// integral (kernel) for the exterior facet domain type
+  /// @brief List of (cell_index, local_facet_index) pairs for the ith
+  /// integral (kernel) for the exterior facet domain type.
   /// @param[in] i Integral ID, i.e. (sub)domain index
   /// @return List of (cell_index, local_facet_index) pairs. This data is
   /// flattened with row-major layout, shape=(num_facets, 2)

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -24,22 +24,10 @@
 namespace dolfinx::fem::impl
 {
 
-/// The matrix A must already be initialised. The matrix may be a proxy,
-/// i.e. a view into a larger matrix, and assembly is performed using
-/// local indices. Rows (bc0) and columns (bc1) with Dirichlet
-/// conditions are zeroed. Markers (bc0 and bc1) can be empty if not bcs
-/// are applied. Matrix is not finalised.
-template <typename T, typename U>
-void assemble_matrix(
-    U mat_set_values, const Form<T>& a, std::span<const T> constants,
-    const std::map<std::pair<IntegralType, int>,
-                   std::pair<std::span<const T>, int>>& coefficients,
-    std::span<const std::int8_t> bc0, std::span<const std::int8_t> bc1);
-
 /// Execute kernel over cells and accumulate result in matrix
-template <typename T, typename U>
+template <typename T>
 void assemble_cells(
-    U mat_set, const mesh::Geometry& geometry,
+    la::MatSet<T> auto mat_set, const mesh::Geometry& geometry,
     std::span<const std::int32_t> cells,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
@@ -50,19 +38,16 @@ void assemble_cells(
                              std::int32_t, int)>& dof_transform_to_transpose,
     const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
     std::span<const std::int8_t> bc0, std::span<const std::int8_t> bc1,
-    const std::function<void(T*, const T*, const T*,
-                             const scalar_value_type_t<T>*, const int*,
-                             const std::uint8_t*)>& kernel,
-    std::span<const T> coeffs, int cstride, std::span<const T> constants,
-    std::span<const std::uint32_t> cell_info)
+    FEkernel<T> auto kernel, std::span<const T> coeffs, int cstride,
+    std::span<const T> constants, std::span<const std::uint32_t> cell_info)
 {
   if (cells.empty())
     return;
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
   const std::size_t num_dofs_g = geometry.cmap().dim();
-  std::span<const double> x_g = geometry.x();
+  std::span<const double> x = geometry.x();
 
   // Iterate over active cells
   const int num_dofs0 = dofmap0.links(0).size();
@@ -82,7 +67,7 @@ void assemble_cells(
     auto x_dofs = x_dofmap.links(c);
     for (std::size_t i = 0; i < x_dofs.size(); ++i)
     {
-      std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3,
+      std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
                   std::next(coordinate_dofs.begin(), 3 * i));
     }
 
@@ -135,9 +120,10 @@ void assemble_cells(
 }
 
 /// Execute kernel over exterior facets and  accumulate result in Mat
-template <typename T, typename U>
+template <typename T>
 void assemble_exterior_facets(
-    U mat_set, const mesh::Mesh& mesh, std::span<const std::int32_t> facets,
+    la::MatSet<T> auto mat_set, const mesh::Mesh& mesh,
+    std::span<const std::int32_t> facets,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
@@ -147,19 +133,16 @@ void assemble_exterior_facets(
                              std::int32_t, int)>& dof_transform_to_transpose,
     const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
     std::span<const std::int8_t> bc0, std::span<const std::int8_t> bc1,
-    const std::function<void(T*, const T*, const T*,
-                             const scalar_value_type_t<T>*, const int*,
-                             const std::uint8_t*)>& kernel,
-    std::span<const T> coeffs, int cstride, std::span<const T> constants,
-    std::span<const std::uint32_t> cell_info)
+    FEkernel<T> auto kernel, std::span<const T> coeffs, int cstride,
+    std::span<const T> constants, std::span<const std::uint32_t> cell_info)
 {
   if (facets.empty())
     return;
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
   const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
-  std::span<const double> x_g = mesh.geometry().x();
+  std::span<const double> x = mesh.geometry().x();
 
   // Data structures used in assembly
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
@@ -179,7 +162,7 @@ void assemble_exterior_facets(
     auto x_dofs = x_dofmap.links(cell);
     for (std::size_t i = 0; i < x_dofs.size(); ++i)
     {
-      std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3,
+      std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
                   std::next(coordinate_dofs.begin(), 3 * i));
     }
 
@@ -231,9 +214,10 @@ void assemble_exterior_facets(
 }
 
 /// Execute kernel over interior facets and  accumulate result in Mat
-template <typename T, typename U>
+template <typename T>
 void assemble_interior_facets(
-    U mat_set, const mesh::Mesh& mesh, std::span<const std::int32_t> facets,
+    la::MatSet<T> auto mat_set, const mesh::Mesh& mesh,
+    std::span<const std::int32_t> facets,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
@@ -242,10 +226,7 @@ void assemble_interior_facets(
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform_to_transpose,
     const DofMap& dofmap1, int bs1, std::span<const std::int8_t> bc0,
-    std::span<const std::int8_t> bc1,
-    const std::function<void(T*, const T*, const T*,
-                             const scalar_value_type_t<T>*, const int*,
-                             const std::uint8_t*)>& kernel,
+    std::span<const std::int8_t> bc1, FEkernel<T> auto kernel,
     std::span<const T> coeffs, int cstride, std::span<const int> offsets,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info,
     const std::function<std::uint8_t(std::size_t)>& get_perm)
@@ -258,7 +239,7 @@ void assemble_interior_facets(
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
   const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
-  std::span<const double> x_g = mesh.geometry().x();
+  std::span<const double> x = mesh.geometry().x();
 
   // Data structures used in assembly
   using X = scalar_value_type_t<T>;
@@ -286,13 +267,13 @@ void assemble_interior_facets(
     auto x_dofs0 = x_dofmap.links(cells[0]);
     for (std::size_t i = 0; i < x_dofs0.size(); ++i)
     {
-      std::copy_n(std::next(x_g.begin(), 3 * x_dofs0[i]), 3,
+      std::copy_n(std::next(x.begin(), 3 * x_dofs0[i]), 3,
                   std::next(cdofs0.begin(), 3 * i));
     }
     auto x_dofs1 = x_dofmap.links(cells[1]);
     for (std::size_t i = 0; i < x_dofs1.size(); ++i)
     {
-      std::copy_n(std::next(x_g.begin(), 3 * x_dofs1[i]), 3,
+      std::copy_n(std::next(x.begin(), 3 * x_dofs1[i]), 3,
                   std::next(cdofs1.begin(), 3 * i));
     }
 
@@ -377,9 +358,14 @@ void assemble_interior_facets(
   }
 }
 
-template <typename T, typename U>
+/// The matrix A must already be initialised. The matrix may be a proxy,
+/// i.e. a view into a larger matrix, and assembly is performed using
+/// local indices. Rows (bc0) and columns (bc1) with Dirichlet
+/// conditions are zeroed. Markers (bc0 and bc1) can be empty if not bcs
+/// are applied. Matrix is not finalised.
+template <typename T>
 void assemble_matrix(
-    U mat_set, const Form<T>& a, std::span<const T> constants,
+    la::MatSet<T> auto mat_set, const Form<T>& a, std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients,
     std::span<const std::int8_t> bc0, std::span<const std::int8_t> bc1)