Use geometry type independently

cpp/dolfinx/fem/assemble_expression_impl.h

@@ -61,11 +61,11 @@ namespace dolfinx::fem::impl
 /// @param[in] perms Entity permutation information for use in `fn`.
 template <dolfinx::scalar T, std::floating_point U>
 void tabulate_expression(
-    std::span<T> values, fem::FEkernel<T> auto fn,
+    std::span<T> values, fem::FEkernel<T, U> auto fn,
     std::array<std::size_t, 2> Xshape, std::size_t value_size,
     std::size_t num_argument_dofs,
     md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> x_dofmap,
-    std::span<const scalar_value_t<T>> x,
+    std::span<const U> x,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, fem::MDSpan2 auto entities,
     std::span<const std::uint32_t> cell_info,
@@ -153,7 +153,7 @@ void tabulate_expression(
 /// expression values at the evaluation points.
 template <dolfinx::scalar T, std::floating_point U>
 void tabulate_expression(
-    std::span<T> values, fem::FEkernel<T> auto fn,
+    std::span<T> values, fem::FEkernel<T, U> auto fn,
     std::array<std::size_t, 2> Xshape, std::size_t value_size,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, const mesh::Mesh<U>& mesh,

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -62,18 +62,16 @@ using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
 /// function mesh.
 /// @param cell_info1 Cell permutation information for the trial
 /// function mesh.
-template <dolfinx::scalar T, bool LiftingMode = false>
+template <dolfinx::scalar T, std::floating_point U, bool LiftingMode = false>
 void assemble_cells_matrix(
     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 U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     std::span<const std::int32_t> cells,
     std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
     fem::DofTransformKernel<T> auto P0,
     std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
     fem::DofTransformKernel<T> auto P1T, std::span<const std::int8_t> bc0,
-    std::span<const std::int8_t> bc1, FEkernel<T> auto kernel,
+    std::span<const std::int8_t> bc1, FEkernel<T, U> auto kernel,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1)
@@ -90,7 +88,7 @@ void assemble_cells_matrix(
   const int ndim0 = bs0 * num_dofs0;
   const int ndim1 = bs1 * num_dofs1;
   std::vector<T> Ae(ndim0 * ndim1);
-  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+  std::vector<U> cdofs(3 * x_dofmap.extent(1));
 
   // Iterate over active cells
   assert(cells0.size() == cells.size());
@@ -226,12 +224,10 @@ void assemble_cells_matrix(
 /// function mesh.
 /// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
-template <dolfinx::scalar T, bool LiftingMode = false>
+template <dolfinx::scalar T, std::floating_point U, bool LiftingMode = false>
 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 U, 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>>
         entities,
@@ -245,7 +241,7 @@ void assemble_entities(
                           std::extents<std::size_t, md::dynamic_extent, 2>>>
         dofmap1,
     fem::DofTransformKernel<T> auto P1T, std::span<const std::int8_t> bc0,
-    std::span<const std::int8_t> bc1, FEkernel<T> auto kernel,
+    std::span<const std::int8_t> bc1, FEkernel<T, U> auto kernel,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1,
@@ -258,7 +254,7 @@ void assemble_entities(
   const auto [dmap1, bs1, entities1] = dofmap1;
 
   // Data structures used in assembly
-  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+  std::vector<U> cdofs(3 * x_dofmap.extent(1));
   const int num_dofs0 = dmap0.extent(1);
   const int num_dofs1 = dmap1.extent(1);
   const int ndim0 = bs0 * num_dofs0;
@@ -393,12 +389,10 @@ void assemble_entities(
 /// function mesh.
 /// @param[in] perms Facet permutation integer. Empty if facet
 /// permutations are not required.
-template <dolfinx::scalar T, bool LiftingMode = false>
+template <dolfinx::scalar T, std::floating_point U, bool LiftingMode = false>
 void assemble_interior_facets(
     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 U, 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, 2>>
         facets,
@@ -412,7 +406,7 @@ void assemble_interior_facets(
                           std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
         dofmap1,
     fem::DofTransformKernel<T> auto P1T, std::span<const std::int8_t> bc0,
-    std::span<const std::int8_t> bc1, FEkernel<T> auto kernel,
+    std::span<const std::int8_t> bc1, FEkernel<T, U> auto kernel,
     md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                     md::dynamic_extent>>
         coeffs,
@@ -427,7 +421,7 @@ void assemble_interior_facets(
   const auto [dmap1, bs1, facets1] = dofmap1;
 
   // Data structures used in assembly
-  using X = scalar_value_t<T>;
+  using X = U;
   std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
   std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
   std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
@@ -612,9 +606,7 @@ void assemble_interior_facets(
 template <dolfinx::scalar T, std::floating_point U, bool LiftingMode = false>
 void assemble_matrix(
     la::MatSet<T> auto mat_set, const Form<T, U>& a,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients,
@@ -687,7 +679,7 @@ void assemble_matrix(
       std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx);
       auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
       assert(cells.size() * cstride == coeffs.size());
-      impl::assemble_cells_matrix<T, LiftingMode>(
+      impl::assemble_cells_matrix<T, U, LiftingMode>(
           mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0}, P0,
           {dofs1, bs1, cells1}, P1T, bc0, bc1, fn,
           md::mdspan(coeffs.data(), cells.size(), cstride), constants,
@@ -732,7 +724,7 @@ void assemble_matrix(
       std::span facets0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0);
       std::span facets1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0);
       assert((facets.size() / 4) * 2 * cstride == coeffs.size());
-      impl::assemble_interior_facets<T, LiftingMode>(
+      impl::assemble_interior_facets<T, U, LiftingMode>(
           mat_set, x_dofmap, x,
           mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
           {*dofmap0, bs0,
@@ -777,7 +769,7 @@ void assemble_matrix(
         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<T, LiftingMode>(
+        impl::assemble_entities<T, U, LiftingMode>(
             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,

cpp/dolfinx/fem/assemble_scalar_impl.h

@@ -17,20 +17,20 @@
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/Topology.h>
 #include <memory>
+#include <type_traits>
 #include <vector>
 
 namespace dolfinx::fem::impl
 {
 /// Assemble functional over cells
-template <dolfinx::scalar T>
-T assemble_cells(mdspan2_t x_dofmap,
-                 md::mdspan<const scalar_value_t<T>,
-                            md::extents<std::size_t, md::dynamic_extent, 3>>
-                     x,
-                 std::span<const std::int32_t> cells, FEkernel<T> auto fn,
-                 std::span<const T> constants,
-                 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
-                 std::span<scalar_value_t<T>> cdofs_b)
+template <dolfinx::scalar T, std::floating_point U>
+T assemble_cells(
+    mdspan2_t x_dofmap,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
+    std::span<const std::int32_t> cells, FEkernel<T, U> auto fn,
+    std::span<const T> constants,
+    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
+    std::span<std::type_identity_t<U>> cdofs_b)
 {
   T value(0);
   if (cells.empty())
@@ -65,19 +65,17 @@ T assemble_cells(mdspan2_t x_dofmap,
 /// 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.
-template <dolfinx::scalar T>
+template <dolfinx::scalar T, std::floating_point U>
 T assemble_entities(
     mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     md::mdspan<const std::int32_t,
                md::extents<std::size_t, md::dynamic_extent, 2>>
         entities,
-    FEkernel<T> auto fn, std::span<const T> constants,
+    FEkernel<T, U> auto fn, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
-    std::span<scalar_value_t<T>> cdofs_b)
+    std::span<std::type_identity_t<U>> cdofs_b)
 {
   T value(0);
   if (entities.empty())
@@ -106,21 +104,19 @@ T assemble_entities(
 }
 
 /// Assemble functional over interior facets
-template <dolfinx::scalar T>
+template <dolfinx::scalar T, std::floating_point U>
 T assemble_interior_facets(
     mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     md::mdspan<const std::int32_t,
                md::extents<std::size_t, md::dynamic_extent, 2, 2>>
         facets,
-    FEkernel<T> auto fn, std::span<const T> constants,
+    FEkernel<T, U> auto fn, std::span<const T> constants,
     md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                     md::dynamic_extent>>
         coeffs,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms,
-    std::span<scalar_value_t<T>> cdofs_b)
+    std::span<std::type_identity_t<U>> cdofs_b)
 {
   T value(0);
   if (facets.empty())
@@ -160,9 +156,7 @@ T assemble_interior_facets(
 template <dolfinx::scalar T, std::floating_point U>
 T assemble_scalar(
     const fem::Form<T, U>& M, mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients,
@@ -171,7 +165,7 @@ T assemble_scalar(
   std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
   assert(mesh);
 
-  std::vector<scalar_value_t<T>> cdofs_b(2 * 3 * x_dofmap.extent(1));
+  std::vector<U> cdofs_b(2 * 3 * x_dofmap.extent(1));
 
   T value = 0;
   for (int i = 0; i < M.num_integrals(IntegralType::cell, cell_type_idx); ++i)

cpp/dolfinx/fem/assemble_vector_impl.h

@@ -61,18 +61,16 @@ using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
 /// coefficient for cell `i`.
 /// @param[in] cell_info0 Cell permutation information for the test
 /// function mesh.
-template <int _bs = -1, typename V,
+template <int _bs = -1, typename V, std::floating_point U,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
 
 void assemble_cells(
     fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     std::span<const std::int32_t> cells,
     std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
-    FEkernel<T> auto kernel, std::span<const T> constants,
+    FEkernel<T, U> auto kernel, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const std::uint32_t> cell_info0)
 {
@@ -83,7 +81,7 @@ void assemble_cells(
   assert(_bs < 0 or _bs == bs);
 
   // Create data structures used in assembly
-  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+  std::vector<U> cdofs(3 * x_dofmap.extent(1));
   std::vector<T> be(bs * dmap.extent(1));
 
   // Iterate over active cells
@@ -155,22 +153,20 @@ void assemble_cells(
 /// function mesh.
 /// @param[in] perms Entity permutation integer. Empty if entity
 /// permutations are not required.
-template <int _bs = -1, typename V,
+template <int _bs = -1, typename V, std::floating_point U,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
 void assemble_entities(
     fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, 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>>
         entities,
     std::tuple<mdspan2_t, int,
                md::mdspan<const std::int32_t,
                           std::extents<std::size_t, md::dynamic_extent, 2>>>
         dofmap,
-    FEkernel<T> auto kernel, std::span<const T> constants,
+    FEkernel<T, U> auto kernel, std::span<const T> constants,
     md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
     std::span<const std::uint32_t> cell_info0,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
@@ -183,7 +179,7 @@ void assemble_entities(
 
   // Create data structures used in assembly
   const int num_dofs = dmap.extent(1);
-  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+  std::vector<U> cdofs(3 * x_dofmap.extent(1));
   std::vector<T> be(bs * num_dofs);
   assert(entities0.size() == entities.size());
   for (std::size_t f = 0; f < entities.extent(0); ++f)
@@ -251,29 +247,27 @@ void assemble_entities(
 /// function mesh.
 /// @param[in] perms Facet permutation integer. Empty if facet
 /// permutations are not required.
-template <int _bs = -1, typename V,
+template <int _bs = -1, typename V, std::floating_point U,
           dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
 void assemble_interior_facets(
     fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, 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, 2>>
         facets,
     std::tuple<const DofMap&, int,
                md::mdspan<const std::int32_t,
                           std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
         dofmap,
-    FEkernel<T> auto kernel, std::span<const T> constants,
+    FEkernel<T, U> auto kernel, std::span<const T> constants,
     md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                     md::dynamic_extent>>
         coeffs,
     std::span<const std::uint32_t> cell_info0,
     md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
 {
-  using X = scalar_value_t<T>;
+  using X = U;
 
   if (facets.empty())
     return;
@@ -567,9 +561,7 @@ template <typename V, std::floating_point U,
   requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
 void assemble_vector(
     V&& b, const Form<T, U>& L,
-    md::mdspan<const scalar_value_t<T>,
-               md::extents<std::size_t, md::dynamic_extent, 3>>
-        x,
+    md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>> x,
     std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients)
@@ -765,22 +757,12 @@ void assemble_vector(
                    std::pair<std::span<const T>, int>>& coefficients)
 {
   using mdspanx3_t
-      = md::mdspan<const scalar_value_t<T>,
-                   md::extents<std::size_t, md::dynamic_extent, 3>>;
+      = md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>>;
 
   std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
   assert(mesh);
   auto x = mesh->geometry().x();
-  if constexpr (std::is_same_v<U, scalar_value_t<T>>)
-  {
-    impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3),
-                          constants, coefficients);
-  }
-  else
-  {
-    std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
-    impl::assemble_vector(b, L, mdspanx3_t(_x.data(), _x.size() / 3, 3),
-                          constants, coefficients);
-  }
+  impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3), constants,
+                        coefficients);
 }
 } // namespace dolfinx::fem::impl