Put FEM 'packing' functions into separate header (#3621)

* Simplify coefficient packing for Expressions

* Work on packing

* Merge pack functions

* Simplify

* Simplify

* Break up fem/utils.h

* Add pack.h

* Add concept

* Fix docstring

* Simplifications

* Simplifications

* Small edit

Author: garth-wells

cpp/dolfinx/fem/CMakeLists.txt

@@ -18,6 +18,7 @@ set(HEADERS_fem
     ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_fem.h
     ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.h
     ${CMAKE_CURRENT_SOURCE_DIR}/petsc.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/pack.h
     ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.h
     ${CMAKE_CURRENT_SOURCE_DIR}/traits.h
     ${CMAKE_CURRENT_SOURCE_DIR}/utils.h

cpp/dolfinx/fem/Expression.h

@@ -8,6 +8,8 @@
 
 #include "Constant.h"
 #include "Function.h"
+#include "pack.h"
+#include "utils.h"
 #include <algorithm>
 #include <array>
 #include <concepts>
@@ -20,9 +22,6 @@
 
 namespace dolfinx::fem
 {
-template <dolfinx::scalar T>
-class Constant;
-
 /// @brief Represents a mathematical expression evaluated at a
 /// pre-defined set of points on the reference cell.
 ///
@@ -167,8 +166,21 @@ class Expression
       throw std::runtime_error("Invalid dimension of evaluation points.");
 
     // Prepare coefficients and constants
-    auto [coeffs, cstride] = pack_coefficients(*this, entities, estride);
-    std::vector<scalar_type> constant_data = pack_constants(*this);
+    std::vector<T> coeffs((entities.size() / estride)
+                          * this->coefficient_offsets().back());
+    int cstride = this->coefficient_offsets().back();
+    {
+      std::vector<
+          std::reference_wrapper<const Function<scalar_type, geometry_type>>>
+          c;
+      std::ranges::transform(this->coefficients(), std::back_inserter(c),
+                             [](auto c) -> const Function<T, U>&
+                             { return *c; });
+      fem::pack_coefficients(c, this->coefficient_offsets(), entities, estride,
+                             std::span(coeffs));
+    }
+    std::vector<scalar_type> constant_data = fem::pack_constants(*this);
+
     auto fn = this->get_tabulate_expression();
 
     // Prepare cell geometry

cpp/dolfinx/fem/assembler.h

@@ -9,6 +9,7 @@
 #include "assemble_matrix_impl.h"
 #include "assemble_scalar_impl.h"
 #include "assemble_vector_impl.h"
+#include "pack.h"
 #include "traits.h"
 #include "utils.h"
 #include <algorithm>
@@ -38,10 +39,10 @@ make_coefficients_span(const std::map<std::pair<IntegralType, int>,
 {
   using Key = typename std::remove_reference_t<decltype(coeffs)>::key_type;
   std::map<Key, std::pair<std::span<const T>, int>> c;
-  std::ranges::transform(coeffs, std::inserter(c, c.end()),
-                         [](auto& e) -> typename decltype(c)::value_type {
-                           return {e.first, {e.second.first, e.second.second}};
-                         });
+  std::ranges::transform(
+      coeffs, std::inserter(c, c.end()),
+      [](auto& e) -> typename decltype(c)::value_type
+      { return {e.first, {e.second.first, e.second.second}}; });
   return c;
 }