diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
index a4532d1cd43..efe72dd1dae 100644
--- a/cpp/demo/poisson/main.cpp
+++ b/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});
 
diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp
index 004dfee2732..0e1f2a05707 100644
--- a/cpp/demo/poisson_matrix_free/main.cpp
+++ b/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
diff --git a/cpp/dolfinx/fem/DirichletBC.h b/cpp/dolfinx/fem/DirichletBC.h
index 9e16a4f9e29..88f9fc7f200 100644
--- a/cpp/dolfinx/fem/DirichletBC.h
+++ b/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))
     {
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index e5aefd0bf08..b4fe30ceb40 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/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);
     }
   }
 }
diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h
index a968f0986ac..a72db267ead 100644
--- a/cpp/dolfinx/fem/assemble_scalar_impl.h
+++ b/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<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);
     }
   }
 
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index bbf42cdadc9..2af569f7e68 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -37,8 +37,8 @@ namespace dolfinx::fem::impl
 /// @tparam _bs1 The block size of the trial function dof map.
 template <typename T, int _bs0 = -1, int _bs1 = -1>
 void _lift_bc_cells(
-    std::span<T> b, const mesh::Geometry& geometry, FEkernel<T> auto kernel,
-    std::span<const std::int32_t> cells,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    FEkernel<T> auto kernel, std::span<const std::int32_t> cells,
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
@@ -49,8 +49,7 @@ void _lift_bc_cells(
     const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
     std::span<const T> constants, std::span<const T> coeffs, int cstride,
     std::span<const std::uint32_t> cell_info, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
-    double scale)
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T scale)
 {
   assert(_bs0 < 0 or _bs0 == bs0);
   assert(_bs1 < 0 or _bs1 == bs1);
@@ -61,13 +60,11 @@ void _lift_bc_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();
 
   // Data structures used in bc application
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
   std::vector<T> Ae, be;
-  const scalar_value_type_t<T> _scale
-      = static_cast<scalar_value_type_t<T>>(scale);
   for (std::size_t index = 0; index < cells.size(); ++index)
   {
     std::int32_t c = cells[index];
@@ -148,7 +145,7 @@ void _lift_bc_cells(
             // const T _x0 = 0.0;
             // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
             for (int m = 0; m < num_rows; ++m)
-              be[m] -= Ae[m * num_cols + _bs1 * j + k] * _scale * (bc - _x0);
+              be[m] -= Ae[m * num_cols + _bs1 * j + k] * scale * (bc - _x0);
           }
         }
       }
@@ -164,7 +161,7 @@ void _lift_bc_cells(
             const T _x0 = x0.empty() ? 0.0 : x0[jj];
             // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
             for (int m = 0; m < num_rows; ++m)
-              be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
+              be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
           }
         }
       }
@@ -194,8 +191,8 @@ void _lift_bc_cells(
 /// @tparam _bs1 The block size of the trial function dof map.
 template <typename T, int _bs = -1>
 void _lift_bc_exterior_facets(
-    std::span<T> b, const mesh::Geometry& geometry, FEkernel<T> auto kernel,
-    std::span<const std::int32_t> facets,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    FEkernel<T> auto kernel, 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,
@@ -206,8 +203,7 @@ void _lift_bc_exterior_facets(
     const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
     std::span<const T> constants, std::span<const T> coeffs, int cstride,
     std::span<const std::uint32_t> cell_info, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
-    double scale)
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T scale)
 {
   if (facets.empty())
     return;
@@ -215,14 +211,12 @@ void _lift_bc_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 bc application
   std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
   std::vector<T> Ae, be;
   assert(facets.size() % 2 == 0);
-  const scalar_value_type_t<T> _scale
-      = static_cast<scalar_value_type_t<T>>(scale);
   for (std::size_t index = 0; index < facets.size(); index += 2)
   {
     std::int32_t cell = facets[index];
@@ -283,7 +277,7 @@ void _lift_bc_exterior_facets(
           const T _x0 = x0.empty() ? 0.0 : x0[jj];
           // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
           for (int m = 0; m < num_rows; ++m)
-            be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
+            be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
         }
       }
     }
@@ -302,8 +296,9 @@ void _lift_bc_exterior_facets(
 /// @tparam _bs1 The block size of the trial function dof map.
 template <typename T, int _bs = -1>
 void _lift_bc_interior_facets(
-    std::span<T> b, const mesh::Geometry& geometry, int num_cell_facets,
-    FEkernel<T> auto kernel, std::span<const std::int32_t> facets,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    int num_cell_facets, FEkernel<T> auto kernel,
+    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,
@@ -316,7 +311,7 @@ void _lift_bc_interior_facets(
     std::span<const std::uint32_t> cell_info,
     const std::function<std::uint8_t(std::size_t)>& get_perm,
     std::span<const T> bc_values1, std::span<const std::int8_t> bc_markers1,
-    std::span<const T> x0, double scale)
+    std::span<const T> x0, T scale)
 {
   if (facets.empty())
     return;
@@ -324,7 +319,7 @@ void _lift_bc_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>;
@@ -337,8 +332,6 @@ void _lift_bc_interior_facets(
   std::vector<std::int32_t> dmapjoint0, dmapjoint1;
   assert(facets.size() % 4 == 0);
 
-  const scalar_value_type_t<T> _scale
-      = static_cast<scalar_value_type_t<T>>(scale);
   for (std::size_t index = 0; index < facets.size(); index += 4)
   {
     std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
@@ -443,7 +436,7 @@ void _lift_bc_interior_facets(
           const T _x0 = x0.empty() ? 0.0 : x0[jj];
           // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
           for (int m = 0; m < num_rows; ++m)
-            be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0);
+            be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
         }
       }
     }
@@ -462,8 +455,8 @@ void _lift_bc_interior_facets(
           // be -= Ae.col(offset + bs1 * j + k) * scale * (bc - x0[jj]);
           for (int m = 0; m < num_rows; ++m)
           {
-            be[m] -= Ae[m * num_cols + offset + bs1 * j + k] * _scale
-                     * (bc - _x0);
+            be[m]
+                -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0);
           }
         }
       }
@@ -490,7 +483,7 @@ void assemble_cells(
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
-    std::span<T> b, const mesh::Geometry& geometry,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
     std::span<const std::int32_t> cells,
     const graph::AdjacencyList<std::int32_t>& dofmap, int bs,
     FEkernel<T> auto kernel, std::span<const T> constants,
@@ -505,7 +498,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();
 
   // FIXME: Add proper interface for num_dofs
   // Create data structures used in assembly
@@ -561,7 +554,7 @@ void assemble_exterior_facets(
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
-    std::span<T> b, const mesh::Geometry& geometry,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
     std::span<const std::int32_t> facets,
     const graph::AdjacencyList<std::int32_t>& dofmap, int bs,
     FEkernel<T> auto fn, std::span<const T> constants,
@@ -576,7 +569,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();
 
   // FIXME: Add proper interface for num_dofs
   // Create data structures used in assembly
@@ -633,17 +626,17 @@ void assemble_interior_facets(
     const std::function<void(const std::span<T>&,
                              const std::span<const std::uint32_t>&,
                              std::int32_t, int)>& dof_transform,
-    std::span<T> b, const mesh::Geometry& geometry, int num_cell_facets,
-    std::span<const std::int32_t> facets, const fem::DofMap& dofmap,
-    FEkernel<T> auto fn, std::span<const T> constants,
-    std::span<const T> coeffs, int cstride,
+    std::span<T> b, const mesh::Geometry<scalar_value_type_t<T>>& geometry,
+    int num_cell_facets, std::span<const std::int32_t> facets,
+    const fem::DofMap& dofmap, FEkernel<T> auto fn,
+    std::span<const T> constants, std::span<const T> coeffs, int cstride,
     std::span<const std::uint32_t> cell_info,
     const std::function<std::uint8_t(std::size_t)>& get_perm)
 {
   // 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>;
@@ -722,6 +715,7 @@ void assemble_interior_facets(
 ///
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear form that generates A
+/// @param[in] geometry The mesh geometry
 /// @param[in] constants Constants that appear in `a`
 /// @param[in] coefficients Coefficients that appear in `a`
 /// @param[in] bc_values1 The boundary condition 'values'
@@ -731,12 +725,14 @@ void assemble_interior_facets(
 /// solution' in a Newton method
 /// @param[in] scale Scaling to apply
 template <typename T>
-void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
+void lift_bc(std::span<T> b, 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 T> bc_values1,
              std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
-             double scale)
+             T scale)
 {
   std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
   assert(mesh);
@@ -784,22 +780,22 @@ void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
     const std::vector<std::int32_t>& cells = a.cell_domains(i);
     if (bs0 == 1 and bs1 == 1)
     {
-      _lift_bc_cells<T, 1, 1>(b, mesh->geometry(), kernel, cells, dof_transform,
+      _lift_bc_cells<T, 1, 1>(b, geometry, kernel, cells, dof_transform,
                               dofmap0, bs0, dof_transform_to_transpose, dofmap1,
                               bs1, constants, coeffs, cstride, cell_info,
                               bc_values1, bc_markers1, x0, scale);
     }
     else if (bs0 == 3 and bs1 == 3)
     {
-      _lift_bc_cells<T, 3, 3>(b, mesh->geometry(), kernel, cells, dof_transform,
+      _lift_bc_cells<T, 3, 3>(b, geometry, kernel, cells, dof_transform,
                               dofmap0, bs0, dof_transform_to_transpose, dofmap1,
                               bs1, constants, coeffs, cstride, cell_info,
                               bc_values1, bc_markers1, x0, scale);
     }
     else
     {
-      _lift_bc_cells(b, mesh->geometry(), kernel, cells, dof_transform, dofmap0,
-                     bs0, dof_transform_to_transpose, dofmap1, bs1, constants,
+      _lift_bc_cells(b, geometry, kernel, cells, dof_transform, dofmap0, bs0,
+                     dof_transform_to_transpose, dofmap1, bs1, constants,
                      coeffs, cstride, cell_info, bc_values1, bc_markers1, x0,
                      scale);
     }
@@ -811,7 +807,7 @@ void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
     const auto& [coeffs, cstride]
         = coefficients.at({IntegralType::exterior_facet, i});
     const std::vector<std::int32_t>& facets = a.exterior_facet_domains(i);
-    _lift_bc_exterior_facets(b, mesh->geometry(), kernel, facets, dof_transform,
+    _lift_bc_exterior_facets(b, geometry, kernel, facets, dof_transform,
                              dofmap0, bs0, dof_transform_to_transpose, dofmap1,
                              bs1, constants, coeffs, cstride, cell_info,
                              bc_values1, bc_markers1, x0, scale);
@@ -839,11 +835,10 @@ void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
       const auto& [coeffs, cstride]
           = coefficients.at({IntegralType::interior_facet, i});
       const std::vector<std::int32_t>& facets = a.interior_facet_domains(i);
-      _lift_bc_interior_facets(b, mesh->geometry(), num_cell_facets, kernel,
-                               facets, dof_transform, dofmap0, bs0,
-                               dof_transform_to_transpose, dofmap1, bs1,
-                               constants, coeffs, cstride, cell_info, get_perm,
-                               bc_values1, bc_markers1, x0, scale);
+      _lift_bc_interior_facets(
+          b, geometry, num_cell_facets, kernel, facets, dof_transform, dofmap0,
+          bs0, dof_transform_to_transpose, dofmap1, bs1, constants, coeffs,
+          cstride, cell_info, get_perm, bc_values1, bc_markers1, x0, scale);
     }
   }
 }
@@ -860,6 +855,7 @@ void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear forms, where a[j] is the form that
 /// generates A_j
+/// @param[in] geometry The mesh geometry
 /// @param[in] constants Constants that appear in `a`
 /// @param[in] coeffs Coefficients that appear in `a`
 /// @param[in] bcs1 List of boundary conditions for each block, i.e.
@@ -869,12 +865,13 @@ void lift_bc(std::span<T> b, const Form<T>& a, std::span<const T> constants,
 /// @param[in] scale Scaling to apply
 template <typename T>
 void apply_lifting(
-    std::span<T> b, std::vector<std::shared_ptr<const Form<T>>> a,
+    std::span<T> b, const std::vector<std::shared_ptr<const Form<T>>> a,
+    const mesh::Geometry<scalar_value_type_t<T>>& geometry,
     const std::vector<std::span<const T>>& constants,
     const std::vector<std::map<std::pair<IntegralType, int>,
                                std::pair<std::span<const T>, int>>>& coeffs,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
-    const std::vector<std::span<const T>>& x0, double scale)
+    const std::vector<std::span<const T>>& x0, T scale)
 {
   // FIXME: make changes to reactivate this check
   if (!x0.empty() and x0.size() != a.size())
@@ -913,13 +910,13 @@ void apply_lifting(
 
       if (!x0.empty())
       {
-        lift_bc<T>(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1,
-                   x0[j], scale);
+        lift_bc<T>(b, *a[j], geometry, constants[j], coeffs[j], bc_values1,
+                   bc_markers1, x0[j], scale);
       }
       else
       {
-        lift_bc<T>(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1,
-                   std::span<const T>(), scale);
+        lift_bc<T>(b, *a[j], geometry, constants[j], coeffs[j], bc_values1,
+                   bc_markers1, std::span<const T>(), scale);
       }
     }
   }
@@ -929,11 +926,14 @@ void apply_lifting(
 /// @param[in,out] b The vector to be assembled. It will not be zeroed
 /// before assembly.
 /// @param[in] L The linear forms to assemble into b
+/// @param[in] geometry The mesh geometry
 /// @param[in] constants Packed constants that appear in `L`
 /// @param[in] coefficients Packed coefficients that appear in `L`
 template <typename T>
 void assemble_vector(
-    std::span<T> b, const Form<T>& L, std::span<const T> constants,
+    std::span<T> b, const Form<T>& L,
+    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)
 {
@@ -971,20 +971,18 @@ void assemble_vector(
     const std::vector<std::int32_t>& cells = L.cell_domains(i);
     if (bs == 1)
     {
-      impl::assemble_cells<T, 1>(dof_transform, b, mesh->geometry(), cells,
-                                 dofs, bs, fn, constants, coeffs, cstride,
-                                 cell_info);
+      impl::assemble_cells<T, 1>(dof_transform, b, geometry, cells, dofs, bs,
+                                 fn, constants, coeffs, cstride, cell_info);
     }
     else if (bs == 3)
     {
-      impl::assemble_cells<T, 3>(dof_transform, b, mesh->geometry(), cells,
-                                 dofs, bs, fn, constants, coeffs, cstride,
-                                 cell_info);
+      impl::assemble_cells<T, 3>(dof_transform, b, geometry, cells, dofs, bs,
+                                 fn, constants, coeffs, cstride, cell_info);
     }
     else
     {
-      impl::assemble_cells(dof_transform, b, mesh->geometry(), cells, dofs, bs,
-                           fn, constants, coeffs, cstride, cell_info);
+      impl::assemble_cells(dof_transform, b, geometry, cells, dofs, bs, fn,
+                           constants, coeffs, cstride, cell_info);
     }
   }
 
@@ -996,20 +994,20 @@ void assemble_vector(
     const std::vector<std::int32_t>& facets = L.exterior_facet_domains(i);
     if (bs == 1)
     {
-      impl::assemble_exterior_facets<T, 1>(dof_transform, b, mesh->geometry(),
-                                           facets, dofs, bs, fn, constants,
-                                           coeffs, cstride, cell_info);
+      impl::assemble_exterior_facets<T, 1>(dof_transform, b, geometry, facets,
+                                           dofs, bs, fn, constants, coeffs,
+                                           cstride, cell_info);
     }
     else if (bs == 3)
     {
-      impl::assemble_exterior_facets<T, 3>(dof_transform, b, mesh->geometry(),
-                                           facets, dofs, bs, fn, constants,
-                                           coeffs, cstride, cell_info);
+      impl::assemble_exterior_facets<T, 3>(dof_transform, b, geometry, facets,
+                                           dofs, bs, fn, constants, coeffs,
+                                           cstride, cell_info);
     }
     else
     {
-      impl::assemble_exterior_facets(dof_transform, b, mesh->geometry(), facets,
-                                     dofs, bs, fn, constants, coeffs, cstride,
+      impl::assemble_exterior_facets(dof_transform, b, geometry, facets, dofs,
+                                     bs, fn, constants, coeffs, cstride,
                                      cell_info);
     }
   }
@@ -1038,22 +1036,45 @@ void assemble_vector(
       if (bs == 1)
       {
         impl::assemble_interior_facets<T, 1>(
-            dof_transform, b, mesh->geometry(), num_cell_facets, facets,
-            *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm);
+            dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn,
+            constants, coeffs, cstride, cell_info, get_perm);
       }
       else if (bs == 3)
       {
         impl::assemble_interior_facets<T, 3>(
-            dof_transform, b, mesh->geometry(), num_cell_facets, facets,
-            *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm);
+            dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn,
+            constants, coeffs, cstride, cell_info, get_perm);
       }
       else
       {
         impl::assemble_interior_facets(
-            dof_transform, b, mesh->geometry(), num_cell_facets, facets,
-            *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm);
+            dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn,
+            constants, coeffs, cstride, cell_info, get_perm);
       }
     }
   }
 }
+
+/// Assemble linear form into a vector
+/// @param[in,out] b The vector to be assembled. It will not be zeroed
+/// before assembly.
+/// @param[in] L The linear forms to assemble into b
+/// @param[in] constants Packed constants that appear in `L`
+/// @param[in] coefficients Packed coefficients that appear in `L`
+template <typename T>
+void assemble_vector(
+    std::span<T> b, const Form<T>& L, std::span<const T> constants,
+    const std::map<std::pair<IntegralType, int>,
+                   std::pair<std::span<const T>, int>>& coefficients)
+{
+  std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
+  assert(mesh);
+  if constexpr (std::is_same_v<double, scalar_value_type_t<T>>)
+    assemble_vector(b, L, mesh->geometry(), constants, coefficients);
+  else
+  {
+    assemble_vector(b, L, mesh->geometry().astype<scalar_value_type_t<T>>(),
+                    constants, coefficients);
+  }
+}
 } // namespace dolfinx::fem::impl
diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
index 4500f5e0c8e..440bdaa3388 100644
--- a/cpp/dolfinx/fem/assembler.h
+++ b/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 "utils.h"
 #include <cstdint>
 #include <memory>
 #include <span>
@@ -58,7 +59,16 @@ T assemble_scalar(
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients)
 {
-  return impl::assemble_scalar(M, constants, coefficients);
+  std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
+  assert(mesh);
+  if constexpr (std::is_same_v<double, impl::scalar_value_type_t<T>>)
+    return impl::assemble_scalar(M, mesh->geometry(), constants, coefficients);
+  else
+  {
+    return impl::assemble_scalar(
+        M, mesh->geometry().astype<impl::scalar_value_type_t<T>>(), constants,
+        coefficients);
+  }
 }
 
 /// Assemble functional into scalar
@@ -136,9 +146,21 @@ void apply_lifting(
     const std::vector<std::map<std::pair<IntegralType, int>,
                                std::pair<std::span<const T>, int>>>& coeffs,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
-    const std::vector<std::span<const T>>& x0, double scale)
+    const std::vector<std::span<const T>>& x0, T scale)
 {
-  impl::apply_lifting(b, a, constants, coeffs, bcs1, x0, scale);
+  std::shared_ptr<const mesh::Mesh> mesh = a[0]->mesh();
+  assert(mesh);
+  if constexpr (std::is_same_v<double, impl::scalar_value_type_t<T>>)
+  {
+    impl::apply_lifting<T>(b, a, mesh->geometry(), constants, coeffs, bcs1, x0,
+                           scale);
+  }
+  else
+  {
+    impl::apply_lifting<T>(
+        b, a, mesh->geometry().astype<impl::scalar_value_type_t<T>>(),
+        constants, coeffs, bcs1, x0, scale);
+  }
 }
 
 /// Modify b such that:
@@ -157,7 +179,7 @@ template <typename T>
 void apply_lifting(
     std::span<T> b, const std::vector<std::shared_ptr<const Form<T>>>& a,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
-    const std::vector<std::span<const T>>& x0, double scale)
+    const std::vector<std::span<const T>>& x0, T scale)
 {
   std::vector<
       std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
@@ -186,11 +208,48 @@ void apply_lifting(
       _coeffs;
   std::transform(coeffs.cbegin(), coeffs.cend(), std::back_inserter(_coeffs),
                  [](auto& c) { return make_coefficients_span(c); });
+
   apply_lifting(b, a, _constants, _coeffs, bcs1, x0, scale);
 }
 
 // -- Matrices ---------------------------------------------------------------
 
+/// @brief Assemble bilinear form into a matrix. Matrix must already be
+/// initialised. Does not zero or finalise the matrix.
+/// @param[in] mat_add The function for adding values into the matrix
+/// @param[in] a The bilinear form to assemble
+/// @param[in] constants Constants that appear in `a`
+/// @param[in] coefficients Coefficients that appear in `a`
+/// @param[in] dof_marker0 Boundary condition markers for the rows. If
+/// bc[i] is true then rows i in A will be zeroed. The index i is a
+/// local index.
+/// @param[in] dof_marker1 Boundary condition markers for the columns.
+/// If bc[i] is true then rows i in A will be zeroed. The index i is a
+/// local index.
+template <typename T>
+void assemble_matrix(
+    auto mat_add, 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> dof_marker0,
+    std::span<const std::int8_t> dof_marker1)
+
+{
+  std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
+  assert(mesh);
+  if constexpr (std::is_same_v<double, impl::scalar_value_type_t<T>>)
+  {
+    impl::assemble_matrix(mat_add, a, mesh->geometry(), constants, coefficients,
+                          dof_marker0, dof_marker1);
+  }
+  else
+  {
+    impl::assemble_matrix(
+        mat_add, a, mesh->geometry().astype<impl::scalar_value_type_t<T>>(),
+        constants, coefficients, dof_marker0, dof_marker1);
+  }
+}
+
 /// Assemble bilinear form into a matrix
 /// @param[in] mat_add The function for adding values into the matrix
 /// @param[in] a The bilinear from to assemble
@@ -235,8 +294,8 @@ void assemble_matrix(
   }
 
   // Assemble
-  impl::assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
-                        dof_marker1);
+  assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
+                  dof_marker1);
 }
 
 /// Assemble bilinear form into a matrix
@@ -259,31 +318,6 @@ void assemble_matrix(
                   make_coefficients_span(coefficients), bcs);
 }
 
-/// @brief Assemble bilinear form into a matrix. Matrix must already be
-/// initialised. Does not zero or finalise the matrix.
-/// @param[in] mat_add The function for adding values into the matrix
-/// @param[in] a The bilinear form to assemble
-/// @param[in] constants Constants that appear in `a`
-/// @param[in] coefficients Coefficients that appear in `a`
-/// @param[in] dof_marker0 Boundary condition markers for the rows. If
-/// bc[i] is true then rows i in A will be zeroed. The index i is a
-/// local index.
-/// @param[in] dof_marker1 Boundary condition markers for the columns.
-/// If bc[i] is true then rows i in A will be zeroed. The index i is a
-/// local index.
-template <typename T>
-void assemble_matrix(
-    auto mat_add, 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> dof_marker0,
-    std::span<const std::int8_t> dof_marker1)
-
-{
-  impl::assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
-                        dof_marker1);
-}
-
 /// @brief Assemble bilinear form into a matrix. Matrix must already be
 /// initialised. Does not zero or finalise the matrix.
 /// @param[in] mat_add The function for adding values into the matrix
@@ -378,7 +412,7 @@ void set_diagonal(auto set_fn, const fem::FunctionSpace& V,
 template <typename T>
 void set_bc(std::span<T> b,
             const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
-            std::span<const T> x0, double scale = 1.0)
+            std::span<const T> x0, T scale = 1)
 {
   if (b.size() > x0.size())
     throw std::runtime_error("Size mismatch between b and x0 vectors.");
@@ -395,7 +429,7 @@ void set_bc(std::span<T> b,
 template <typename T>
 void set_bc(std::span<T> b,
             const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
-            double scale = 1.0)
+            T scale = 1)
 {
   for (const auto& bc : bcs)
   {
diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp
index 6cf4d0a986a..74426c04820 100644
--- a/cpp/dolfinx/fem/petsc.cpp
+++ b/cpp/dolfinx/fem/petsc.cpp
@@ -304,7 +304,7 @@ void fem::petsc::apply_lifting(
         coeffs,
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>>& bcs1,
-    const std::vector<Vec>& x0, double scale)
+    const std::vector<Vec>& x0, PetscScalar scale)
 {
   Vec b_local;
   VecGhostGetLocalForm(b, &b_local);
@@ -350,7 +350,7 @@ void fem::petsc::apply_lifting(
     Vec b, const std::vector<std::shared_ptr<const Form<PetscScalar>>>& a,
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>>& bcs1,
-    const std::vector<Vec>& x0, double scale)
+    const std::vector<Vec>& x0, PetscScalar scale)
 {
   Vec b_local;
   VecGhostGetLocalForm(b, &b_local);
@@ -394,7 +394,7 @@ void fem::petsc::apply_lifting(
 void fem::petsc::set_bc(
     Vec b,
     const std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>& bcs,
-    const Vec x0, double scale)
+    const Vec x0, PetscScalar scale)
 {
   PetscInt n = 0;
   VecGetLocalSize(b, &n);
diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h
index f1b59faf8a0..341355756c5 100644
--- a/cpp/dolfinx/fem/petsc.h
+++ b/cpp/dolfinx/fem/petsc.h
@@ -128,7 +128,7 @@ void apply_lifting(
         coeffs,
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>>& bcs1,
-    const std::vector<Vec>& x0, double scale);
+    const std::vector<Vec>& x0, PetscScalar scale);
 
 // FIXME: clarify how x0 is used
 // FIXME: if bcs entries are set
@@ -152,7 +152,7 @@ void apply_lifting(
     Vec b, const std::vector<std::shared_ptr<const Form<PetscScalar>>>& a,
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>>& bcs1,
-    const std::vector<Vec>& x0, double scale);
+    const std::vector<Vec>& x0, PetscScalar scale);
 
 // -- Setting bcs ------------------------------------------------------------
 
@@ -167,6 +167,6 @@ void apply_lifting(
 void set_bc(
     Vec b,
     const std::vector<std::shared_ptr<const DirichletBC<PetscScalar>>>& bcs,
-    const Vec x0, double scale = 1.0);
+    const Vec x0, PetscScalar scale = 1);
 } // namespace petsc
 } // namespace dolfinx::fem
diff --git a/cpp/dolfinx/geometry/utils.cpp b/cpp/dolfinx/geometry/utils.cpp
index 14c70976999..b40c15610c8 100644
--- a/cpp/dolfinx/geometry/utils.cpp
+++ b/cpp/dolfinx/geometry/utils.cpp
@@ -386,7 +386,7 @@ geometry::shortest_vector(const mesh::Mesh& mesh, int dim,
                           std::span<const double> points)
 {
   const int tdim = mesh.topology().dim();
-  const mesh::Geometry& geometry = mesh.geometry();
+  const mesh::Geometry<double>& geometry = mesh.geometry();
   std::span<const double> geom_dofs = geometry.x();
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
   std::vector<double> shortest_vectors(3 * entities.size());
@@ -509,7 +509,7 @@ int geometry::compute_first_colliding_cell(
   else
   {
     constexpr double eps2 = 1e-20;
-    const mesh::Geometry& geometry = mesh.geometry();
+    const mesh::Geometry<double>& geometry = mesh.geometry();
     std::span<const double> geom_dofs = geometry.x();
     const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
     const std::size_t num_nodes = geometry.cmap().dim();
diff --git a/cpp/dolfinx/io/ADIOS2Writers.cpp b/cpp/dolfinx/io/ADIOS2Writers.cpp
index ec231fc4796..deea238cf27 100644
--- a/cpp/dolfinx/io/ADIOS2Writers.cpp
+++ b/cpp/dolfinx/io/ADIOS2Writers.cpp
@@ -145,7 +145,7 @@ void vtx_write_data(adios2::IO& io, adios2::Engine& engine,
 void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine,
                     const mesh::Mesh& mesh)
 {
-  const mesh::Geometry& geometry = mesh.geometry();
+  const mesh::Geometry<double>& geometry = mesh.geometry();
   const mesh::Topology& topology = mesh.topology();
 
   // "Put" geometry
@@ -411,7 +411,7 @@ std::vector<T> pack_function_data(const fem::Function<T>& u)
   assert(dofmap);
   auto mesh = V->mesh();
   assert(mesh);
-  const mesh::Geometry& geometry = mesh->geometry();
+  const mesh::Geometry<double>& geometry = mesh->geometry();
   const mesh::Topology& topology = mesh->topology();
 
   // The Function and the mesh must have identical element_dof_layouts
@@ -534,7 +534,7 @@ void fides_write_data(adios2::IO& io, adios2::Engine& engine,
 void fides_write_mesh(adios2::IO& io, adios2::Engine& engine,
                       const mesh::Mesh& mesh)
 {
-  const mesh::Geometry& geometry = mesh.geometry();
+  const mesh::Geometry<double>& geometry = mesh.geometry();
   const mesh::Topology& topology = mesh.topology();
 
   // "Put" geometry data
@@ -568,7 +568,7 @@ void fides_write_mesh(adios2::IO& io, adios2::Engine& engine,
 /// @param[in] mesh The mesh
 void fides_initialize_mesh_attributes(adios2::IO& io, const mesh::Mesh& mesh)
 {
-  const mesh::Geometry& geometry = mesh.geometry();
+  const mesh::Geometry<double>& geometry = mesh.geometry();
   const mesh::Topology& topology = mesh.topology();
 
   // Check that mesh is first order mesh
diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp
index e1a9fe4249c..d8d666aebef 100644
--- a/cpp/dolfinx/io/VTKFile.cpp
+++ b/cpp/dolfinx/io/VTKFile.cpp
@@ -421,7 +421,7 @@ void write_function(
     std::vector<std::int64_t> tmp;
     std::tie(tmp, cshape) = io::extract_vtk_connectivity(*mesh0);
     cells.assign(tmp.begin(), tmp.end());
-    const mesh::Geometry& geometry = mesh0->geometry();
+    const mesh::Geometry<double>& geometry = mesh0->geometry();
     x.assign(geometry.x().begin(), geometry.x().end());
     xshape = {geometry.x().size() / 3, 3};
     x_id = geometry.input_global_indices();
@@ -749,7 +749,7 @@ void io::VTKFile::write(const mesh::Mesh& mesh, double time)
 
   // Get mesh data for this rank
   const mesh::Topology& topology = mesh.topology();
-  const mesh::Geometry& geometry = mesh.geometry();
+  const mesh::Geometry<double>& geometry = mesh.geometry();
   auto xmap = geometry.index_map();
   assert(xmap);
   const int tdim = topology.dim();
diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp
index d0bb5fd8292..94f8a0f884a 100644
--- a/cpp/dolfinx/io/XDMFFile.cpp
+++ b/cpp/dolfinx/io/XDMFFile.cpp
@@ -202,7 +202,7 @@ void XDMFFile::write_mesh(const mesh::Mesh& mesh, std::string xpath)
     _xml_doc->save_file(_filename.c_str(), "  ");
 }
 //-----------------------------------------------------------------------------
-void XDMFFile::write_geometry(const mesh::Geometry& geometry, std::string name,
+void XDMFFile::write_geometry(const mesh::Geometry<double>& geometry, std::string name,
                               std::string xpath)
 {
   pugi::xml_node node = _xml_doc->select_node(xpath.c_str()).node();
diff --git a/cpp/dolfinx/io/XDMFFile.h b/cpp/dolfinx/io/XDMFFile.h
index fd6f07714f0..5fab785ca5d 100644
--- a/cpp/dolfinx/io/XDMFFile.h
+++ b/cpp/dolfinx/io/XDMFFile.h
@@ -32,6 +32,7 @@ class Function;
 
 namespace dolfinx::mesh
 {
+template <typename T>
 class Geometry;
 enum class GhostMode : int;
 class Mesh;
@@ -92,7 +93,7 @@ class XDMFFile
   /// @param[in] geometry
   /// @param[in] name
   /// @param[in] xpath XPath of a node where Geometry will be inserted
-  void write_geometry(const mesh::Geometry& geometry, std::string name,
+  void write_geometry(const mesh::Geometry<double>& geometry, std::string name,
                       std::string xpath = "/Xdmf/Domain");
 
   /// Read in Mesh
diff --git a/cpp/dolfinx/io/xdmf_mesh.cpp b/cpp/dolfinx/io/xdmf_mesh.cpp
index 66663ab611a..2f9cf861491 100644
--- a/cpp/dolfinx/io/xdmf_mesh.cpp
+++ b/cpp/dolfinx/io/xdmf_mesh.cpp
@@ -21,7 +21,7 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node,
                                   const hid_t h5_id,
                                   const std::string path_prefix,
                                   const mesh::Topology& topology,
-                                  const mesh::Geometry& geometry, int dim,
+                                  const mesh::Geometry<double>& geometry, int dim,
                                   std::span<const std::int32_t> entities)
 {
   LOG(INFO) << "Adding topology data to node \"" << xml_node.path('/') << "\"";
@@ -149,7 +149,7 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node,
 void xdmf_mesh::add_geometry_data(MPI_Comm comm, pugi::xml_node& xml_node,
                                   const hid_t h5_id,
                                   const std::string path_prefix,
-                                  const mesh::Geometry& geometry)
+                                  const mesh::Geometry<double>& geometry)
 {
   LOG(INFO) << "Adding geometry data to node \"" << xml_node.path('/') << "\"";
   auto map = geometry.index_map();
diff --git a/cpp/dolfinx/io/xdmf_mesh.h b/cpp/dolfinx/io/xdmf_mesh.h
index 2e4f150013b..1c6a4391168 100644
--- a/cpp/dolfinx/io/xdmf_mesh.h
+++ b/cpp/dolfinx/io/xdmf_mesh.h
@@ -24,6 +24,7 @@ namespace dolfinx
 
 namespace mesh
 {
+template <typename T>
 class Geometry;
 class Mesh;
 class Topology;
@@ -53,13 +54,13 @@ void add_mesh(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id,
 void add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node,
                        const hid_t h5_id, const std::string path_prefix,
                        const mesh::Topology& topology,
-                       const mesh::Geometry& geometry, int cell_dim,
+                       const mesh::Geometry<double>& geometry, int cell_dim,
                        std::span<const std::int32_t> entities);
 
 /// Add Geometry xml node
 void add_geometry_data(MPI_Comm comm, pugi::xml_node& xml_node,
                        const hid_t h5_id, const std::string path_prefix,
-                       const mesh::Geometry& geometry);
+                       const mesh::Geometry<double>& geometry);
 
 /// @brief Read geometry (coordinate) data.
 ///
diff --git a/cpp/dolfinx/mesh/Geometry.cpp b/cpp/dolfinx/mesh/Geometry.cpp
index f19731a4095..d87ff24ac2b 100644
--- a/cpp/dolfinx/mesh/Geometry.cpp
+++ b/cpp/dolfinx/mesh/Geometry.cpp
@@ -16,32 +16,7 @@ using namespace dolfinx;
 using namespace dolfinx::mesh;
 
 //-----------------------------------------------------------------------------
-int Geometry::dim() const { return _dim; }
-//-----------------------------------------------------------------------------
-const graph::AdjacencyList<std::int32_t>& Geometry::dofmap() const
-{
-  return _dofmap;
-}
-//-----------------------------------------------------------------------------
-std::shared_ptr<const common::IndexMap> Geometry::index_map() const
-{
-  return _index_map;
-}
-//-----------------------------------------------------------------------------
-std::span<double> Geometry::x() { return _x; }
-//-----------------------------------------------------------------------------
-std::span<const double> Geometry::x() const { return _x; }
-//-----------------------------------------------------------------------------
-const fem::CoordinateElement& Geometry::cmap() const { return _cmap; }
-//-----------------------------------------------------------------------------
-const std::vector<std::int64_t>& Geometry::input_global_indices() const
-{
-  return _input_global_indices;
-}
-//-----------------------------------------------------------------------------
-
-//-----------------------------------------------------------------------------
-mesh::Geometry mesh::create_geometry(
+mesh::Geometry<double> mesh::create_geometry(
     MPI_Comm comm, const Topology& topology,
     const fem::CoordinateElement& element,
     const graph::AdjacencyList<std::int64_t>& cell_nodes,
@@ -116,13 +91,13 @@ mesh::Geometry mesh::create_geometry(
                 std::next(xg.begin(), 3 * i));
   }
 
-  return Geometry(dof_index_map, std::move(dofmap), element, std::move(xg), dim,
-                  std::move(igi));
+  return Geometry<double>(dof_index_map, std::move(dofmap), element,
+                          std::move(xg), dim, std::move(igi));
 }
 //-----------------------------------------------------------------------------
-std::pair<mesh::Geometry, std::vector<int32_t>>
-mesh::create_subgeometry(const Topology& topology, const Geometry& geometry,
-                         int dim,
+std::pair<mesh::Geometry<double>, std::vector<int32_t>>
+mesh::create_subgeometry(const Topology& topology,
+                         const Geometry<double>& geometry, int dim,
                          std::span<const std::int32_t> subentity_to_entity)
 {
   // Get the geometry dofs in the sub-geometry based on the entities in
@@ -237,8 +212,9 @@ mesh::create_subgeometry(const Topology& topology, const Geometry& geometry,
                  [&igi](std::int32_t sub_x_dof) { return igi[sub_x_dof]; });
 
   // Create geometry
-  return {Geometry(sub_x_dof_index_map, std::move(sub_x_dofmap), sub_coord_ele,
-                   std::move(sub_x), geometry.dim(), std::move(sub_igi)),
+  return {Geometry<double>(sub_x_dof_index_map, std::move(sub_x_dofmap),
+                           sub_coord_ele, std::move(sub_x), geometry.dim(),
+                           std::move(sub_igi)),
           std::move(subx_to_x_dofmap)};
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h
index 59120ecb6d1..eba9202eaaf 100644
--- a/cpp/dolfinx/mesh/Geometry.h
+++ b/cpp/dolfinx/mesh/Geometry.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2006-2020 Anders Logg and Garth N. Wells
+// Copyright (C) 2006-2022 Anders Logg and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -26,6 +26,7 @@ namespace dolfinx::mesh
 class Topology;
 
 /// @brief Geometry stores the geometry imposed on a mesh.
+template <typename T>
 class Geometry
 {
 public:
@@ -42,7 +43,7 @@ class Geometry
   /// point, commonly from a mesh input file. The type is
   /// `std:vector<std::int64_t>`.
   template <std::convertible_to<graph::AdjacencyList<std::int32_t>> U,
-            std::convertible_to<std::vector<double>> V,
+            std::convertible_to<std::vector<T>> V,
             std::convertible_to<std::vector<std::int64_t>> W>
   Geometry(std::shared_ptr<const common::IndexMap> index_map, U&& dofmap,
            const fem::CoordinateElement& element, V&& x, int dim,
@@ -71,35 +72,50 @@ class Geometry
   /// Move Assignment
   Geometry& operator=(Geometry&&) = default;
 
+  /// Copy constructor
+  template <typename U>
+  Geometry<U> astype() const
+  {
+    return Geometry<U>(_index_map, _dofmap, _cmap,
+                       std::vector<U>(_x.begin(), _x.end()), _dim,
+                       _input_global_indices);
+  }
+
   /// Return Euclidean dimension of coordinate system
-  int dim() const;
+  int dim() const { return _dim; }
 
   /// DOF map
-  const graph::AdjacencyList<std::int32_t>& dofmap() const;
+  const graph::AdjacencyList<std::int32_t>& dofmap() const { return _dofmap; }
 
   /// Index map
-  std::shared_ptr<const common::IndexMap> index_map() const;
+  std::shared_ptr<const common::IndexMap> index_map() const
+  {
+    return _index_map;
+  }
 
   /// @brief Access geometry degrees-of-freedom data (const version).
   ///
   /// @return The flattened row-major geometry data, where the shape is
   /// (num_points, 3)
-  std::span<const double> x() const;
+  std::span<const T> x() const { return _x; }
 
   /// @brief Access geometry degrees-of-freedom data (non-const
   /// version).
   ///
   /// @return The flattened row-major geometry data, where the shape is
   /// (num_points, 3)
-  std::span<double> x();
+  std::span<T> x() { return _x; }
 
   /// @brief The element that describes the geometry map.
   ///
   /// @return The coordinate/geometry element
-  const fem::CoordinateElement& cmap() const;
+  const fem::CoordinateElement& cmap() const { return _cmap; }
 
   /// Global user indices
-  const std::vector<std::int64_t>& input_global_indices() const;
+  const std::vector<std::int64_t>& input_global_indices() const
+  {
+    return _input_global_indices;
+  }
 
 private:
   // Geometric dimension
@@ -116,7 +132,7 @@ class Geometry
 
   // Coordinates for all points stored as a contiguous array (row-major,
   // column size = 3)
-  std::vector<double> _x;
+  std::vector<T> _x;
 
   // Global indices as provided on Geometry creation
   std::vector<std::int64_t> _input_global_indices;
@@ -125,8 +141,8 @@ class Geometry
 /// @brief Build Geometry from input data.
 ///
 /// This function should be called after the mesh topology is built. It
-/// distributes the 'node' coordinate data to the required MPI process and then
-/// creates a mesh::Geometry object.
+/// distributes the 'node' coordinate data to the required MPI process
+/// and then creates a mesh::Geometry object.
 ///
 /// @param[in] comm The MPI communicator to build the Geometry on
 /// @param[in] topology The mesh topology
@@ -142,7 +158,7 @@ class Geometry
 /// @param[in] dim The geometric dimension (1, 2, or 3)
 /// @param[in] reorder_fn Function for re-ordering the degree-of-freedom
 /// map associated with the geometry data
-mesh::Geometry
+mesh::Geometry<double>
 create_geometry(MPI_Comm comm, const Topology& topology,
                 const fem::CoordinateElement& element,
                 const graph::AdjacencyList<std::int64_t>& cells,
@@ -159,7 +175,7 @@ create_geometry(MPI_Comm comm, const Topology& topology,
 /// entity in the parent topology
 /// @return A sub-geometry and a map from sub-geometry coordinate
 /// degree-of-freedom to the coordinate degree-of-freedom in `geometry`.
-std::pair<mesh::Geometry, std::vector<int32_t>>
-create_subgeometry(const Topology& topology, const Geometry& geometry, int dim,
-                   std::span<const std::int32_t> subentity_to_entity);
+std::pair<mesh::Geometry<double>, std::vector<int32_t>>
+create_subgeometry(const Topology& topology, const Geometry<double>& geometry,
+                   int dim, std::span<const std::int32_t> subentity_to_entity);
 } // namespace dolfinx::mesh
diff --git a/cpp/dolfinx/mesh/Mesh.cpp b/cpp/dolfinx/mesh/Mesh.cpp
index 656eda44ed2..c94e2226e1e 100644
--- a/cpp/dolfinx/mesh/Mesh.cpp
+++ b/cpp/dolfinx/mesh/Mesh.cpp
@@ -216,9 +216,9 @@ const Topology& Mesh::topology() const { return _topology; }
 //-----------------------------------------------------------------------------
 Topology& Mesh::topology_mutable() const { return _topology; }
 //-----------------------------------------------------------------------------
-Geometry& Mesh::geometry() { return _geometry; }
+Geometry<double>& Mesh::geometry() { return _geometry; }
 //-----------------------------------------------------------------------------
-const Geometry& Mesh::geometry() const { return _geometry; }
+const Geometry<double>& Mesh::geometry() const { return _geometry; }
 //-----------------------------------------------------------------------------
 MPI_Comm Mesh::comm() const { return _comm.comm(); }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/mesh/Mesh.h b/cpp/dolfinx/mesh/Mesh.h
index bf81f3a1fc9..0c139733b1a 100644
--- a/cpp/dolfinx/mesh/Mesh.h
+++ b/cpp/dolfinx/mesh/Mesh.h
@@ -37,7 +37,8 @@ class Mesh
   /// @param[in] comm MPI Communicator
   /// @param[in] topology Mesh topology
   /// @param[in] geometry Mesh geometry
-  template <std::convertible_to<Topology> U, std::convertible_to<Geometry> V>
+  template <std::convertible_to<Topology> U,
+            std::convertible_to<Geometry<double>> V>
   Mesh(MPI_Comm comm, U&& topology, V&& geometry)
       : _topology(std::forward<U>(topology)),
         _geometry(std::forward<V>(geometry)), _comm(comm)
@@ -80,17 +81,17 @@ class Mesh
 
   /// Get mesh geometry
   /// @return The geometry object associated with the mesh
-  Geometry& geometry();
+  Geometry<double>& geometry();
 
   /// Get mesh geometry (const version)
   /// @return The geometry object associated with the mesh
-  const Geometry& geometry() const;
+  const Geometry<double>& geometry() const;
 
   /// Mesh MPI communicator
   /// @return The communicator on which the mesh is distributed
   MPI_Comm comm() const;
 
-  /// Name
+  /// Nameg
   std::string name = "mesh";
 
 private:
@@ -102,7 +103,7 @@ class Mesh
   mutable Topology _topology;
 
   // Mesh geometry
-  Geometry _geometry;
+  Geometry<double> _geometry;
 
   // MPI communicator
   dolfinx::MPI::Comm _comm;
diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp
index 770643e5700..58eaea83050 100644
--- a/cpp/dolfinx/mesh/utils.cpp
+++ b/cpp/dolfinx/mesh/utils.cpp
@@ -500,7 +500,7 @@ mesh::entities_to_geometry(const Mesh& mesh, int dim,
   if (orient and (cell_type != CellType::tetrahedron or dim != 2))
     throw std::runtime_error("Can only orient facets of a tetrahedral mesh");
 
-  const Geometry& geometry = mesh.geometry();
+  const Geometry<double>& geometry = mesh.geometry();
   auto x = geometry.x();
 
   const Topology& topology = mesh.topology();
diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp
index bca527c919d..a2f880e2253 100644
--- a/python/dolfinx/wrappers/mesh.cpp
+++ b/python/dolfinx/wrappers/mesh.cpp
@@ -250,25 +250,27 @@ void mesh(py::module& m)
       .value("shared_vertex", dolfinx::mesh::GhostMode::shared_vertex);
 
   // dolfinx::mesh::Geometry class
-  py::class_<dolfinx::mesh::Geometry, std::shared_ptr<dolfinx::mesh::Geometry>>(
+  py::class_<dolfinx::mesh::Geometry<double>,
+             std::shared_ptr<dolfinx::mesh::Geometry<double>>>(
       m, "Geometry", "Geometry object")
-      .def_property_readonly("dim", &dolfinx::mesh::Geometry::dim,
+      .def_property_readonly("dim", &dolfinx::mesh::Geometry<double>::dim,
                              "Geometric dimension")
-      .def_property_readonly("dofmap", &dolfinx::mesh::Geometry::dofmap)
-      .def("index_map", &dolfinx::mesh::Geometry::index_map)
+      .def_property_readonly("dofmap", &dolfinx::mesh::Geometry<double>::dofmap)
+      .def("index_map", &dolfinx::mesh::Geometry<double>::index_map)
       .def_property_readonly(
           "x",
-          [](const dolfinx::mesh::Geometry& self)
+          [](const dolfinx::mesh::Geometry<double>& self)
           {
             std::array<std::size_t, 2> shape = {self.x().size() / 3, 3};
             return py::array_t<double>(shape, self.x().data(), py::cast(self));
           },
           "Return coordinates of all geometry points. Each row is the "
           "coordinate of a point.")
-      .def_property_readonly("cmap", &dolfinx::mesh::Geometry::cmap,
+      .def_property_readonly("cmap", &dolfinx::mesh::Geometry<double>::cmap,
                              "The coordinate map")
-      .def_property_readonly("input_global_indices",
-                             &dolfinx::mesh::Geometry::input_global_indices);
+      .def_property_readonly(
+          "input_global_indices",
+          &dolfinx::mesh::Geometry<double>::input_global_indices);
 
   // dolfinx::mesh::TopologyComputation
   m.def(
@@ -340,7 +342,7 @@ void mesh(py::module& m)
       .def(py::init(
                [](const MPICommWrapper comm,
                   const dolfinx::mesh::Topology& topology,
-                  dolfinx::mesh::Geometry& geometry)
+                  dolfinx::mesh::Geometry<double>& geometry)
                { return dolfinx::mesh::Mesh(comm.get(), topology, geometry); }),
            py::arg("comm"), py::arg("topology"), py::arg("geometry"))
       .def_property_readonly(