diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
index 481efb4b46..ae761140e9 100644
--- a/cpp/demo/biharmonic/main.cpp
+++ b/cpp/demo/biharmonic/main.cpp
@@ -134,6 +134,7 @@
 #include <dolfinx/common/types.h>
 #include <dolfinx/fem/Constant.h>
 #include <dolfinx/fem/petsc.h>
+#include <memory>
 #include <numbers>
 #include <utility>
 #include <vector>
@@ -169,8 +170,9 @@ int main(int argc, char* argv[])
         basix::element::dpc_variant::unset, false);
 
     //  Create function space
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element, 1)));
 
     // The source function $f$ and the penalty term $\alpha$ are
     // declared:
diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp
index eae77f3671..3789952fae 100644
--- a/cpp/demo/codim_0_assembly/main.cpp
+++ b/cpp/demo/codim_0_assembly/main.cpp
@@ -41,8 +41,9 @@ int main(int argc, char* argv[])
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element, 1)));
 
     // Next we find all cells of the mesh with y<0.5
     const int tdim = mesh->topology()->dim();
@@ -85,8 +86,9 @@ int main(int argc, char* argv[])
     }
 
     // We create the function space used for the trial space
-    auto W = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(submesh, element, {}));
+    auto W
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            submesh, std::make_shared<fem::FiniteElement<U>>(element, 1)));
 
     // A mixed-domain form has functions defined over different meshes.
     // The mesh associated with the measure (dx, ds, etc.) is called the
diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp
index 95f5bb9bed..f589c81902 100644
--- a/cpp/demo/custom_kernel/main.cpp
+++ b/cpp/demo/custom_kernel/main.cpp
@@ -209,8 +209,8 @@ void assemble(MPI_Comm comm)
   mdspand_t<const T, 2> X(X_b.data(), weights.size(), 2);
 
   // Create a scalar function space
-  auto V = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+      mesh, std::make_shared<fem::FiniteElement<T>>(e, 1)));
 
   // Build list of cells to assembler over (all cells owned by this
   // rank)
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 470e9f9ba0..d4d728205e 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -166,8 +166,9 @@ int main(int argc, char* argv[])
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element, std::vector<std::size_t>{3}));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element, 3)));
 
     auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
     auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
@@ -272,8 +273,10 @@ int main(int argc, char* argv[])
     basix::FiniteElement S_element = basix::create_element<U>(
         family, cell_type, k, basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, discontinuous);
-    auto S = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
-        mesh, S_element, std::vector<std::size_t>{3, 3}));
+    auto S
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(S_element, 1),
+            std::vector<std::size_t>{3, 3}));
 
     fem::Function<T> sigma(S);
     sigma.name = "cauchy_stress";
diff --git a/cpp/demo/interpolation-io/main.cpp b/cpp/demo/interpolation-io/main.cpp
index 4392a50627..87611af941 100644
--- a/cpp/demo/interpolation-io/main.cpp
+++ b/cpp/demo/interpolation-io/main.cpp
@@ -46,8 +46,8 @@ void interpolate_scalar(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a scalar function space
-  auto V = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      mesh, std::make_shared<fem::FiniteElement<U>>(e, 1)));
 
   // Create a finite element Function
   auto u = std::make_shared<fem::Function<T>>(V);
@@ -98,8 +98,8 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a Nedelec function space
-  auto V = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      mesh, std::make_shared<fem::FiniteElement<U>>(e, 1)));
 
   // Create a Nedelec finite element Function
   auto u = std::make_shared<fem::Function<T>>(V);
@@ -167,8 +167,9 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, true);
 
   // Create a function space
-  auto V_l = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e_l, std::vector<std::size_t>{2}));
+  auto V_l
+      = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+          mesh, std::make_shared<fem::FiniteElement<U>>(e_l, 2)));
 
   auto u_l = std::make_shared<fem::Function<T>>(V_l);
 
diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index c83cd0852f..75ff553f48 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -15,80 +15,82 @@
 using namespace dolfinx;
 using T = double;
 
-int main(int argc, char* argv[])
+// int main(int argc, char* argv[])
+int main()
 {
-  init_logging(argc, argv);
-  MPI_Init(&argc, &argv);
-  {
-    MPI_Comm comm = MPI_COMM_WORLD;
+  /*
+init_logging(argc, argv);
+MPI_Init(&argc, &argv);
+{
+  MPI_Comm comm = MPI_COMM_WORLD;
 
-    // Create a tetrahedral mesh
-    auto mesh_tet = std::make_shared<mesh::Mesh<double>>(
-        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {20, 20, 20},
-                         mesh::CellType::tetrahedron));
+  // Create a tetrahedral mesh
+  auto mesh_tet = std::make_shared<mesh::Mesh<double>>(
+      mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {20, 20, 20},
+                       mesh::CellType::tetrahedron));
 
-    // Create a hexahedral mesh
-    auto mesh_hex = std::make_shared<mesh::Mesh<double>>(
-        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {15, 15, 15},
-                         mesh::CellType::hexahedron));
+  // Create a hexahedral mesh
+  auto mesh_hex = std::make_shared<mesh::Mesh<double>>(
+      mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {15, 15, 15},
+                       mesh::CellType::hexahedron));
 
-    basix::FiniteElement element_tet = basix::element::create_lagrange<double>(
-        mesh::cell_type_to_basix_type(mesh_tet->topology()->cell_type()), 1,
-        basix::element::lagrange_variant::equispaced, false);
-    auto V_tet = std::make_shared<fem::FunctionSpace<double>>(
-        fem::create_functionspace(mesh_tet, element_tet,
-                                  std::vector<std::size_t>{3}));
+  basix::FiniteElement element_tet = basix::element::create_lagrange<double>(
+      mesh::cell_type_to_basix_type(mesh_tet->topology()->cell_type()), 1,
+      basix::element::lagrange_variant::equispaced, false);
+  auto V_tet = std::make_shared<fem::FunctionSpace<double>>(
+      fem::create_functionspace(mesh_tet, element_tet,
+                                std::vector<std::size_t>{3}));
 
-    basix::FiniteElement element_hex = basix::element::create_lagrange<double>(
-        mesh::cell_type_to_basix_type(mesh_hex->topology()->cell_type()), 2,
-        basix::element::lagrange_variant::equispaced, false);
-    auto V_hex = std::make_shared<fem::FunctionSpace<double>>(
-        fem::create_functionspace(mesh_hex, element_hex,
-                                  std::vector<std::size_t>{3}));
+  basix::FiniteElement element_hex = basix::element::create_lagrange<double>(
+      mesh::cell_type_to_basix_type(mesh_hex->topology()->cell_type()), 2,
+      basix::element::lagrange_variant::equispaced, false);
+  auto V_hex = std::make_shared<fem::FunctionSpace<double>>(
+      fem::create_functionspace(mesh_hex, element_hex,
+                                std::vector<std::size_t>{3}));
 
-    auto u_tet = std::make_shared<fem::Function<T>>(V_tet);
-    auto u_hex = std::make_shared<fem::Function<T>>(V_hex);
+  auto u_tet = std::make_shared<fem::Function<T>>(V_tet);
+  auto u_hex = std::make_shared<fem::Function<T>>(V_hex);
 
-    auto fun = [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
+  auto fun = [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
+  {
+    std::vector<T> fdata(3 * x.extent(1), 0.0);
+    using dextent = MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>;
+    MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<double, dextent> f(fdata.data(), 3,
+                                                              x.extent(1));
+    for (std::size_t i = 0; i < x.extent(1); ++i)
     {
-      std::vector<T> fdata(3 * x.extent(1), 0.0);
-      using dextent = MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>;
-      MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<double, dextent> f(fdata.data(), 3,
-                                                                x.extent(1));
-      for (std::size_t i = 0; i < x.extent(1); ++i)
-      {
-        f(0, i) = std::cos(10 * x(0, i)) * std::sin(10 * x(2, i));
-        f(1, i) = std::sin(10 * x(0, i)) * std::sin(10 * x(2, i));
-        f(2, i) = std::cos(10 * x(0, i)) * std::cos(10 * x(2, i));
-      }
-      return {std::move(fdata), {3, x.extent(1)}};
-    };
+      f(0, i) = std::cos(10 * x(0, i)) * std::sin(10 * x(2, i));
+      f(1, i) = std::sin(10 * x(0, i)) * std::sin(10 * x(2, i));
+      f(2, i) = std::cos(10 * x(0, i)) * std::cos(10 * x(2, i));
+    }
+    return {std::move(fdata), {3, x.extent(1)}};
+  };
 
-    // Interpolate an expression into u_tet
-    u_tet->interpolate(fun);
+  // Interpolate an expression into u_tet
+  u_tet->interpolate(fun);
 
-    // Interpolate from u_tet to u_hex
-    auto cell_map
-        = mesh_hex->topology()->index_map(mesh_hex->topology()->dim());
-    assert(cell_map);
-    std::vector<std::int32_t> cells(
-        cell_map->size_local() + cell_map->num_ghosts(), 0);
-    std::iota(cells.begin(), cells.end(), 0);
-    geometry::PointOwnershipData<T> interpolation_data
-        = fem::create_interpolation_data(
-            u_hex->function_space()->mesh()->geometry(),
-            *u_hex->function_space()->element(),
-            *u_tet->function_space()->mesh(), std::span(cells), 1e-8);
-    u_hex->interpolate(*u_tet, cells, interpolation_data);
+  // Interpolate from u_tet to u_hex
+  auto cell_map
+      = mesh_hex->topology()->index_map(mesh_hex->topology()->dim());
+  assert(cell_map);
+  std::vector<std::int32_t> cells(
+      cell_map->size_local() + cell_map->num_ghosts(), 0);
+  std::iota(cells.begin(), cells.end(), 0);
+  geometry::PointOwnershipData<T> interpolation_data
+      = fem::create_interpolation_data(
+          u_hex->function_space()->mesh()->geometry(),
+          *u_hex->function_space()->element(),
+          *u_tet->function_space()->mesh(), std::span(cells), 1e-8);
+  u_hex->interpolate(*u_tet, cells, interpolation_data);
 
 #ifdef HAS_ADIOS2
-    io::VTXWriter<double> write_tet(mesh_tet->comm(), "u_tet.bp", {u_tet});
-    write_tet.write(0.0);
-    io::VTXWriter<double> write_hex(mesh_hex->comm(), "u_hex.bp", {u_hex});
-    write_hex.write(0.0);
+  io::VTXWriter<double> write_tet(mesh_tet->comm(), "u_tet.bp", {u_tet});
+  write_tet.write(0.0);
+  io::VTXWriter<double> write_hex(mesh_hex->comm(), "u_hex.bp", {u_hex});
+  write_hex.write(0.0);
 #endif
-  }
-  MPI_Finalize();
-
+}
+MPI_Finalize();
+*/
   return 0;
 }
diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp
index cb6373ce28..7897d2c552 100644
--- a/cpp/demo/mixed_poisson/main.cpp
+++ b/cpp/demo/mixed_poisson/main.cpp
@@ -67,13 +67,12 @@ int main(int argc, char* argv[])
     // auto V = std::make_shared<fem::FunctionSpace<U>>(mesh, ME, dofmap, vs);
     // auto V = std::make_shared<fem::FunctionSpace<U>>(
     //     fem::create_functionspace<U>(mesh, ME));
-    std::vector<
-        std::tuple<std::reference_wrapper<const basix::FiniteElement<U>>,
-                   std::size_t, bool>>
-        ME1{{RT, 1, false}, {P0, 1, false}};
+    // std::vector<
+    //     std::tuple<std::reference_wrapper<const basix::FiniteElement<U>>,
+    //                std::size_t, bool>>
+    //     ME1{{RT, 1, false}, {P0, 1, false}};
     auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace<U>(mesh, ME1));
-    /*
+        fem::create_functionspace<U>(mesh, ME));
 
     // Get subspaces (views into V)
     auto V0 = std::make_shared<fem::FunctionSpace<U>>(V->sub({0}));
@@ -232,7 +231,6 @@ int main(int argc, char* argv[])
     io::VTXWriter<U> vtx(MPI_COMM_WORLD, "u.bp", {u_soln}, "bp4");
     vtx.write(0);
 #endif
-*/
   }
 
   PetscFinalize();
diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
index f392c49b67..c85a0e4ff2 100644
--- a/cpp/demo/poisson/main.cpp
+++ b/cpp/demo/poisson/main.cpp
@@ -124,8 +124,9 @@ int main(int argc, char* argv[])
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element, 1)));
 
     //  Next, we define the variational formulation by initializing the
     //  bilinear and linear forms ($a$, $L$) using the previously
diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp
index 426d62f1f5..7187dba821 100644
--- a/cpp/demo/poisson_matrix_free/main.cpp
+++ b/cpp/demo/poisson_matrix_free/main.cpp
@@ -143,8 +143,8 @@ void solver(MPI_Comm comm)
       basix::element::family::P, basix::cell::type::triangle, 2,
       basix::element::lagrange_variant::unset,
       basix::element::dpc_variant::unset, false);
-  auto V = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, element, {}));
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      mesh, std::make_shared<fem::FiniteElement<U>>(element, 1)));
 
   // Prepare and set Constants for the bilinear form
   auto f = std::make_shared<fem::Constant<T>>(-6.0);
diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h
index 8613c7e983..a1a0a042c5 100644
--- a/cpp/dolfinx/fem/assemble_matrix_impl.h
+++ b/cpp/dolfinx/fem/assemble_matrix_impl.h
@@ -24,6 +24,7 @@
 
 namespace dolfinx::fem::impl
 {
+/// @brief Typedef
 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
     const std::int32_t,
     MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 1e73834258..3c91eac8b6 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -784,42 +784,42 @@ Form<T, U> create_form(
   return L;
 }
 
-/// @brief Create a function space from a Basix finite element.
-/// @param[in] mesh Mesh.
-/// @param[in] e Basix finite element.
-/// @param[in] value_shape Value shape for 'blocked' elements, e.g.
-/// vector-valued Lagrange elements where each component for the vector
-/// field is a Lagrange element. For example, a vector-valued element in
-/// 3D will have `value_shape` equal to `{3}`, and for a second-order
-/// tensor element in 2D `value_shape` equal to `{2, 2}`.
-/// @param[in] reorder_fn The graph reordering function to call on the
-/// dofmap. If `nullptr`, the default re-ordering is used.
-/// @return The created function space.
-template <std::floating_point T>
-FunctionSpace<T> create_functionspace(
-    std::shared_ptr<mesh::Mesh<T>> mesh, const basix::FiniteElement<T>& e,
-    std::optional<std::vector<std::size_t>> value_shape = std::nullopt,
-    std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
-        reorder_fn
-    = nullptr)
-{
-  if (mesh::cell_type_from_basix_type(e.cell_type())
-      != mesh->topology()->cell_type())
-  {
-    throw std::runtime_error("Cell type of element and mesh must match.");
-  }
-
-  // Create a DOLFINx element
-  std::size_t bs
-      = value_shape.has_value()
-            ? std::accumulate(value_shape->begin(), value_shape->end(), 1,
-                              std::multiplies{})
-            : 1;
-  auto _e = std::make_shared<const FiniteElement<T>>(e, bs);
-  assert(_e);
-
-  return create_functionspace(mesh, _e, std::nullopt, reorder_fn);
-}
+// /// @brief Create a function space from a Basix finite element.
+// /// @param[in] mesh Mesh.
+// /// @param[in] e Basix finite element.
+// /// @param[in] value_shape Value shape for 'blocked' elements, e.g.
+// /// vector-valued Lagrange elements where each component for the vector
+// /// field is a Lagrange element. For example, a vector-valued element in
+// /// 3D will have `value_shape` equal to `{3}`, and for a second-order
+// /// tensor element in 2D `value_shape` equal to `{2, 2}`.
+// /// @param[in] reorder_fn The graph reordering function to call on the
+// /// dofmap. If `nullptr`, the default re-ordering is used.
+// /// @return The created function space.
+// template <std::floating_point T>
+// FunctionSpace<T> create_functionspace(
+//     std::shared_ptr<mesh::Mesh<T>> mesh, const basix::FiniteElement<T>& e,
+//     std::optional<std::vector<std::size_t>> value_shape = std::nullopt,
+//     std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
+//         reorder_fn
+//     = nullptr)
+// {
+//   if (mesh::cell_type_from_basix_type(e.cell_type())
+//       != mesh->topology()->cell_type())
+//   {
+//     throw std::runtime_error("Cell type of element and mesh must match.");
+//   }
+
+//   // Create a DOLFINx element
+//   std::size_t bs
+//       = value_shape.has_value()
+//             ? std::accumulate(value_shape->begin(), value_shape->end(), 1,
+//                               std::multiplies{})
+//             : 1;
+//   auto _e = std::make_shared<const FiniteElement<T>>(e, bs);
+//   assert(_e);
+
+//   return create_functionspace(mesh, _e, std::nullopt, reorder_fn);
+// }
 
 /// @brief NEW Create a function space from a fem::FiniteElement.
 template <std::floating_point T>
@@ -869,20 +869,20 @@ FunctionSpace<T> create_functionspace(
 /// @param[in] reorder_fn The graph reordering function to call on the
 /// dofmap. If `nullptr`, the default re-ordering is used.
 /// @return The created function space.
-template <std::floating_point T>
-FunctionSpace<T> create_functionspace(
-    std::shared_ptr<mesh::Mesh<T>> mesh,
-    std::vector<
-        std::tuple<std::reference_wrapper<const basix::FiniteElement<T>>,
-                   std::size_t, bool>>
-        e,
-    std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
-        reorder_fn
-    = nullptr)
-{
-  auto _e = std::make_shared<fem::FiniteElement<T>>(e);
-  return create_functionspace<T>(mesh, _e, std::nullopt, reorder_fn);
-}
+// template <std::floating_point T>
+// FunctionSpace<T> create_functionspace(
+//     std::shared_ptr<mesh::Mesh<T>> mesh,
+//     std::vector<
+//         std::tuple<std::reference_wrapper<const basix::FiniteElement<T>>,
+//                    std::size_t, bool>>
+//         e,
+//     std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
+//         reorder_fn
+//     = nullptr)
+// {
+//   auto _e = std::make_shared<fem::FiniteElement<T>>(e);
+//   return create_functionspace<T>(mesh, _e, std::nullopt, reorder_fn);
+// }
 
 /// @private
 namespace impl