Template mesh::Geometry over scalar type

cpp/demo/poisson/main.cpp

@@ -229,7 +229,7 @@ int main(int argc, char* argv[])
 
     b.set(0.0);
     fem::assemble_vector(b.mutable_array(), *L);
-    fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, 1.0);
+    fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, T(1));
     b.scatter_rev(std::plus<T>());
     fem::set_bc(b.mutable_array(), {bc});
 

cpp/demo/poisson_matrix_free/main.cpp

@@ -164,14 +164,14 @@ int main(int argc, char* argv[])
 
     // Apply lifting to account for Dirichlet boundary condition
     // b <- b - A * x_bc
-    fem::set_bc(ui->x()->mutable_array(), {bc}, -1.0);
+    fem::set_bc(ui->x()->mutable_array(), {bc}, T(-1));
     fem::assemble_vector(b.mutable_array(), *M);
 
     // Communicate ghost values
     b.scatter_rev(std::plus<T>());
 
     // Set BC dofs to zero (effectively zeroes columns of A)
-    fem::set_bc(b.mutable_array(), {bc}, 0.0);
+    fem::set_bc(b.mutable_array(), {bc}, T(0));
 
     b.scatter_fwd();
 
@@ -196,7 +196,7 @@ int main(int argc, char* argv[])
                            fem::make_coefficients_span(coeff));
 
       // Set BC dofs to zero (effectively zeroes rows of A)
-      fem::set_bc(y.mutable_array(), {bc}, 0.0);
+      fem::set_bc(y.mutable_array(), {bc}, T(0));
 
       // Accumuate ghost values
       y.scatter_rev(std::plus<T>());
@@ -210,7 +210,7 @@ int main(int argc, char* argv[])
     int num_it = linalg::cg(*u->x(), b, action, 200, 1e-6);
 
     // Set BC values in the solution vectors
-    fem::set_bc(u->x()->mutable_array(), {bc}, 1.0);
+    fem::set_bc(u->x()->mutable_array(), {bc}, T(1));
 
     // Compute L2 error (squared) of the solution vector e = (u - u_d, u
     // - u_d)*dx

cpp/dolfinx/fem/DirichletBC.h

@@ -348,7 +348,7 @@ class DirichletBC
   /// of the array @p x should be equal to the number of dofs owned by
   /// this rank.
   /// @param[in] scale The scaling value to apply
-  void set(std::span<T> x, double scale = 1.0) const
+  void set(std::span<T> x, T scale = 1) const
   {
     if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
     {
@@ -385,7 +385,7 @@ class DirichletBC
   /// @param[in] x The array in which to set `scale * (x0 - x_bc)`
   /// @param[in] x0 The array used in compute the value to set
   /// @param[in] scale The scaling value to apply
-  void set(std::span<T> x, std::span<const T> x0, double scale = 1.0) const
+  void set(std::span<T> x, std::span<const T> x0, T scale = 1) const
   {
     if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
     {

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -27,7 +27,8 @@ namespace dolfinx::fem::impl
 /// Execute kernel over cells and accumulate result in matrix
 template <typename T>
 void assemble_cells(
-    la::MatSet<T> auto mat_set, const mesh::Geometry& geometry,
+    la::MatSet<T> auto mat_set,
+    const mesh::Geometry<scalar_value_type_t<T>> geometry,
     std::span<const std::int32_t> cells,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
@@ -47,7 +48,7 @@ void assemble_cells(
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Iterate over active cells
   const int num_dofs0 = dofmap0.links(0).size();
@@ -122,7 +123,8 @@ void assemble_cells(
 /// Execute kernel over exterior facets and  accumulate result in Mat
 template <typename T>
 void assemble_exterior_facets(
-    la::MatSet<T> auto mat_set, const mesh::Geometry& geometry,
+    la::MatSet<T> auto mat_set,
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry,
     std::span<const std::int32_t> facets,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
@@ -142,7 +144,7 @@ void assemble_exterior_facets(
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Data structures used in assembly
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
@@ -216,8 +218,9 @@ void assemble_exterior_facets(
 /// Execute kernel over interior facets and  accumulate result in Mat
 template <typename T>
 void assemble_interior_facets(
-    la::MatSet<T> auto mat_set, const mesh::Geometry& geometry,
-    int num_cell_facets, std::span<const std::int32_t> facets,
+    la::MatSet<T> auto mat_set,
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry, int num_cell_facets,
+    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,
@@ -237,7 +240,7 @@ void assemble_interior_facets(
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Data structures used in assembly
   using X = scalar_value_type_t<T>;
@@ -360,7 +363,9 @@ void assemble_interior_facets(
 /// are applied. Matrix is not finalised.
 template <typename T>
 void assemble_matrix(
-    la::MatSet<T> auto mat_set, const Form<T>& a, std::span<const T> constants,
+    la::MatSet<T> auto mat_set, const Form<T>& a,
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    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)
@@ -409,9 +414,9 @@ void assemble_matrix(
     const auto& fn = a.kernel(IntegralType::cell, i);
     const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
     const std::vector<std::int32_t>& cells = a.cell_domains(i);
-    impl::assemble_cells(mat_set, mesh->geometry(), cells, dof_transform, dofs0,
-                         bs0, dof_transform_to_transpose, dofs1, bs1, bc0, bc1,
-                         fn, coeffs, cstride, constants, cell_info);
+    impl::assemble_cells(mat_set, geometry, cells, dof_transform, dofs0, bs0,
+                         dof_transform_to_transpose, dofs1, bs1, bc0, bc1, fn,
+                         coeffs, cstride, constants, cell_info);
   }
 
   for (int i : a.integral_ids(IntegralType::exterior_facet))
@@ -420,10 +425,10 @@ void assemble_matrix(
     const auto& [coeffs, cstride]
         = coefficients.at({IntegralType::exterior_facet, i});
     const std::vector<std::int32_t>& facets = a.exterior_facet_domains(i);
-    impl::assemble_exterior_facets(
-        mat_set, mesh->geometry(), facets, dof_transform, dofs0, bs0,
-        dof_transform_to_transpose, dofs1, bs1, bc0, bc1, fn, coeffs, cstride,
-        constants, cell_info);
+    impl::assemble_exterior_facets(mat_set, geometry, facets, dof_transform,
+                                   dofs0, bs0, dof_transform_to_transpose,
+                                   dofs1, bs1, bc0, bc1, fn, coeffs, cstride,
+                                   constants, cell_info);
   }
 
   if (a.num_integrals(IntegralType::interior_facet) > 0)
@@ -449,9 +454,9 @@ void assemble_matrix(
           = coefficients.at({IntegralType::interior_facet, i});
       const std::vector<std::int32_t>& facets = a.interior_facet_domains(i);
       impl::assemble_interior_facets(
-          mat_set, mesh->geometry(), num_cell_facets, facets, dof_transform,
-          *dofmap0, bs0, dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1,
-          fn, coeffs, cstride, c_offsets, constants, cell_info, get_perm);
+          mat_set, geometry, num_cell_facets, facets, dof_transform, *dofmap0,
+          bs0, dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1, fn, coeffs,
+          cstride, c_offsets, constants, cell_info, get_perm);
     }
   }
 }

cpp/dolfinx/fem/assemble_scalar_impl.h

@@ -23,10 +23,10 @@ namespace dolfinx::fem::impl
 
 /// Assemble functional over cells
 template <typename T>
-T assemble_cells(const mesh::Geometry& geometry,
-                 std::span<const std::int32_t> cells, FEkernel<T> auto fn,
-                 std::span<const T> constants, std::span<const T> coeffs,
-                 int cstride)
+T assemble_cells(
+    const mesh::Geometry<scalar_value_type_t<scalar_value_type_t<T>>>& geometry,
+    std::span<const std::int32_t> cells, FEkernel<T> auto fn,
+    std::span<const T> constants, std::span<const T> coeffs, int cstride)
 {
   T value(0);
   if (cells.empty())
@@ -35,7 +35,7 @@ T assemble_cells(const mesh::Geometry& geometry,
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Create data structures used in assembly
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
@@ -63,10 +63,10 @@ T assemble_cells(const mesh::Geometry& geometry,
 
 /// Execute kernel over exterior facets and accumulate result
 template <typename T>
-T assemble_exterior_facets(const mesh::Geometry& geometry,
-                           std::span<const std::int32_t> facets,
-                           FEkernel<T> auto fn, std::span<const T> constants,
-                           std::span<const T> coeffs, int cstride)
+T assemble_exterior_facets(
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    std::span<const std::int32_t> facets, FEkernel<T> auto fn,
+    std::span<const T> constants, std::span<const T> coeffs, int cstride)
 {
   T value(0);
   if (facets.empty())
@@ -75,7 +75,7 @@ T assemble_exterior_facets(const mesh::Geometry& geometry,
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Create data structures used in assembly
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
@@ -105,12 +105,11 @@ T assemble_exterior_facets(const mesh::Geometry& geometry,
 
 /// Assemble functional over interior facets
 template <typename T>
-T assemble_interior_facets(const mesh::Geometry& geometry, int num_cell_facets,
-                           std::span<const std::int32_t> facets,
-                           FEkernel<T> auto fn, std::span<const T> constants,
-                           std::span<const T> coeffs, int cstride,
-                           std::span<const int> offsets,
-                           std::span<const std::uint8_t> perms)
+T assemble_interior_facets(
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry, int num_cell_facets,
+    std::span<const std::int32_t> facets, FEkernel<T> auto fn,
+    std::span<const T> constants, std::span<const T> coeffs, int cstride,
+    std::span<const int> offsets, std::span<const std::uint8_t> perms)
 {
   T value(0);
   if (facets.empty())
@@ -119,7 +118,7 @@ T assemble_interior_facets(const mesh::Geometry& geometry, int num_cell_facets,
   // 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 = geometry.x();
+  auto x = geometry.x();
 
   // Create data structures used in assembly
   using X = scalar_value_type_t<T>;
@@ -161,10 +160,12 @@ T assemble_interior_facets(const mesh::Geometry& geometry, int num_cell_facets,
   return value;
 }
 
-/// Assemble functional into an scalar
+/// Assemble functional into an scalar with provided mesh geometry.
 template <typename T>
 T assemble_scalar(
-    const fem::Form<T>& M, std::span<const T> constants,
+    const fem::Form<T>& M,
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients)
 {
@@ -177,8 +178,8 @@ T assemble_scalar(
     const auto& fn = M.kernel(IntegralType::cell, i);
     const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
     const std::vector<std::int32_t>& cells = M.cell_domains(i);
-    value += impl::assemble_cells(mesh->geometry(), cells, fn, constants,
-                                  coeffs, cstride);
+    value += impl::assemble_cells(geometry, cells, fn, constants, coeffs,
+                                  cstride);
   }
 
   for (int i : M.integral_ids(IntegralType::exterior_facet))
@@ -187,8 +188,8 @@ T assemble_scalar(
     const auto& [coeffs, cstride]
         = coefficients.at({IntegralType::exterior_facet, i});
     const std::vector<std::int32_t>& facets = M.exterior_facet_domains(i);
-    value += impl::assemble_exterior_facets(mesh->geometry(), facets, fn,
-                                            constants, coeffs, cstride);
+    value += impl::assemble_exterior_facets(geometry, facets, fn, constants,
+                                            coeffs, cstride);
   }
 
   if (M.num_integrals(IntegralType::interior_facet) > 0)
@@ -207,9 +208,9 @@ T assemble_scalar(
       const auto& [coeffs, cstride]
           = coefficients.at({IntegralType::interior_facet, i});
       const std::vector<std::int32_t>& facets = M.interior_facet_domains(i);
-      value += impl::assemble_interior_facets(mesh->geometry(), num_cell_facets,
-                                              facets, fn, constants, coeffs,
-                                              cstride, c_offsets, perms);
+      value += impl::assemble_interior_facets(geometry, num_cell_facets, facets,
+                                              fn, constants, coeffs, cstride,
+                                              c_offsets, perms);
     }
   }