diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
index 04b2477c6a4..8db8c934b95 100644
--- a/cpp/demo/biharmonic/main.cpp
+++ b/cpp/demo/biharmonic/main.cpp
@@ -128,6 +128,7 @@
 // convenience we also include the DOLFINx namespace.
 
 #include "biharmonic.h"
+#include <basix/finite-element.h>
 #include <cmath>
 #include <dolfinx.h>
 #include <dolfinx/common/types.h>
@@ -162,9 +163,14 @@ int main(int argc, char* argv[])
     //    A function space object, which is defined in the generated code,
     //    is created:
 
+    auto element = basix::create_element<U>(
+        basix::element::family::P, basix::cell::type::triangle, 2,
+        basix::element::lagrange_variant::unset,
+        basix::element::dpc_variant::unset, false);
+
     //  Create function space
     auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(functionspace_form_biharmonic_a, "u", mesh));
+        fem::create_functionspace(mesh, element));
 
     // The source function $f$ and the penalty term $\alpha$ are
     // declared:
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 35da5a18895..666220dd414 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -153,8 +153,13 @@ int main(int argc, char* argv[])
         mesh::CellType::tetrahedron,
         mesh::create_cell_partitioner(mesh::GhostMode::none)));
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
-        functionspace_form_hyperelasticity_F_form, "u", mesh));
+    auto element = basix::create_element<U>(
+        basix::element::family::P, basix::cell::type::tetrahedron, 1,
+        basix::element::lagrange_variant::unset,
+        basix::element::dpc_variant::unset, false);
+
+    auto V = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace(mesh, 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});
diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
index 0a296adb73f..a92d60a2d8d 100644
--- a/cpp/demo/poisson/main.cpp
+++ b/cpp/demo/poisson/main.cpp
@@ -80,6 +80,7 @@
 // namespace.
 
 #include "poisson.h"
+#include <basix/finite-element.h>
 #include <cmath>
 #include <dolfinx.h>
 #include <dolfinx/fem/Constant.h>
@@ -114,8 +115,13 @@ int main(int argc, char* argv[])
         mesh::create_rectangle<U>(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}},
                                   {32, 16}, mesh::CellType::triangle, part));
 
+    auto element = basix::create_element<U>(
+        basix::element::family::P, basix::cell::type::triangle, 1,
+        basix::element::lagrange_variant::unset,
+        basix::element::dpc_variant::unset, false);
+
     auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
+        fem::create_functionspace(mesh, element, {}));
 
     //  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 3281b664ce4..2488c9fcb2c 100644
--- a/cpp/demo/poisson_matrix_free/main.cpp
+++ b/cpp/demo/poisson_matrix_free/main.cpp
@@ -42,6 +42,7 @@
 // ## C++ program
 
 #include "poisson.h"
+#include <basix/finite-element.h>
 #include <cmath>
 #include <complex>
 #include <concepts>
@@ -137,8 +138,12 @@ void solver(MPI_Comm comm)
   auto mesh = std::make_shared<mesh::Mesh<U>>(mesh::create_rectangle<U>(
       comm, {{{0.0, 0.0}, {1.0, 1.0}}}, {10, 10}, mesh::CellType::triangle,
       mesh::create_cell_partitioner(mesh::GhostMode::none)));
+  auto element = basix::create_element<U>(
+      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(functionspace_form_poisson_M, "ui", mesh));
+      fem::create_functionspace(mesh, element, {}));
 
   // Prepare and set Constants for the bilinear form
   auto f = std::make_shared<fem::Constant<T>>(-6.0);
diff --git a/cpp/dolfinx/fem/FiniteElement.cpp b/cpp/dolfinx/fem/FiniteElement.cpp
index af34bf25799..00d2f05996e 100644
--- a/cpp/dolfinx/fem/FiniteElement.cpp
+++ b/cpp/dolfinx/fem/FiniteElement.cpp
@@ -13,7 +13,6 @@
 #include <dolfinx/common/log.h>
 #include <functional>
 #include <numeric>
-#include <ufcx.h>
 #include <utility>
 #include <vector>
 
@@ -22,64 +21,11 @@ using namespace dolfinx::fem;
 
 namespace
 {
-//-----------------------------------------------------------------------------
-
-/// Check if an element is a Basix element (or a blocked element
-/// containing a Basix element)
-bool is_basix_element(const ufcx_finite_element& element)
-{
-  if (element.element_type == ufcx_basix_element)
-    return true;
-  else if (element.block_size != 1)
-  {
-    // TODO: what should happen if the element is a blocked element
-    // containing a blocked element containing a Basix element?
-    return element.sub_elements[0]->element_type == ufcx_basix_element;
-  }
-  else
-    return false;
-}
-//-----------------------------------------------------------------------------
-
-/// Check if an element is a quadrature element (or a blocked element
-/// containing a quadrature element)
-bool is_quadrature_element(const ufcx_finite_element& element)
-{
-  if (element.element_type == ufcx_quadrature_element)
-    return true;
-  else if (element.block_size != 1)
-  {
-    // TODO: what should happen if the element is a blocked element
-    // containing a blocked element containing a quadrature element?
-    return element.sub_elements[0]->element_type == ufcx_quadrature_element;
-  }
-  else
-    return false;
-}
-//-----------------------------------------------------------------------------
-
-/// Check if an element is a custom Basix element (or a blocked element
-/// containing a custom Basix element)
-bool is_basix_custom_element(const ufcx_finite_element& element)
-{
-  if (element.element_type == ufcx_basix_custom_element)
-    return true;
-  else if (element.block_size != 1)
-  {
-    // TODO: what should happen if the element is a blocked element
-    // containing a blocked element containing a Basix element?
-    return element.sub_elements[0]->element_type == ufcx_basix_custom_element;
-  }
-  else
-    return false;
-}
-//-----------------------------------------------------------------------------
-
 /// Recursively extract sub finite element
 template <std::floating_point T>
 std::shared_ptr<const FiniteElement<T>>
 _extract_sub_element(const FiniteElement<T>& finite_element,
-                     const std::vector<int>& component)
+                     std::span<const int> component)
 {
   // Check that a sub system has been specified
   if (component.empty())
@@ -115,180 +61,35 @@ _extract_sub_element(const FiniteElement<T>& finite_element,
 
   return _extract_sub_element(*sub_element, sub_component);
 }
-//-----------------------------------------------------------------------------
-
 } // namespace
 
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-FiniteElement<T>::FiniteElement(const ufcx_finite_element& e)
-    : _signature(e.signature), _space_dim(e.space_dimension),
-      _reference_value_shape(e.reference_value_shape,
-                             e.reference_value_shape + e.reference_value_rank),
-      _bs(e.block_size), _is_mixed(e.element_type == ufcx_mixed_element),
-      _symmetric(e.symmetric)
-{
-  const ufcx_shape _shape = e.cell_shape;
-  switch (_shape)
-  {
-  case interval:
-    _cell_shape = mesh::CellType::interval;
-    break;
-  case triangle:
-    _cell_shape = mesh::CellType::triangle;
-    break;
-  case quadrilateral:
-    _cell_shape = mesh::CellType::quadrilateral;
-    break;
-  case tetrahedron:
-    _cell_shape = mesh::CellType::tetrahedron;
-    break;
-  case prism:
-    _cell_shape = mesh::CellType::prism;
-    break;
-  case hexahedron:
-    _cell_shape = mesh::CellType::hexahedron;
-    break;
-  default:
-    throw std::runtime_error(
-        "Unknown UFC cell type when building FiniteElement.");
-  }
-  assert(mesh::cell_dim(_cell_shape) == e.topological_dimension);
-
-  static const std::map<ufcx_shape, std::string> ufcx_to_cell
-      = {{vertex, "point"},         {interval, "interval"},
-         {triangle, "triangle"},    {tetrahedron, "tetrahedron"},
-         {prism, "prism"},          {quadrilateral, "quadrilateral"},
-         {hexahedron, "hexahedron"}};
-  const std::string cell_shape = ufcx_to_cell.at(e.cell_shape);
-
-  _needs_dof_transformations = false;
-  _needs_dof_permutations = false;
-  // Create all sub-elements
-  for (int i = 0; i < e.num_sub_elements; ++i)
-  {
-    ufcx_finite_element* ufcx_sub_element = e.sub_elements[i];
-    _sub_elements.push_back(
-        std::make_shared<FiniteElement<T>>(*ufcx_sub_element));
-    if (_sub_elements[i]->needs_dof_permutations())
-      _needs_dof_permutations = true;
-    if (_sub_elements[i]->needs_dof_transformations())
-      _needs_dof_transformations = true;
-  }
-
-  if (is_basix_custom_element(e))
-  {
-    // Recreate the custom Basix element using information written into
-    // the generated code
-    ufcx_basix_custom_finite_element* ce = e.custom_element;
-    const basix::cell::type cell_type
-        = static_cast<basix::cell::type>(ce->cell_type);
-
-    const std::vector<std::size_t> value_shape(
-        ce->value_shape, ce->value_shape + ce->value_shape_length);
-    const std::size_t value_size = std::reduce(
-        value_shape.begin(), value_shape.end(), 1, std::multiplies{});
-
-    const int nderivs = ce->interpolation_nderivs;
-    const std::size_t nderivs_dim = basix::polyset::nderivs(cell_type, nderivs);
-
-    using array2_t = std::pair<std::vector<T>, std::array<std::size_t, 2>>;
-    using array4_t = std::pair<std::vector<T>, std::array<std::size_t, 4>>;
-    std::array<std::vector<array2_t>, 4> x;
-    std::array<std::vector<array4_t>, 4> M;
-    { // scope
-      int pt_n = 0;
-      int p_e = 0;
-      int m_e = 0;
-      const std::size_t dim = static_cast<std::size_t>(
-          basix::cell::topological_dimension(cell_type));
-      for (std::size_t d = 0; d <= dim; ++d)
-      {
-        const int num_entities = basix::cell::num_sub_entities(cell_type, d);
-        for (int entity = 0; entity < num_entities; ++entity)
-        {
-          std::size_t npts = ce->npts[pt_n + entity];
-          std::size_t ndofs = ce->ndofs[pt_n + entity];
-
-          std::array pshape = {npts, dim};
-          auto& pts
-              = x[d].emplace_back(std::vector<T>(pshape[0] * pshape[1]), pshape)
-                    .first;
-          std::copy_n(ce->x + p_e, pts.size(), pts.begin());
-          p_e += pts.size();
-
-          std::array mshape = {ndofs, value_size, npts, nderivs_dim};
-          std::size_t msize
-              = std::reduce(mshape.begin(), mshape.end(), 1, std::multiplies{});
-          auto& mat = M[d].emplace_back(std::vector<T>(msize), mshape).first;
-          std::copy_n(ce->M + m_e, mat.size(), mat.begin());
-          m_e += mat.size();
-        }
-
-        pt_n += num_entities;
-      }
-    }
-
-    using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
-        const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
-    using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
-        const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
-
-    std::array<std::vector<cmdspan2_t>, 4> _x;
-    for (std::size_t i = 0; i < x.size(); ++i)
-      for (auto& [buffer, shape] : x[i])
-        _x[i].push_back(cmdspan2_t(buffer.data(), shape));
-
-    std::array<std::vector<cmdspan4_t>, 4> _M;
-    for (std::size_t i = 0; i < M.size(); ++i)
-      for (auto& [buffer, shape] : M[i])
-        _M[i].push_back(cmdspan4_t(buffer.data(), shape));
-
-    std::vector<T> wcoeffs_b(ce->wcoeffs_rows * ce->wcoeffs_cols);
-    cmdspan2_t wcoeffs(wcoeffs_b.data(), ce->wcoeffs_rows, ce->wcoeffs_cols);
-    std::copy_n(ce->wcoeffs, wcoeffs_b.size(), wcoeffs_b.begin());
-
-    _element = std::make_unique<basix::FiniteElement<T>>(
-        basix::create_custom_element(
-            cell_type, value_shape, wcoeffs, _x, _M, nderivs,
-            static_cast<basix::maps::type>(ce->map_type),
-            static_cast<basix::sobolev::space>(ce->sobolev_space),
-            ce->discontinuous, ce->embedded_subdegree, ce->embedded_superdegree,
-            static_cast<basix::polyset::type>(ce->polyset_type)));
-    _needs_dof_transformations
-        = !_element->dof_transformations_are_identity()
-          and !_element->dof_transformations_are_permutations();
-
-    _needs_dof_permutations
-        = !_element->dof_transformations_are_identity()
-          and _element->dof_transformations_are_permutations();
-  }
-  else if (is_basix_element(e))
+FiniteElement<T>::FiniteElement(mesh::CellType cell_type,
+                                std::span<const geometry_type> points,
+                                std::array<std::size_t, 2> pshape,
+                                const std::size_t block_size)
+    : _signature("Quadrature element " + std::to_string(pshape[0]) + " "
+                 + std::to_string(block_size)),
+      _space_dim(pshape[0] * block_size), _reference_value_shape({}),
+      _bs(block_size), _is_mixed(false), _symmetric(false),
+      _needs_dof_permutations(false), _needs_dof_transformations(false),
+      _entity_dofs(mesh::cell_dim(cell_type) + 1),
+      _entity_closure_dofs(mesh::cell_dim(cell_type) + 1),
+      _points(std::vector<T>(points.begin(), points.end()), pshape)
+{
+  const int tdim = mesh::cell_dim(cell_type);
+  for (int d = 0; d <= tdim; ++d)
   {
-    _element
-        = std::make_unique<basix::FiniteElement<T>>(basix::create_element<T>(
-            static_cast<basix::element::family>(e.basix_family),
-            static_cast<basix::cell::type>(e.basix_cell), e.degree,
-            static_cast<basix::element::lagrange_variant>(e.lagrange_variant),
-            static_cast<basix::element::dpc_variant>(e.dpc_variant),
-            e.discontinuous));
-    _needs_dof_transformations
-        = !_element->dof_transformations_are_identity()
-          and !_element->dof_transformations_are_permutations();
-    _needs_dof_permutations
-        = !_element->dof_transformations_are_identity()
-          and _element->dof_transformations_are_permutations();
+    int num_entities = mesh::cell_num_entities(cell_type, d);
+    _entity_dofs[d].resize(num_entities);
+    _entity_closure_dofs[d].resize(num_entities);
   }
 
-  if (is_quadrature_element(e))
+  for (std::size_t i = 0; i < pshape[0]; ++i)
   {
-    assert(e.custom_quadrature);
-    ufcx_quadrature_rule* qr = e.custom_quadrature;
-    std::size_t npts = qr->npts;
-    std::size_t tdim = qr->topological_dimension;
-    std::array qshape = {npts, tdim};
-    _points = std::make_pair(std::vector<T>(qshape[0] * qshape[1]), qshape);
-    std::copy_n(qr->points, _points.first.size(), _points.first.begin());
+    _entity_dofs[tdim][0].push_back(i);
+    _entity_closure_dofs[tdim][0].push_back(i);
   }
 }
 //-----------------------------------------------------------------------------
@@ -296,29 +97,29 @@ template <std::floating_point T>
 FiniteElement<T>::FiniteElement(const basix::FiniteElement<T>& element,
                                 const std::size_t block_size,
                                 const bool symmetric)
-    : _reference_value_shape(element.value_shape()), _bs(block_size),
-      _is_mixed(false), _symmetric(symmetric)
+    : _space_dim(block_size * element.dim()),
+      _reference_value_shape(element.value_shape()), _bs(block_size),
+      _is_mixed(false),
+      _element(std::make_unique<basix::FiniteElement<T>>(element)),
+      _symmetric(symmetric),
+      _needs_dof_permutations(
+          !_element->dof_transformations_are_identity()
+          and _element->dof_transformations_are_permutations()),
+      _needs_dof_transformations(
+          !_element->dof_transformations_are_identity()
+          and !_element->dof_transformations_are_permutations()),
+      _entity_dofs(element.entity_dofs()),
+      _entity_closure_dofs(element.entity_closure_dofs())
 {
-  _space_dim = _bs * element.dim();
-
   // Create all sub-elements
   if (_bs > 1)
   {
-    for (int i = 0; i < _bs; ++i)
-    {
-      _sub_elements.push_back(std::make_shared<FiniteElement<T>>(element, 1));
-    }
+    _sub_elements
+        = std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>(
+            _bs, std::make_shared<FiniteElement<T>>(element, 1));
     _reference_value_shape = {block_size};
   }
 
-  _element = std::make_unique<basix::FiniteElement<T>>(element);
-  assert(_element);
-  _needs_dof_transformations
-      = !_element->dof_transformations_are_identity()
-        and !_element->dof_transformations_are_permutations();
-  _needs_dof_permutations
-      = !_element->dof_transformations_are_identity()
-        and _element->dof_transformations_are_permutations();
   std::string family;
   switch (_element->family())
   {
@@ -337,6 +138,65 @@ FiniteElement<T>::FiniteElement(const basix::FiniteElement<T>& element,
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
+FiniteElement<T>::FiniteElement(
+    const std::vector<std::shared_ptr<const FiniteElement<T>>>& elements)
+    : _space_dim(0), _sub_elements(elements), _bs(1), _is_mixed(true),
+      _symmetric(false), _needs_dof_permutations(false),
+      _needs_dof_transformations(false)
+{
+  std::size_t vsize = 0;
+  _signature = "Mixed element (";
+
+  const std::vector<std::vector<std::vector<int>>>& ed
+      = elements[0]->entity_dofs();
+  _entity_dofs.resize(ed.size());
+  _entity_closure_dofs.resize(ed.size());
+  for (std::size_t i = 0; i < ed.size(); ++i)
+  {
+    _entity_dofs[i].resize(ed[i].size());
+    _entity_closure_dofs[i].resize(ed[i].size());
+  }
+
+  int dof_offset = 0;
+  for (auto& e : elements)
+  {
+    vsize += e->reference_value_size();
+    _signature += e->signature() + ", ";
+
+    if (e->needs_dof_permutations())
+      _needs_dof_permutations = true;
+    if (e->needs_dof_transformations())
+      _needs_dof_transformations = true;
+
+    const std::size_t sub_bs = e->block_size();
+    for (std::size_t i = 0; i < _entity_dofs.size(); ++i)
+    {
+      for (std::size_t j = 0; j < _entity_dofs[i].size(); ++j)
+      {
+        const std::vector<int> sub_ed = e->entity_dofs()[i][j];
+        const std::vector<int> sub_ecd = e->entity_closure_dofs()[i][j];
+        for (auto k : sub_ed)
+        {
+          for (std::size_t b = 0; b < sub_bs; ++b)
+            _entity_dofs[i][j].push_back(dof_offset + k * sub_bs + b);
+        }
+        for (auto k : sub_ecd)
+        {
+          for (std::size_t b = 0; b < sub_bs; ++b)
+            _entity_closure_dofs[i][j].push_back(dof_offset + k * sub_bs + b);
+        }
+      }
+    }
+
+    dof_offset += e->space_dimension();
+  }
+
+  _space_dim = dof_offset;
+  _reference_value_shape = {vsize};
+  _signature += ")";
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point T>
 bool FiniteElement<T>::operator==(const FiniteElement& e) const
 {
   if (!_element or !e._element)
@@ -361,21 +221,30 @@ std::string FiniteElement<T>::signature() const noexcept
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-mesh::CellType FiniteElement<T>::cell_shape() const noexcept
+int FiniteElement<T>::space_dimension() const noexcept
 {
-  return _cell_shape;
+  return _space_dim;
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-int FiniteElement<T>::space_dimension() const noexcept
+std::span<const std::size_t>
+FiniteElement<T>::reference_value_shape() const noexcept
 {
-  return _space_dim;
+  return _reference_value_shape;
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-std::span<const std::size_t> FiniteElement<T>::reference_value_shape() const
+const std::vector<std::vector<std::vector<int>>>&
+FiniteElement<T>::entity_dofs() const noexcept
 {
-  return _reference_value_shape;
+  return _entity_dofs;
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point T>
+const std::vector<std::vector<std::vector<int>>>&
+FiniteElement<T>::entity_closure_dofs() const noexcept
+{
+  return _entity_closure_dofs;
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h
index 86cbec6ad8d..9f704b8c034 100644
--- a/cpp/dolfinx/fem/FiniteElement.h
+++ b/cpp/dolfinx/fem/FiniteElement.h
@@ -18,8 +18,6 @@
 #include <utility>
 #include <vector>
 
-struct ufcx_finite_element;
-
 namespace dolfinx::fem
 {
 /// DOF transformation type
@@ -42,17 +40,28 @@ class FiniteElement
   /// Geometry type of the Mesh that the FunctionSpace is defined on.
   using geometry_type = T;
 
-  /// @brief Create finite element from UFC finite element.
-  /// @param[in] e UFC finite element.
-  explicit FiniteElement(const ufcx_finite_element& e);
-
-  /// @brief Create finite element from a Basix finite element.
+  /// @brief Create a finite element from a Basix finite element.
   /// @param[in] element Basix finite element
   /// @param[in] block_size The block size for the element
   /// @param[in] symmetric Is the element a symmetric tensor?
   FiniteElement(const basix::FiniteElement<geometry_type>& element,
                 const std::size_t block_size, const bool symmetric = false);
 
+  /// @brief Create mixed finite element from a list of finite elements.
+  /// @param[in] elements Basix finite elements
+  FiniteElement(
+      const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
+          elements);
+
+  /// @brief Create a quadrature element
+  /// @param[in] cell_type The cell type
+  /// @param[in] points Quadrature points
+  /// @param[in] pshape Shape of points array
+  /// @param[in] block_size The block size for the element
+  FiniteElement(mesh::CellType cell_type, std::span<const geometry_type> points,
+                std::array<std::size_t, 2> pshape,
+                const std::size_t block_size);
+
   /// Copy constructor
   FiniteElement(const FiniteElement& element) = delete;
 
@@ -88,10 +97,6 @@ class FiniteElement
   /// properties.
   std::string signature() const noexcept;
 
-  /// Cell shape
-  /// @return Element cell shape
-  mesh::CellType cell_shape() const noexcept;
-
   /// Dimension of the finite element function space (the number of
   /// degrees-of-freedom for the element)
   /// @return Dimension of the finite element space
@@ -109,7 +114,15 @@ class FiniteElement
   int reference_value_size() const;
 
   /// The reference value shape
-  std::span<const std::size_t> reference_value_shape() const;
+  std::span<const std::size_t> reference_value_shape() const noexcept;
+
+  /// The local DOFs associated with each subentity of the cell
+  const std::vector<std::vector<std::vector<int>>>&
+  entity_dofs() const noexcept;
+
+  /// The local DOFs associated with the closure of each subentity of the cell
+  const std::vector<std::vector<std::vector<int>>>&
+  entity_closure_dofs() const noexcept;
 
   /// Does the element represent a symmetric 2-tensor?
   bool symmetric() const;
@@ -315,25 +328,25 @@ class FiniteElement
         // Mixed element
         std::vector<std::function<void(
             std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
-            sub_element_functions;
+            sub_element_fns;
         std::vector<int> dims;
         for (std::size_t i = 0; i < _sub_elements.size(); ++i)
         {
-          sub_element_functions.push_back(
+          sub_element_fns.push_back(
               _sub_elements[i]->template dof_transformation_fn<U>(ttype));
           dims.push_back(_sub_elements[i]->space_dimension());
         }
 
-        return [dims, sub_element_functions](
-                   std::span<U> data, std::span<const std::uint32_t> cell_info,
-                   std::int32_t cell, int block_size)
+        return [dims, sub_element_fns](std::span<U> data,
+                                       std::span<const std::uint32_t> cell_info,
+                                       std::int32_t cell, int block_size)
         {
           std::size_t offset = 0;
-          for (std::size_t e = 0; e < sub_element_functions.size(); ++e)
+          for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
           {
             const std::size_t width = dims[e] * block_size;
-            sub_element_functions[e](data.subspan(offset, width), cell_info,
-                                     cell, block_size);
+            sub_element_fns[e](data.subspan(offset, width), cell_info, cell,
+                               block_size);
             offset += width;
           }
         };
@@ -343,13 +356,12 @@ class FiniteElement
         // Blocked element
         std::function<void(std::span<U>, std::span<const std::uint32_t>,
                            std::int32_t, int)>
-            sub_function
-            = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
+            sub_fn = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
         const int ebs = _bs;
-        return [ebs, sub_function](std::span<U> data,
-                                   std::span<const std::uint32_t> cell_info,
-                                   std::int32_t cell, int data_block_size)
-        { sub_function(data, cell_info, cell, ebs * data_block_size); };
+        return [ebs, sub_fn](std::span<U> data,
+                             std::span<const std::uint32_t> cell_info,
+                             std::int32_t cell, int data_block_size)
+        { sub_fn(data, cell_info, cell, ebs * data_block_size); };
       }
     }
 
@@ -420,22 +432,22 @@ class FiniteElement
         // Mixed element
         std::vector<std::function<void(
             std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
-            sub_element_functions;
+            sub_element_fns;
         for (std::size_t i = 0; i < _sub_elements.size(); ++i)
         {
-          sub_element_functions.push_back(
+          sub_element_fns.push_back(
               _sub_elements[i]->template dof_transformation_right_fn<U>(ttype));
         }
 
-        return [this, sub_element_functions](
-                   std::span<U> data, std::span<const std::uint32_t> cell_info,
-                   std::int32_t cell, int block_size)
+        return [this, sub_element_fns](std::span<U> data,
+                                       std::span<const std::uint32_t> cell_info,
+                                       std::int32_t cell, int block_size)
         {
           std::size_t offset = 0;
-          for (std::size_t e = 0; e < sub_element_functions.size(); ++e)
+          for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
           {
-            sub_element_functions[e](data.subspan(offset, data.size() - offset),
-                                     cell_info, cell, block_size);
+            sub_element_fns[e](data.subspan(offset, data.size() - offset),
+                               cell_info, cell, block_size);
             offset += _sub_elements[e]->space_dimension();
           }
         };
@@ -449,18 +461,17 @@ class FiniteElement
         // transformation from the left to data using xxxyyyzzz ordering
         std::function<void(std::span<U>, std::span<const std::uint32_t>,
                            std::int32_t, int)>
-            sub_function
-            = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
-        return [this, sub_function](std::span<U> data,
-                                    std::span<const std::uint32_t> cell_info,
-                                    std::int32_t cell, int data_block_size)
+            sub_fn = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
+        return [this, sub_fn](std::span<U> data,
+                              std::span<const std::uint32_t> cell_info,
+                              std::int32_t cell, int data_block_size)
         {
           const int ebs = block_size();
           const std::size_t dof_count = data.size() / data_block_size;
           for (int block = 0; block < data_block_size; ++block)
           {
-            sub_function(data.subspan(block * dof_count, dof_count), cell_info,
-                         cell, ebs);
+            sub_fn(data.subspan(block * dof_count, dof_count), cell_info, cell,
+                   ebs);
           }
         };
       }
@@ -706,8 +717,6 @@ class FiniteElement
 private:
   std::string _signature;
 
-  mesh::CellType _cell_shape;
-
   int _space_dim;
 
   // List of sub-elements (if any)
@@ -724,6 +733,9 @@ class FiniteElement
   // Indicate whether this is a mixed element
   bool _is_mixed;
 
+  // Basix Element (nullptr for mixed elements)
+  std::unique_ptr<basix::FiniteElement<geometry_type>> _element;
+
   // Indicate whether this element represents a symmetric 2-tensor
   bool _symmetric;
 
@@ -731,8 +743,8 @@ class FiniteElement
   bool _needs_dof_permutations;
   bool _needs_dof_transformations;
 
-  // Basix Element (nullptr for mixed elements)
-  std::unique_ptr<basix::FiniteElement<geometry_type>> _element;
+  std::vector<std::vector<std::vector<int>>> _entity_dofs;
+  std::vector<std::vector<std::vector<int>>> _entity_closure_dofs;
 
   // Quadrature points of a quadrature element (0 dimensional array for
   // all elements except quadrature elements)
diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp
index 3fd0a2673f7..dad04f63ba3 100644
--- a/cpp/dolfinx/fem/utils.cpp
+++ b/cpp/dolfinx/fem/utils.cpp
@@ -27,73 +27,6 @@
 
 using namespace dolfinx;
 
-//-----------------------------------------------------------------------------
-fem::ElementDofLayout
-fem::create_element_dof_layout(const ufcx_dofmap& dofmap,
-                               const mesh::CellType cell_type,
-                               const std::vector<int>& parent_map)
-{
-  const int element_block_size = dofmap.block_size;
-
-  // Fill entity dof indices
-  const int tdim = mesh::cell_dim(cell_type);
-  std::vector<std::vector<std::vector<int>>> entity_dofs(tdim + 1);
-  std::vector<std::vector<std::vector<int>>> entity_closure_dofs(tdim + 1);
-  {
-    int* offset0 = dofmap.entity_dof_offsets;
-    int* offset1 = dofmap.entity_closure_dof_offsets;
-    for (int d = 0; d <= tdim; ++d)
-    {
-      int num_entities = mesh::cell_num_entities(cell_type, d);
-      entity_dofs[d].resize(num_entities);
-      entity_closure_dofs[d].resize(num_entities);
-      for (int i = 0; i < num_entities; ++i)
-      {
-        std::copy(dofmap.entity_dofs + *offset0,
-                  dofmap.entity_dofs + *(offset0 + 1),
-                  std::back_inserter(entity_dofs[d][i]));
-        std::copy(dofmap.entity_closure_dofs + *offset1,
-                  dofmap.entity_closure_dofs + *(offset1 + 1),
-                  std::back_inserter(entity_closure_dofs[d][i]));
-        ++offset0;
-        ++offset1;
-      }
-    }
-  }
-
-  // TODO: UFC dofmaps just use simple offset for each field but this
-  // could be different for custom dofmaps. This data should come
-  // directly from the UFC interface in place of the implicit
-  // assumption.
-
-  // Create UFC subdofmaps and compute offset
-  std::vector<int> offsets(1, 0);
-  std::vector<ElementDofLayout> sub_doflayout;
-  for (int i = 0; i < dofmap.num_sub_dofmaps; ++i)
-  {
-    ufcx_dofmap* ufcx_sub_dofmap = dofmap.sub_dofmaps[i];
-    if (element_block_size == 1)
-    {
-      offsets.push_back(offsets.back()
-                        + ufcx_sub_dofmap->num_element_support_dofs
-                              * ufcx_sub_dofmap->block_size);
-    }
-    else
-      offsets.push_back(offsets.back() + 1);
-
-    std::vector<int> parent_map_sub(ufcx_sub_dofmap->num_element_support_dofs
-                                    * ufcx_sub_dofmap->block_size);
-    for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
-      parent_map_sub[j] = offsets[i] + element_block_size * j;
-    sub_doflayout.push_back(
-        create_element_dof_layout(*ufcx_sub_dofmap, cell_type, parent_map_sub));
-  }
-
-  // Check for "block structure". This should ultimately be replaced,
-  // but keep for now to mimic existing code.
-  return ElementDofLayout(element_block_size, entity_dofs, entity_closure_dofs,
-                          parent_map, sub_doflayout);
-}
 //-----------------------------------------------------------------------------
 fem::DofMap fem::create_dofmap(
     MPI_Comm comm, const ElementDofLayout& layout, mesh::Topology& topology,
diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h
index 69ce8d1c3bb..714de9c0804 100644
--- a/cpp/dolfinx/fem/utils.h
+++ b/cpp/dolfinx/fem/utils.h
@@ -239,11 +239,40 @@ la::SparsityPattern create_sparsity_pattern(const Form<T, U>& a)
   return pattern;
 }
 
-/// Create an ElementDofLayout from a ufcx_dofmap
-ElementDofLayout create_element_dof_layout(const ufcx_dofmap& dofmap,
-                                           const mesh::CellType cell_type,
+/// Create an ElementDofLayout from a FiniteElement
+template <std::floating_point T>
+ElementDofLayout create_element_dof_layout(const fem::FiniteElement<T>& element,
                                            const std::vector<int>& parent_map
-                                           = {});
+                                           = {})
+{
+  // Create subdofmaps and compute offset
+  std::vector<int> offsets(1, 0);
+  std::vector<dolfinx::fem::ElementDofLayout> sub_doflayout;
+  int bs = element.block_size();
+  for (int i = 0; i < element.num_sub_elements(); ++i)
+  {
+    // The ith sub-element. For mixed elements this is subelements()[i]. For
+    // blocked elements, the sub-element will always be the same, so we'll use
+    // sub_elements()[0]
+    std::shared_ptr<const fem::FiniteElement<T>> sub_e
+        = element.sub_elements()[bs > 1 ? 0 : i];
+
+    // In a mixed element DOFs are ordered element by element, so the offset to
+    // the next sub-element is sub_e->space_dimension(). Blocked elements use
+    // xxyyzz ordering, so the offset to the next sub-element is 1
+
+    std::vector<int> parent_map_sub(sub_e->space_dimension(), offsets.back());
+    for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
+      parent_map_sub[j] += bs * j;
+    offsets.push_back(offsets.back() + (bs > 1 ? 1 : sub_e->space_dimension()));
+    sub_doflayout.push_back(
+        dolfinx::fem::create_element_dof_layout(*sub_e, parent_map_sub));
+  }
+
+  return ElementDofLayout(bs, element.entity_dofs(),
+                          element.entity_closure_dofs(), parent_map,
+                          sub_doflayout);
+}
 
 /// @brief Create a dof map on mesh
 /// @param[in] comm MPI communicator
@@ -335,14 +364,7 @@ Form<T, U> create_form_factory(
   for (std::size_t i = 0; i < spaces.size(); ++i)
   {
     assert(spaces[i]->element());
-    ufcx_finite_element* ufcx_element = ufcx_form.finite_elements[i];
-    assert(ufcx_element);
-    if (std::string(ufcx_element->signature)
-        != spaces[i]->element()->signature())
-    {
-      throw std::runtime_error(
-          "Cannot create form. Wrong type of function space for argument.");
-    }
+    // TODO: Hash check of element
   }
 #endif
 
@@ -760,68 +782,6 @@ FunctionSpace<T> create_functionspace(
   return FunctionSpace(mesh, _e, dofmap, _value_shape);
 }
 
-/// @brief Create a FunctionSpace from UFC data.
-/// @param[in] fptr Pointer to a ufcx_function_space_create function.
-/// @param[in] function_name Name of a function whose function space is to
-/// create. Function name is the name of the Python variable for
-/// ufl.Coefficient, ufl.TrialFunction or ufl.TestFunction as defined in
-/// the UFL file.
-/// @param[in] mesh Mesh
-/// @param[in] reorder_fn 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(
-    ufcx_function_space* (*fptr)(const char*), const std::string& function_name,
-    std::shared_ptr<mesh::Mesh<T>> mesh,
-    std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
-        reorder_fn
-    = nullptr)
-{
-  ufcx_function_space* space = fptr(function_name.c_str());
-  if (!space)
-  {
-    throw std::runtime_error(
-        "Could not create UFC function space with function name "
-        + function_name);
-  }
-
-  ufcx_finite_element* ufcx_element = space->finite_element;
-  assert(ufcx_element);
-  std::vector<std::size_t> value_shape(space->value_shape,
-                                       space->value_shape + space->value_rank);
-
-  const auto& geometry = mesh->geometry();
-  auto& cmap = geometry.cmap();
-  if (space->geometry_degree != cmap.degree()
-      or static_cast<basix::cell::type>(space->geometry_basix_cell)
-             != mesh::cell_type_to_basix_type(cmap.cell_shape())
-      or static_cast<basix::element::lagrange_variant>(
-             space->geometry_basix_variant)
-             != cmap.variant())
-  {
-    throw std::runtime_error("UFL mesh and CoordinateElement do not match.");
-  }
-
-  auto element = std::make_shared<FiniteElement<T>>(*ufcx_element);
-  assert(element);
-  ufcx_dofmap* ufcx_map = space->dofmap;
-  assert(ufcx_map);
-  const auto topology = mesh->topology();
-  assert(topology);
-  ElementDofLayout layout
-      = create_element_dof_layout(*ufcx_map, topology->cell_type());
-
-  std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv;
-  if (element->needs_dof_permutations())
-    permute_inv = element->dof_permutation_fn(true, true);
-  return FunctionSpace(
-      mesh, element,
-      std::make_shared<DofMap>(create_dofmap(mesh->comm(), layout, *topology,
-                                             permute_inv, reorder_fn)),
-      value_shape);
-}
-
 /// @private
 namespace impl
 {
diff --git a/cpp/test/matrix.cpp b/cpp/test/matrix.cpp
index c406903c13a..c2c00b8c84e 100644
--- a/cpp/test/matrix.cpp
+++ b/cpp/test/matrix.cpp
@@ -111,8 +111,13 @@ la::MatrixCSR<double> create_operator(MPI_Comm comm)
   auto mesh = std::make_shared<mesh::Mesh<double>>(
       mesh::create_box(comm, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {12, 12, 12},
                        mesh::CellType::tetrahedron, part));
+  auto element = basix::create_element<double>(
+      basix::element::family::P, basix::cell::type::tetrahedron, 2,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
   auto V = std::make_shared<fem::FunctionSpace<double>>(
-      fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
+      fem::create_functionspace(mesh, element, {}));
 
   // Prepare and set Constants for the bilinear form
   auto kappa = std::make_shared<fem::Constant<double>>(2.0);
@@ -144,8 +149,13 @@ la::MatrixCSR<double> create_operator(MPI_Comm comm)
       mesh::create_box(comm, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {12, 12, 12},
                        mesh::CellType::tetrahedron, part));
 
+  auto element = basix::create_element<double>(
+      basix::element::family::P, basix::cell::type::tetrahedron, 2,
+      basix::element::lagrange_variant::unset,
+      basix::element::dpc_variant::unset, false);
+
   auto V = std::make_shared<fem::FunctionSpace<double>>(
-      fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
+      fem::create_functionspace(mesh, element, {}));
 
   // Prepare and set Constants for the bilinear form
   auto kappa = std::make_shared<fem::Constant<double>>(2.0);
diff --git a/python/demo/demo_axis/demo_axis.py b/python/demo/demo_axis/demo_axis.py
index defe96b7cd2..d9530190569 100644
--- a/python/demo/demo_axis/demo_axis.py
+++ b/python/demo/demo_axis/demo_axis.py
@@ -25,7 +25,7 @@
 
 import ufl
 from basix.ufl import element, mixed_element
-from dolfinx import default_scalar_type, fem, io, mesh, plot
+from dolfinx import default_real_type, default_scalar_type, fem, io, mesh, plot
 from dolfinx.fem.petsc import LinearProblem
 
 try:
@@ -403,8 +403,8 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # will use Lagrange elements:
 
 degree = 3
-curl_el = element("N1curl", msh.basix_cell(), degree)
-lagr_el = element("Lagrange", msh.basix_cell(), degree)
+curl_el = element("N1curl", msh.basix_cell(), degree, dtype=default_real_type)
+lagr_el = element("Lagrange", msh.basix_cell(), degree, dtype=default_real_type)
 V = fem.functionspace(msh, mixed_element([curl_el, lagr_el]))
 
 # The integration domains of our problem are the following:
@@ -678,7 +678,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # assert err_ext < 0.01
 
 if has_vtx:
-    v_dg_el = element("DG", msh.basix_cell(), degree, shape=(3,))
+    v_dg_el = element("DG", msh.basix_cell(), degree, shape=(3,), dtype=default_real_type)
     W = fem.functionspace(msh, v_dg_el)
     Es_dg = fem.Function(W)
     Es_expr = fem.Expression(Esh, W.element.interpolation_points())
diff --git a/python/demo/demo_cahn-hilliard.py b/python/demo/demo_cahn-hilliard.py
index 0951b2c989e..e75a06211fb 100644
--- a/python/demo/demo_cahn-hilliard.py
+++ b/python/demo/demo_cahn-hilliard.py
@@ -163,7 +163,7 @@
 # using a pair of linear Lagrange elements.
 
 msh = create_unit_square(MPI.COMM_WORLD, 96, 96, CellType.triangle)
-P1 = element("Lagrange", msh.basix_cell(), 1)
+P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
 ME = functionspace(msh, mixed_element([P1, P1]))
 
 # Trial and test functions of the space `ME` are now defined:
diff --git a/python/demo/demo_lagrange_variants.py b/python/demo/demo_lagrange_variants.py
index 2fedd7821f5..b30bfba62de 100644
--- a/python/demo/demo_lagrange_variants.py
+++ b/python/demo/demo_lagrange_variants.py
@@ -25,7 +25,7 @@
 import basix
 import basix.ufl
 import ufl  # type: ignore
-from dolfinx import fem, mesh
+from dolfinx import default_real_type, fem, mesh
 from ufl import dx
 
 # -
@@ -46,7 +46,11 @@
 
 # +
 element = basix.ufl.element(
-    basix.ElementFamily.P, basix.CellType.interval, 10, basix.LagrangeVariant.equispaced
+    basix.ElementFamily.P,
+    basix.CellType.interval,
+    10,
+    basix.LagrangeVariant.equispaced,
+    dtype=default_real_type,
 )
 lattice = basix.create_lattice(basix.CellType.interval, 200, basix.LatticeType.equispaced, True)
 values = element.tabulate(0, lattice)[0, :, :]
@@ -75,7 +79,11 @@
 
 # +
 element = basix.ufl.element(
-    basix.ElementFamily.P, basix.CellType.interval, 10, basix.LagrangeVariant.gll_warped
+    basix.ElementFamily.P,
+    basix.CellType.interval,
+    10,
+    basix.LagrangeVariant.gll_warped,
+    dtype=default_real_type,
 )
 values = element.tabulate(0, lattice)[0, :, :]
 if MPI.COMM_WORLD.size == 1:  # Skip this plotting in parallel
@@ -118,7 +126,9 @@ def saw_tooth(x):
 u_exact = saw_tooth(x[0])
 
 for variant in [basix.LagrangeVariant.equispaced, basix.LagrangeVariant.gll_warped]:
-    ufl_element = basix.ufl.element(basix.ElementFamily.P, basix.CellType.interval, 10, variant)
+    ufl_element = basix.ufl.element(
+        basix.ElementFamily.P, basix.CellType.interval, 10, variant, dtype=default_real_type
+    )
     V = fem.functionspace(msh, ufl_element)
     uh = fem.Function(V)
     uh.interpolate(lambda x: saw_tooth(x[0]))
@@ -157,7 +167,9 @@ def saw_tooth(x):
 
 # +
 for variant in [basix.LagrangeVariant.equispaced, basix.LagrangeVariant.gll_warped]:
-    ufl_element = basix.ufl.element(basix.ElementFamily.P, basix.CellType.interval, 10, variant)
+    ufl_element = basix.ufl.element(
+        basix.ElementFamily.P, basix.CellType.interval, 10, variant, dtype=default_real_type
+    )
     V = fem.functionspace(msh, ufl_element)
     uh = fem.Function(V)
     uh.interpolate(lambda x: saw_tooth(x[0]))
diff --git a/python/demo/demo_mixed-poisson.py b/python/demo/demo_mixed-poisson.py
index 23daf1a14cb..413bbbbd151 100644
--- a/python/demo/demo_mixed-poisson.py
+++ b/python/demo/demo_mixed-poisson.py
@@ -99,15 +99,15 @@
 import numpy as np
 
 from basix.ufl import element, mixed_element
-from dolfinx import fem, io, mesh
+from dolfinx import default_real_type, fem, io, mesh
 from dolfinx.fem.petsc import LinearProblem
 from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner
 
 msh = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
 
 k = 1
-Q_el = element("BDMCF", msh.basix_cell(), k)
-P_el = element("DG", msh.basix_cell(), k - 1)
+Q_el = element("BDMCF", msh.basix_cell(), k, dtype=default_real_type)
+P_el = element("DG", msh.basix_cell(), k - 1, dtype=default_real_type)
 V_el = mixed_element([Q_el, P_el])
 V = fem.functionspace(msh, V_el)
 
diff --git a/python/demo/demo_pml/demo_pml.py b/python/demo/demo_pml/demo_pml.py
index 4d78764ce97..6cb5dd36f21 100644
--- a/python/demo/demo_pml/demo_pml.py
+++ b/python/demo/demo_pml/demo_pml.py
@@ -26,7 +26,7 @@
 
 import ufl
 from basix.ufl import element
-from dolfinx import default_scalar_type, fem, mesh, plot
+from dolfinx import default_real_type, default_scalar_type, fem, mesh, plot
 from dolfinx.fem.petsc import LinearProblem
 from dolfinx.io import VTXWriter, gmshio
 
@@ -243,7 +243,7 @@ def pml_coordinates(x: ufl.indexed.Indexed, alpha: float, k0: complex, l_dom: fl
 # element to represent the electric field:
 
 degree = 3
-curl_el = element("N1curl", msh.basix_cell(), degree)
+curl_el = element("N1curl", msh.basix_cell(), degree, dtype=default_real_type)
 V = fem.functionspace(msh, curl_el)
 
 # Next, we interpolate $\mathbf{E}_b$ into the function space $V$,
diff --git a/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py b/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py
index d7aea7a937b..8c274626d75 100644
--- a/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py
+++ b/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py
@@ -39,7 +39,7 @@
 
 import ufl
 from basix.ufl import element
-from dolfinx import default_scalar_type, fem, io, plot
+from dolfinx import default_real_type, default_scalar_type, fem, io, plot
 from dolfinx.fem.petsc import LinearProblem
 
 try:
@@ -283,7 +283,7 @@ def curl_2d(f: fem.Function):
 # represent the electric field
 
 degree = 3
-curl_el = element("N1curl", domain.basix_cell(), degree)
+curl_el = element("N1curl", domain.basix_cell(), degree, dtype=default_real_type)
 V = fem.functionspace(domain, curl_el)
 
 # Next, we can interpolate $\mathbf{E}_b$ into the function space $V$:
diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py
index 33ced9fc832..ecc71ee6aa0 100644
--- a/python/demo/demo_stokes.py
+++ b/python/demo/demo_stokes.py
@@ -92,7 +92,7 @@
 
 import ufl
 from basix.ufl import element, mixed_element
-from dolfinx import fem, la
+from dolfinx import default_real_type, fem, la
 from dolfinx.fem import (
     Constant,
     Function,
@@ -141,8 +141,8 @@ def lid_velocity_expression(x):
 # piecewise linear basis (scalar).
 
 
-P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
-P1 = element("Lagrange", msh.basix_cell(), 1)
+P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
+P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
 V, Q = functionspace(msh, P2), functionspace(msh, P1)
 
 # Boundary conditions for the velocity field are defined:
@@ -282,7 +282,9 @@ def nested_iterative_solver():
     # `scatter_forward`.
     with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
         u.x.scatter_forward()
-        P1 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
+        P1 = element(
+            "Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,), dtype=default_real_type
+        )
         u1 = Function(functionspace(msh, P1))
         u1.interpolate(u)
         ufile_xdmf.write_mesh(msh)
diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 7aee1db32a7..988fcdb5bfe 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -29,7 +29,7 @@
 
 import basix
 import basix.ufl
-from dolfinx import fem, mesh
+from dolfinx import default_real_type, fem, mesh
 from dolfinx.fem.petsc import LinearProblem
 from ufl import SpatialCoordinate, TestFunction, TrialFunction, cos, div, dx, grad, inner, sin
 
@@ -146,6 +146,7 @@
     False,
     1,
     2,
+    dtype=default_real_type,
 )
 
 # ## Creating higher degree TNT elements
@@ -228,6 +229,7 @@ def create_tnt_quad(degree):
         False,
         degree,
         degree + 1,
+        dtype=default_real_type,
     )
 
 
diff --git a/python/demo/demo_waveguide/demo_half_loaded_waveguide.py b/python/demo/demo_waveguide/demo_half_loaded_waveguide.py
index 39947b3de94..91ae040e0b0 100644
--- a/python/demo/demo_waveguide/demo_half_loaded_waveguide.py
+++ b/python/demo/demo_waveguide/demo_half_loaded_waveguide.py
@@ -48,7 +48,7 @@
 
 import ufl
 from basix.ufl import element, mixed_element
-from dolfinx import default_scalar_type, fem, io, plot
+from dolfinx import default_real_type, default_scalar_type, fem, io, plot
 from dolfinx.fem.petsc import assemble_matrix
 from dolfinx.mesh import CellType, create_rectangle, exterior_facet_indices, locate_entities
 
@@ -193,8 +193,8 @@ def Omega_v(x):
 # `mixed_element`:
 
 degree = 1
-RTCE = element("RTCE", msh.basix_cell(), degree)
-Q = element("Lagrange", msh.basix_cell(), degree)
+RTCE = element("RTCE", msh.basix_cell(), degree, dtype=default_real_type)
+Q = element("Lagrange", msh.basix_cell(), degree, dtype=default_real_type)
 V = fem.functionspace(msh, mixed_element([RTCE, Q]))
 
 # Now we can define our weak form:
diff --git a/python/dolfinx/fem/dofmap.py b/python/dolfinx/fem/dofmap.py
index 5cdd2399ba4..21ccbc1c43a 100644
--- a/python/dolfinx/fem/dofmap.py
+++ b/python/dolfinx/fem/dofmap.py
@@ -11,7 +11,7 @@ class DofMap:
     """Degree-of-freedom map
 
     This class handles the mapping of degrees of freedom. It builds
-    a dof map based on a ufcx_dofmap on a specific mesh.
+    a dof map based on a FiniteElement on a specific mesh.
     """
 
     _cpp_object: _cpp.fem.DofMap
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 25eecbbe7f2..fc375eee5d9 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -609,6 +609,38 @@ class ElementMetaData(typing.NamedTuple):
     symmetry: typing.Optional[bool] = None
 
 
+def _create_dolfinx_element(
+    comm: _MPI.Intracomm,
+    cell_type: _cpp.mesh.CellType,
+    ufl_e: ufl.FiniteElementBase,
+    dtype: np.dtype,
+) -> typing.Union[_cpp.fem.FiniteElement_float32, _cpp.fem.FiniteElement_float64]:
+    """Create a DOLFINx element from a basix.ufl element."""
+    if np.issubdtype(dtype, np.float32):
+        CppElement = _cpp.fem.FiniteElement_float32
+    elif np.issubdtype(dtype, np.float64):
+        CppElement = _cpp.fem.FiniteElement_float64
+    else:
+        raise ValueError(f"Unsupported dtype: {dtype}")
+
+    if ufl_e.is_mixed:
+        elements = [
+            _create_dolfinx_element(
+                comm,
+                cell_type,
+                e,
+                dtype,
+            )
+            for e in ufl_e.sub_elements
+        ]
+        return CppElement(elements)
+    elif ufl_e.is_quadrature:
+        return CppElement(cell_type, ufl_e.custom_quadrature()[0], ufl_e.block_size)
+    else:
+        basix_e = ufl_e.basix_element._e
+        return CppElement(basix_e, ufl_e.block_size, ufl_e.is_symmetric)
+
+
 def functionspace(
     mesh: Mesh,
     element: typing.Union[ufl.FiniteElementBase, ElementMetaData, tuple[str, int, tuple, bool]],
@@ -628,6 +660,7 @@ def functionspace(
 
     """
     # Create UFL element
+    dtype = mesh.geometry.x.dtype
     try:
         e = ElementMetaData(*element)
         ufl_e = basix.ufl.element(
@@ -636,6 +669,7 @@ def functionspace(
             e.degree,
             shape=e.shape,
             symmetry=e.symmetry,
+            dtype=dtype,
         )
     except TypeError:
         ufl_e = element  # type: ignore
@@ -648,30 +682,19 @@ def functionspace(
     value_shape = ufl_space.value_shape
 
     # Compile dofmap and element and create DOLFINx objects
-    dtype = mesh.geometry.x.dtype
     if form_compiler_options is None:
         form_compiler_options = dict()
     form_compiler_options["scalar_type"] = dtype
-    (ufcx_element, ufcx_dofmap), module, code = jit.ffcx_jit(
+
+    cpp_element = _create_dolfinx_element(
         mesh.comm,
+        mesh.topology.cell_type,
         ufl_e,
-        form_compiler_options=form_compiler_options,
-        jit_options=jit_options,
+        dtype,
     )
 
-    ffi = module.ffi
-    if np.issubdtype(dtype, np.float32):
-        cpp_element = _cpp.fem.FiniteElement_float32(
-            ffi.cast("uintptr_t", ffi.addressof(ufcx_element))
-        )
-    elif np.issubdtype(dtype, np.float64):
-        cpp_element = _cpp.fem.FiniteElement_float64(
-            ffi.cast("uintptr_t", ffi.addressof(ufcx_element))
-        )
-
     cpp_dofmap = _cpp.fem.create_dofmap(
         mesh.comm,
-        ffi.cast("uintptr_t", ffi.addressof(ufcx_dofmap)),
         mesh.topology,
         cpp_element,
     )
diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py
index 06627ca7f43..d1f6ed4b756 100644
--- a/python/dolfinx/fem/petsc.py
+++ b/python/dolfinx/fem/petsc.py
@@ -760,11 +760,17 @@ def __init__(
                                                                        "mumps"})
         """
         self._a = _create_form(
-            a, form_compiler_options=form_compiler_options, jit_options=jit_options
+            a,
+            dtype=PETSc.ScalarType,
+            form_compiler_options=form_compiler_options,
+            jit_options=jit_options,
         )
         self._A = create_matrix(self._a)
         self._L = _create_form(
-            L, form_compiler_options=form_compiler_options, jit_options=jit_options
+            L,
+            dtype=PETSc.ScalarType,
+            form_compiler_options=form_compiler_options,
+            jit_options=jit_options,
         )
         self._b = create_vector(self._L)
 
diff --git a/python/dolfinx/jit.py b/python/dolfinx/jit.py
index 0e031b96189..6f4d1c80532 100644
--- a/python/dolfinx/jit.py
+++ b/python/dolfinx/jit.py
@@ -199,8 +199,6 @@ def ffcx_jit(
     # Switch on type and compile, returning cffi object
     if isinstance(ufl_object, ufl.Form):
         r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
-    elif isinstance(ufl_object, ufl.AbstractFiniteElement):
-        r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)
     elif isinstance(ufl_object, ufl.Mesh):
         r = ffcx.codegeneration.jit.compile_coordinate_maps([ufl_object], options=p_ffcx, **p_jit)
     elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 1bb5bf62740..6879a1ab64e 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -107,13 +107,32 @@ void declare_function_space(nb::module_& m, std::string type)
         .def(
             "__init__",
             [](dolfinx::fem::FiniteElement<T>* self,
-               std::uintptr_t ufcx_element)
+               basix::FiniteElement<T>& element, std::size_t block_size,
+               bool symmetric) {
+              new (self) dolfinx::fem::FiniteElement<T>(element, block_size,
+                                                        symmetric);
+            },
+            nb::arg("element"), nb::arg("block_size"), nb::arg("symmetric"))
+        .def(
+            "__init__",
+            [](dolfinx::fem::FiniteElement<T>* self,
+               std::vector<
+                   std::shared_ptr<const dolfinx::fem::FiniteElement<T>>>
+                   elements)
+            { new (self) dolfinx::fem::FiniteElement<T>(elements); },
+            nb::arg("elements"))
+        .def(
+            "__init__",
+            [](dolfinx::fem::FiniteElement<T>* self, mesh::CellType cell_type,
+               nb::ndarray<T, nb::ndim<2>, nb::numpy> points,
+               std::size_t block_size)
             {
-              ufcx_finite_element* p
-                  = reinterpret_cast<ufcx_finite_element*>(ufcx_element);
-              new (self) dolfinx::fem::FiniteElement<T>(*p);
+              std::span<T> pdata(points.data(), points.size());
+              new (self) dolfinx::fem::FiniteElement<T>(
+                  cell_type, pdata, {points.shape(0), points.shape(1)},
+                  block_size);
             },
-            nb::arg("ufcx_element"))
+            nb::arg("cell_type"), nb::arg("points"), nb::arg("block_size"))
         .def("__eq__", &dolfinx::fem::FiniteElement<T>::operator==)
         .def_prop_ro("dtype", [](const dolfinx::fem::FiniteElement<T>&)
                      { return dolfinx_wrappers::numpy_dtype<T>(); })
@@ -898,16 +917,21 @@ void declare_cmap(nb::module_& m, std::string type)
 template <typename T>
 void declare_real_functions(nb::module_& m)
 {
+  m.def(
+      "create_element_dof_layout",
+      [](const dolfinx::fem::FiniteElement<T>& element,
+         const std::vector<int>& parent_map)
+      { return dolfinx::fem::create_element_dof_layout(element, parent_map); },
+      nb::arg("element"), nb::arg("parent_map"),
+      "Create ElementDofLayout object from a ufc dofmap.");
   m.def(
       "create_dofmap",
-      [](const dolfinx_wrappers::MPICommWrapper comm, std::uintptr_t dofmap,
+      [](const dolfinx_wrappers::MPICommWrapper comm,
          dolfinx::mesh::Topology& topology,
          const dolfinx::fem::FiniteElement<T>& element)
       {
-        ufcx_dofmap* p = reinterpret_cast<ufcx_dofmap*>(dofmap);
-        assert(p);
         dolfinx::fem::ElementDofLayout layout
-            = dolfinx::fem::create_element_dof_layout(*p, topology.cell_type());
+            = dolfinx::fem::create_element_dof_layout(element);
 
         std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv
             = nullptr;
@@ -916,33 +940,30 @@ void declare_real_functions(nb::module_& m)
         return dolfinx::fem::create_dofmap(comm.get(), layout, topology,
                                            permute_inv, nullptr);
       },
-      nb::arg("comm"), nb::arg("dofmap"), nb::arg("topology"),
-      nb::arg("element"),
-      "Create DofMap object from a pointer to ufcx_dofmap.");
-
+      nb::arg("comm"), nb::arg("topology"), nb::arg("element"),
+      "Create DofMap object from an element.");
   m.def(
       "create_dofmaps",
       [](const dolfinx_wrappers::MPICommWrapper comm,
-         std::vector<std::uintptr_t> ufcx_dofmaps,
-         dolfinx::mesh::Topology& topology)
+         dolfinx::mesh::Topology& topology,
+         std::vector<std::shared_ptr<const dolfinx::fem::FiniteElement<T>>>
+             elements)
       {
         std::vector<dolfinx::fem::ElementDofLayout> layouts;
         int D = topology.dim();
-        assert(ufcx_dofmaps.size() == topology.entity_types(D).size());
-        for (std::size_t i = 0; i < ufcx_dofmaps.size(); ++i)
+        assert(elements.size() == topology.entity_types(D).size());
+        for (std::size_t i = 0; i < elements.size(); ++i)
         {
-          ufcx_dofmap* p = reinterpret_cast<ufcx_dofmap*>(ufcx_dofmaps[i]);
-          assert(p);
-          layouts.push_back(dolfinx::fem::create_element_dof_layout(
-              *p, topology.entity_types(D)[i]));
+          layouts.push_back(
+              dolfinx::fem::create_element_dof_layout(*elements[i]));
         }
 
         return dolfinx::fem::create_dofmaps(comm.get(), layouts, topology,
                                             nullptr, nullptr);
       },
-      nb::arg("comm"), nb::arg("dofmap"), nb::arg("topology"),
+      nb::arg("comm"), nb::arg("topology"), nb::arg("elements"),
       "Create DofMap objects on a mixed topology mesh from pointers to "
-      "ufcx_dofmaps.");
+      "FiniteElements.");
 
   m.def(
       "locate_dofs_topological",
@@ -1104,17 +1125,6 @@ void fem(nb::module_& m)
   declare_cmap<float>(m, "float32");
   declare_cmap<double>(m, "float64");
 
-  m.def(
-      "create_element_dof_layout",
-      [](std::uintptr_t dofmap, const dolfinx::mesh::CellType cell_type,
-         const std::vector<int>& parent_map)
-      {
-        ufcx_dofmap* p = reinterpret_cast<ufcx_dofmap*>(dofmap);
-        return dolfinx::fem::create_element_dof_layout(*p, cell_type,
-                                                       parent_map);
-      },
-      nb::arg("dofmap"), nb::arg("cell_type"), nb::arg("parent_map"),
-      "Create ElementDofLayout object from a ufc dofmap.");
   m.def(
       "build_dofmap",
       [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology,
diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py
index 4b248060528..5b9ec719278 100644
--- a/python/test/unit/fem/test_assembler.py
+++ b/python/test/unit/fem/test_assembler.py
@@ -283,9 +283,9 @@ def test_matrix_assembly_block(mode):
     structures"""
     mesh = create_unit_square(MPI.COMM_WORLD, 4, 8, ghost_mode=mode)
     p0, p1 = 1, 2
-    P0 = element("Lagrange", mesh.basix_cell(), p0)
-    P1 = element("Lagrange", mesh.basix_cell(), p1)
-    P2 = element("Lagrange", mesh.basix_cell(), p0)
+    P0 = element("Lagrange", mesh.basix_cell(), p0, dtype=default_real_type)
+    P1 = element("Lagrange", mesh.basix_cell(), p1, dtype=default_real_type)
+    P2 = element("Lagrange", mesh.basix_cell(), p0, dtype=default_real_type)
     V0 = functionspace(mesh, P0)
     V1 = functionspace(mesh, P1)
     V2 = functionspace(mesh, P2)
@@ -410,7 +410,7 @@ def test_assembly_solve_block(mode):
     """Solve a two-field mass-matrix like problem with block matrix approaches
     and test that solution is the same"""
     mesh = create_unit_square(MPI.COMM_WORLD, 32, 31, ghost_mode=mode)
-    P = element("Lagrange", mesh.basix_cell(), 1)
+    P = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
     V0 = functionspace(mesh, P)
     V1 = V0.clone()
 
@@ -674,8 +674,10 @@ def blocked_solve():
 
     def monolithic_solve():
         """Monolithic (interleaved) solver"""
-        P2_el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
-        P1_el = element("Lagrange", mesh.basix_cell(), 1)
+        P2_el = element(
+            "Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), dtype=default_real_type
+        )
+        P1_el = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
         TH = mixed_element([P2_el, P1_el])
         W = functionspace(mesh, TH)
         (u, p) = ufl.TrialFunctions(W)
@@ -994,13 +996,17 @@ def test_coefficents_non_constant():
     F = form((ufl.inner(u, v) - ufl.inner(x[0] * x[1] ** 2, v)) * dx)
     b0 = petsc_assemble_vector(F)
     b0.assemble()
-    assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7)
+    assert np.linalg.norm(b0.array) == pytest.approx(
+        0.0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
+    )
 
     # -- Exterior facet integral vector
     F = form((ufl.inner(u, v) - ufl.inner(x[0] * x[1] ** 2, v)) * ds)
     b0 = petsc_assemble_vector(F)
     b0.assemble()
-    assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7)
+    assert np.linalg.norm(b0.array) == pytest.approx(
+        0.0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
+    )
 
     # -- Interior facet integral vector
     V = functionspace(mesh, ("DG", 3))  # degree 3 so that interpolation is exact
@@ -1019,7 +1025,9 @@ def test_coefficents_non_constant():
     F = form(F)
     b0 = petsc_assemble_vector(F)
     b0.assemble()
-    assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7)
+    assert np.linalg.norm(b0.array) == pytest.approx(
+        0.0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
+    )
 
     b0.destroy()
 
diff --git a/python/test/unit/fem/test_bcs.py b/python/test/unit/fem/test_bcs.py
index 1763b34f465..15108674258 100644
--- a/python/test/unit/fem/test_bcs.py
+++ b/python/test/unit/fem/test_bcs.py
@@ -36,8 +36,8 @@ def test_locate_dofs_geometrical():
     spaces, returns the correct degrees of freedom in each space"""
     mesh = create_unit_square(MPI.COMM_WORLD, 4, 8)
     p0, p1 = 1, 2
-    P0 = element("Lagrange", mesh.basix_cell(), p0)
-    P1 = element("Lagrange", mesh.basix_cell(), p1)
+    P0 = element("Lagrange", mesh.basix_cell(), p0, dtype=default_real_type)
+    P1 = element("Lagrange", mesh.basix_cell(), p1, dtype=default_real_type)
 
     W = functionspace(mesh, mixed_element([P0, P1]))
     V = W.sub(0).collapse()[0]
@@ -292,7 +292,10 @@ def test_mixed_constant_bc(mesh_factory):
         mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
     )
     TH = mixed_element(
-        [element("Lagrange", mesh.basix_cell(), 1), element("Lagrange", mesh.basix_cell(), 2)]
+        [
+            element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type),
+            element("Lagrange", mesh.basix_cell(), 2, dtype=default_real_type),
+        ]
     )
     W = functionspace(mesh, TH)
     u = Function(W)
@@ -333,8 +336,14 @@ def test_mixed_blocked_constant():
 
     TH = mixed_element(
         [
-            element("Lagrange", mesh.basix_cell(), 1),
-            element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)),
+            element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type),
+            element(
+                "Lagrange",
+                mesh.basix_cell(),
+                2,
+                shape=(mesh.geometry.dim,),
+                dtype=default_real_type,
+            ),
         ]
     )
     W = functionspace(mesh, TH)
diff --git a/python/test/unit/fem/test_complex_assembler.py b/python/test/unit/fem/test_complex_assembler.py
index fd89dcdc118..75b5512dd55 100644
--- a/python/test/unit/fem/test_complex_assembler.py
+++ b/python/test/unit/fem/test_complex_assembler.py
@@ -24,7 +24,7 @@ def test_complex_assembly(complex_dtype):
 
     real_dtype = np.real(complex_dtype(1.0)).dtype
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, dtype=real_dtype)
-    P2 = element("Lagrange", mesh.basix_cell(), 2)
+    P2 = element("Lagrange", mesh.basix_cell(), 2, dtype=real_dtype)
     V = functionspace(mesh, P2)
     u = ufl.TrialFunction(V)
     v = ufl.TestFunction(V)
@@ -78,7 +78,7 @@ def test_complex_assembly_solve(complex_dtype, cg_solver):
     degree = 3
     real_dtype = np.real(complex_dtype(1.0)).dtype
     mesh = create_unit_square(MPI.COMM_WORLD, 20, 20, dtype=real_dtype)
-    P = element("Lagrange", mesh.basix_cell(), degree)
+    P = element("Lagrange", mesh.basix_cell(), degree, dtype=real_dtype)
     V = functionspace(mesh, P)
 
     x = ufl.SpatialCoordinate(mesh)
diff --git a/python/test/unit/fem/test_custom_basix_element.py b/python/test/unit/fem/test_custom_basix_element.py
index 9e7f8634d65..9f9a1cdbb6a 100644
--- a/python/test/unit/fem/test_custom_basix_element.py
+++ b/python/test/unit/fem/test_custom_basix_element.py
@@ -7,6 +7,7 @@
 import basix
 import basix.ufl
 import ufl
+from dolfinx import default_real_type
 from dolfinx.fem import (
     Function,
     assemble_scalar,
@@ -86,7 +87,11 @@ def run_scalar_test(V, degree):
 @pytest.mark.parametrize("degree", range(1, 6))
 def test_basix_element_wrapper(degree):
     ufl_element = basix.ufl.element(
-        basix.ElementFamily.P, basix.CellType.triangle, degree, basix.LagrangeVariant.gll_isaac
+        basix.ElementFamily.P,
+        basix.CellType.triangle,
+        degree,
+        basix.LagrangeVariant.gll_isaac,
+        dtype=default_real_type,
     )
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
     V = functionspace(mesh, ufl_element)
@@ -116,6 +121,7 @@ def test_custom_element_triangle_degree1():
         False,
         1,
         1,
+        dtype=default_real_type,
     )
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
     V = functionspace(mesh, ufl_element)
@@ -154,6 +160,7 @@ def test_custom_element_triangle_degree4():
         False,
         4,
         4,
+        dtype=default_real_type,
     )
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
     V = functionspace(mesh, ufl_element)
@@ -162,9 +169,9 @@ def test_custom_element_triangle_degree4():
 
 def test_custom_element_triangle_degree4_integral():
     pts, wts = basix.make_quadrature(basix.CellType.interval, 10)
-    tab = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, 2).tabulate(0, pts)[
-        0, :, :, 0
-    ]
+    tab = basix.create_element(
+        basix.ElementFamily.P, basix.CellType.interval, 2, dtype=default_real_type
+    ).tabulate(0, pts)[0, :, :, 0]
     wcoeffs = np.eye(15)
     x = [
         [np.array([[0.0, 0.0]]), np.array([[1.0, 0.0]]), np.array([[0.0, 1.0]])],
@@ -202,6 +209,7 @@ def test_custom_element_triangle_degree4_integral():
         False,
         4,
         4,
+        dtype=default_real_type,
     )
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
     V = functionspace(mesh, ufl_element)
@@ -246,6 +254,7 @@ def test_custom_element_quadrilateral_degree1():
         False,
         1,
         1,
+        dtype=default_real_type,
     )
     mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, CellType.quadrilateral)
     V = functionspace(mesh, ufl_element)
@@ -276,7 +285,9 @@ def test_vector_copy_degree1(cell_type, element_family):
     def func(x):
         return x[:tdim]
 
-    e1 = basix.ufl.element(element_family, getattr(basix.CellType, cell_type.name), 1)
+    e1 = basix.ufl.element(
+        element_family, getattr(basix.CellType, cell_type.name), 1, dtype=default_real_type
+    )
 
     e2 = basix.ufl.custom_element(
         e1._element.cell_type,
@@ -290,6 +301,7 @@ def func(x):
         e1._element.discontinuous,
         e1._element.embedded_subdegree,
         e1._element.embedded_superdegree,
+        dtype=default_real_type,
     )
 
     space1 = functionspace(mesh, e1)
@@ -325,7 +337,9 @@ def test_scalar_copy_degree1(cell_type, element_family):
     def func(x):
         return x[0]
 
-    e1 = basix.ufl.element(element_family, getattr(basix.CellType, cell_type.name), 1)
+    e1 = basix.ufl.element(
+        element_family, getattr(basix.CellType, cell_type.name), 1, dtype=default_real_type
+    )
     e2 = basix.ufl.custom_element(
         e1._element.cell_type,
         e1._element.value_shape,
@@ -338,6 +352,7 @@ def func(x):
         e1._element.discontinuous,
         e1._element.embedded_subdegree,
         e1._element.embedded_superdegree,
+        dtype=default_real_type,
     )
 
     space1 = functionspace(mesh, e1)
diff --git a/python/test/unit/fem/test_dofmap.py b/python/test/unit/fem/test_dofmap.py
index b3e047b3bd4..aa1ee0b4d21 100644
--- a/python/test/unit/fem/test_dofmap.py
+++ b/python/test/unit/fem/test_dofmap.py
@@ -15,6 +15,7 @@
 import dolfinx
 import ufl
 from basix.ufl import element, mixed_element
+from dolfinx import default_real_type
 from dolfinx.fem import functionspace
 from dolfinx.mesh import (
     CellType,
@@ -43,8 +44,10 @@ def mesh():
 def test_tabulate_dofs(mesh_factory):
     func, args = mesh_factory
     mesh = func(*args)
-    W0 = element("Lagrange", mesh.basix_cell(), 1)
-    W1 = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
+    W0 = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
+    W1 = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
     W = functionspace(mesh, W0 * W1)
 
     L0 = W.sub(0)
@@ -160,7 +163,7 @@ def test_block_size():
         create_unit_cube(MPI.COMM_WORLD, 4, 4, 4, CellType.hexahedron),
     ]
     for mesh in meshes:
-        P2 = element("Lagrange", mesh.basix_cell(), 2)
+        P2 = element("Lagrange", mesh.basix_cell(), 2, dtype=default_real_type)
         V = functionspace(mesh, P2)
         assert V.dofmap.bs == 1
 
@@ -180,8 +183,8 @@ def test_block_size():
 @pytest.mark.skip
 def test_block_size_real():
     mesh = create_unit_interval(MPI.COMM_WORLD, 12)
-    V = element("DG", mesh.basix_cell(), 0)
-    R = element("R", mesh.basix_cell(), 0)
+    V = element("DG", mesh.basix_cell(), 0, dtype=default_real_type)
+    R = element("R", mesh.basix_cell(), 0, dtype=default_real_type)
     X = functionspace(mesh, V * R)
     assert X.dofmap.index_map_bs == 1
 
@@ -198,8 +201,10 @@ def test_local_dimension(mesh_factory):
     func, args = mesh_factory
     mesh = func(*args)
 
-    v = element("Lagrange", mesh.basix_cell(), 1)
-    q = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
+    v = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
+    q = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
     w = v * q
 
     V = functionspace(mesh, v)
@@ -316,7 +321,9 @@ def test_readonly_view_local_to_global_unwoned(mesh):
 def test_higher_order_coordinate_map(points, celltype, order):
     """Computes physical coordinates of a cell, based on the coordinate map."""
     cells = np.array([range(len(points))])
-    domain = ufl.Mesh(element("Lagrange", celltype.name, order, shape=(points.shape[1],)))
+    domain = ufl.Mesh(
+        element("Lagrange", celltype.name, order, shape=(points.shape[1],), dtype=default_real_type)
+    )
     mesh = create_mesh(MPI.COMM_WORLD, cells, points, domain)
 
     V = functionspace(mesh, ("Lagrange", 2))
@@ -333,11 +340,11 @@ def test_higher_order_coordinate_map(points, celltype, order):
         i += 1
     x = cmap.push_forward(X, x_coord_new)
 
-    assert np.allclose(x[:, 0], X[:, 0])
-    assert np.allclose(x[:, 1], 2 * X[:, 1])
+    assert np.allclose(x[:, 0], X[:, 0], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
+    assert np.allclose(x[:, 1], 2 * X[:, 1], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
 
     if mesh.geometry.dim == 3:
-        assert np.allclose(x[:, 2], 3 * X[:, 2])
+        assert np.allclose(x[:, 2], 3 * X[:, 2], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
 
 
 @pytest.mark.skip_in_parallel
@@ -390,7 +397,9 @@ def test_higher_order_tetra_coordinate_map(order):
             ]
         )
     cells = np.array([range(len(points))])
-    domain = ufl.Mesh(element("Lagrange", celltype.name, order, shape=(3,)))
+    domain = ufl.Mesh(
+        element("Lagrange", celltype.name, order, shape=(3,), dtype=default_real_type)
+    )
     mesh = create_mesh(MPI.COMM_WORLD, cells, points, domain)
     V = functionspace(mesh, ("Lagrange", order))
     X = V.element.interpolation_points()
@@ -402,9 +411,9 @@ def test_higher_order_tetra_coordinate_map(order):
         x_coord_new[node] = x_g[x_dofs[0, node], : mesh.geometry.dim]
 
     x = mesh.geometry.cmap.push_forward(X, x_coord_new)
-    assert np.allclose(x[:, 0], X[:, 0])
-    assert np.allclose(x[:, 1], 2 * X[:, 1])
-    assert np.allclose(x[:, 2], 3 * X[:, 2])
+    assert np.allclose(x[:, 0], X[:, 0], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
+    assert np.allclose(x[:, 1], 2 * X[:, 1], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
+    assert np.allclose(x[:, 2], 3 * X[:, 2], atol=100 * np.finfo(mesh.geometry.x.dtype).eps)
 
 
 @pytest.mark.skip_in_parallel
diff --git a/python/test/unit/fem/test_element_integrals.py b/python/test/unit/fem/test_element_integrals.py
index 418212e9d78..6ec1be80a15 100644
--- a/python/test/unit/fem/test_element_integrals.py
+++ b/python/test/unit/fem/test_element_integrals.py
@@ -229,7 +229,7 @@ def test_facet_integral(cell_type, dtype):
             a = form(v * ufl.ds(subdomain_data=marker, subdomain_id=j), dtype=dtype)
             result = assemble_scalar(a)
             out.append(result)
-            assert np.isclose(result, out[0])
+            assert np.isclose(result, out[0], atol=np.finfo(dtype).eps)
 
 
 @pytest.mark.skip_in_parallel
@@ -300,7 +300,7 @@ def test_facet_normals(cell_type, dtype):
                     dtype=dtype,
                 )
                 result = assemble_scalar(a)
-                if np.isclose(result, 1):
+                if np.isclose(result, 1, atol=np.finfo(dtype).eps):
                     ones += 1
                 else:
                     assert np.isclose(result, 0, atol=1.0e-6)
@@ -326,7 +326,7 @@ def test_plus_minus(cell_type, space_type, dtype):
                 a = form(v(pm1) * v(pm2) * ufl.dS, dtype=dtype)
                 results.append(assemble_scalar(a))
     for i, j in combinations(results, 2):
-        assert np.isclose(i, j)
+        assert np.isclose(i, j, atol=np.finfo(dtype).eps)
 
 
 @pytest.mark.skip_in_parallel
@@ -378,7 +378,7 @@ def test_plus_minus_simple_vector(cell_type, pm, dtype):
                     # If no matching point found, fail
                     assert False
 
-                assert np.isclose(results[0][dof0], result[dof1])
+                assert np.isclose(results[0][dof0], result[dof1], atol=np.finfo(dtype).eps)
 
 
 @pytest.mark.skip_in_parallel
@@ -490,7 +490,7 @@ def test_plus_minus_matrix(cell_type, pm1, pm2, dtype):
         # For all dof pairs, check that entries in the matrix agree
         for a, b in dof_order:
             for c, d in dof_order:
-                assert np.isclose(results[0][a, c], result[b, d])
+                assert np.isclose(results[0][a, c], result[b, d], atol=np.finfo(dtype).eps)
 
 
 @pytest.mark.skip(
diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py
index e4eb977f0b5..dc9c5266ae1 100644
--- a/python/test/unit/fem/test_expression.py
+++ b/python/test/unit/fem/test_expression.py
@@ -61,7 +61,9 @@ def scatter(vec, array_evaluated, dofmap):
     b2 = Function(vdP1, dtype=dtype)
     b2.interpolate(lambda x: np.vstack((2.0 * x[0], 4.0 * x[1])))
 
-    assert np.allclose(b2.x.array, b.x.array, rtol=1.0e-5, atol=1.0e-5)
+    assert np.allclose(
+        b2.x.array, b.x.array, rtol=np.sqrt(np.finfo(dtype).eps), atol=np.sqrt(np.finfo(dtype).eps)
+    )
 
 
 @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128])
@@ -195,7 +197,12 @@ def exact_grad_f(x):
 
     # Evaluate exact gradient using global points
     grad_f_exact = exact_grad_f(x_evaluated)
-    assert np.allclose(grad_f_evaluated, grad_f_exact, rtol=1.0e-5, atol=1.0e-5)
+    assert np.allclose(
+        grad_f_evaluated,
+        grad_f_exact,
+        rtol=np.sqrt(np.finfo(dtype).eps),
+        atol=np.sqrt(np.finfo(dtype).eps),
+    )
 
 
 @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128])
diff --git a/python/test/unit/fem/test_fem_pipeline.py b/python/test/unit/fem/test_fem_pipeline.py
index c99c4070789..344b5b80e45 100644
--- a/python/test/unit/fem/test_fem_pipeline.py
+++ b/python/test/unit/fem/test_fem_pipeline.py
@@ -258,7 +258,7 @@ def test_curl_curl_eigenvalue(family, order):
         CellType.triangle,
     )
 
-    e = element(family, basix.CellType.triangle, order)
+    e = element(family, basix.CellType.triangle, order, dtype=default_real_type)
     V = functionspace(mesh, e)
 
     u = ufl.TrialFunction(V)
@@ -329,8 +329,8 @@ def test_biharmonic(family):
 
     e = mixed_element(
         [
-            element(family, basix.CellType.triangle, 1),
-            element(basix.ElementFamily.P, basix.CellType.triangle, 2),
+            element(family, basix.CellType.triangle, 1, dtype=default_real_type),
+            element(basix.ElementFamily.P, basix.CellType.triangle, 2, dtype=default_real_type),
         ]
     )
 
diff --git a/python/test/unit/fem/test_function.py b/python/test/unit/fem/test_function.py
index 53c51fbd115..550b0528e7a 100644
--- a/python/test/unit/fem/test_function.py
+++ b/python/test/unit/fem/test_function.py
@@ -128,7 +128,7 @@ def test_interpolation_mismatch_rank1(W):
 
 def test_mixed_element_interpolation():
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)
-    el = element("Lagrange", mesh.basix_cell(), 1)
+    el = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
     V = functionspace(mesh, mixed_element([el, el]))
     u = Function(V)
     with pytest.raises(RuntimeError):
diff --git a/python/test/unit/fem/test_function_space.py b/python/test/unit/fem/test_function_space.py
index b8d7dabab1e..fc9c89438b5 100644
--- a/python/test/unit/fem/test_function_space.py
+++ b/python/test/unit/fem/test_function_space.py
@@ -36,8 +36,10 @@ def W(mesh):
 
 @pytest.fixture
 def Q(mesh):
-    W = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
-    V = element("Lagrange", mesh.basix_cell(), 1)
+    W = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    V = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
     return functionspace(mesh, mixed_element([W, V]))
 
 
@@ -240,7 +242,7 @@ def test_argument_equality(mesh, V, V2, W, W2):
 
 def test_cell_mismatch(mesh):
     """Test that cell mismatch raises early enough from UFL"""
-    e = element("P", "triangle", 1)
+    e = element("P", "triangle", 1, dtype=default_real_type)
     with pytest.raises(BaseException):
         functionspace(mesh, e)
 
@@ -267,7 +269,7 @@ def test_vector_function_space_cell_type():
 
     # Create a mesh containing a single interval living in 2D
     cell = Cell("interval")
-    domain = Mesh(element("Lagrange", "interval", 1, shape=(1,)))
+    domain = Mesh(element("Lagrange", "interval", 1, shape=(1,), dtype=default_real_type))
     cells = np.array([[0, 1]], dtype=np.int64)
     x = np.array([[0.0, 0.0], [1.0, 1.0]])
     mesh = create_mesh(comm, cells, x, domain)
diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py
index 5a9d93205c2..a7c4e63c939 100644
--- a/python/test/unit/fem/test_interpolation.py
+++ b/python/test/unit/fem/test_interpolation.py
@@ -345,8 +345,10 @@ def test_mixed_sub_interpolation():
     def f(x):
         return np.vstack((10 + x[0], -10 - x[1], 25 + x[0]))
 
-    P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
-    P1 = element("Lagrange", mesh.basix_cell(), 1)
+    P2 = element(
+        "Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    P1 = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
     for i, P in enumerate((mixed_element([P2, P1]), mixed_element([P1, P2]))):
         W = functionspace(mesh, P)
         U = Function(W)
@@ -397,8 +399,10 @@ def f(x):
 def test_mixed_interpolation():
     """Test that mixed interpolation raised an exception."""
     mesh = one_cell_mesh(CellType.triangle)
-    A = element("Lagrange", mesh.basix_cell(), 1)
-    B = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
+    A = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
+    B = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
     v = Function(functionspace(mesh, mixed_element([A, B])))
     with pytest.raises(RuntimeError):
         v.interpolate(lambda x: (x[1], 2 * x[0], 3 * x[1]))
@@ -679,7 +683,7 @@ def test_de_rahm_2D(order):
     x = ufl.SpatialCoordinate(mesh)
     g_ex = ufl.as_vector((1 + x[1], 4 * x[1] + x[0]))
     assert np.abs(assemble_scalar(form(ufl.inner(q - g_ex, q - g_ex) * ufl.dx))) == pytest.approx(
-        0, abs=1e-10
+        0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
     )
 
     V = functionspace(mesh, ("BDM", order - 1))
@@ -691,7 +695,7 @@ def curl2D(u):
     v.interpolate(Expression(curl2D(ufl.grad(w)), V.element.interpolation_points()))
     h_ex = ufl.as_vector((1, -1))
     assert np.abs(assemble_scalar(form(ufl.inner(v - h_ex, v - h_ex) * ufl.dx))) == pytest.approx(
-        0, abs=1.0e-6
+        0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
     )
 
 
@@ -765,17 +769,27 @@ def test_interpolate_callable_subset(bound):
 @pytest.mark.parametrize(
     "scalar_element",
     [
-        element("P", "triangle", 1),
-        element("P", "triangle", 2),
-        element("P", "triangle", 3),
-        element("Q", "quadrilateral", 1),
-        element("Q", "quadrilateral", 2),
-        element("Q", "quadrilateral", 3),
-        element("S", "quadrilateral", 1),
-        element("S", "quadrilateral", 2),
-        element("S", "quadrilateral", 3),
-        enriched_element([element("P", "triangle", 1), element("Bubble", "triangle", 3)]),
-        enriched_element([element("P", "quadrilateral", 1), element("Bubble", "quadrilateral", 2)]),
+        element("P", "triangle", 1, dtype=default_real_type),
+        element("P", "triangle", 2, dtype=default_real_type),
+        element("P", "triangle", 3, dtype=default_real_type),
+        element("Q", "quadrilateral", 1, dtype=default_real_type),
+        element("Q", "quadrilateral", 2, dtype=default_real_type),
+        element("Q", "quadrilateral", 3, dtype=default_real_type),
+        element("S", "quadrilateral", 1, dtype=default_real_type),
+        element("S", "quadrilateral", 2, dtype=default_real_type),
+        element("S", "quadrilateral", 3, dtype=default_real_type),
+        enriched_element(
+            [
+                element("P", "triangle", 1, dtype=default_real_type),
+                element("Bubble", "triangle", 3, dtype=default_real_type),
+            ]
+        ),
+        enriched_element(
+            [
+                element("P", "quadrilateral", 1, dtype=default_real_type),
+                element("Bubble", "quadrilateral", 2, dtype=default_real_type),
+            ]
+        ),
     ],
 )
 def test_vector_element_interpolation(scalar_element):
@@ -821,6 +835,7 @@ def test_custom_vector_element():
         False,
         1,
         1,
+        dtype=default_real_type,
     )
 
     V = functionspace(mesh, e)
@@ -846,9 +861,11 @@ def g(x):
     x = ufl.SpatialCoordinate(mesh)
     dgdy = ufl.cos(x[1])
 
-    curl_el = element("N1curl", mesh.basix_cell(), 1)
-    vlag_el = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
-    lagr_el = element("Lagrange", mesh.basix_cell(), order)
+    curl_el = element("N1curl", mesh.basix_cell(), 1, dtype=default_real_type)
+    vlag_el = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    lagr_el = element("Lagrange", mesh.basix_cell(), order, dtype=default_real_type)
 
     V = functionspace(mesh, mixed_element([curl_el, lagr_el]))
     Eb_m = Function(V)
@@ -873,9 +890,9 @@ def test_nonmatching_mesh_interpolation(xtype, cell_type0, cell_type1):
     def f(x):
         return (7 * x[1], 3 * x[0], x[2] + 0.4)
 
-    el0 = element("Lagrange", mesh0.basix_cell(), 1, shape=(3,))
+    el0 = element("Lagrange", mesh0.basix_cell(), 1, shape=(3,), dtype=xtype)
     V0 = functionspace(mesh0, el0)
-    el1 = element("Lagrange", mesh1.basix_cell(), 1, shape=(3,))
+    el1 = element("Lagrange", mesh1.basix_cell(), 1, shape=(3,), dtype=xtype)
     V1 = functionspace(mesh1, el1)
 
     # Interpolate on 3D mesh
@@ -911,7 +928,12 @@ def f(x):
     u1_ex.interpolate(f)
     u1_ex.x.scatter_forward()
 
-    assert np.allclose(u1_ex.x.array, u1.x.array, rtol=1.0e-4, atol=1.0e-6)
+    assert np.allclose(
+        u1_ex.x.array,
+        u1.x.array,
+        rtol=np.sqrt(np.finfo(xtype).eps),
+        atol=np.sqrt(np.finfo(xtype).eps),
+    )
 
     # Interpolate 2D->3D
     u0_2 = Function(V0, dtype=xtype)
diff --git a/python/test/unit/fem/test_mixed_element.py b/python/test/unit/fem/test_mixed_element.py
index fe902dcab01..9663eee525c 100644
--- a/python/test/unit/fem/test_mixed_element.py
+++ b/python/test/unit/fem/test_mixed_element.py
@@ -12,6 +12,7 @@
 import dolfinx
 import ufl
 from basix.ufl import element, mixed_element
+from dolfinx import default_real_type
 from dolfinx.fem import form, functionspace
 from dolfinx.mesh import CellType, GhostMode, create_unit_cube, create_unit_square
 
@@ -32,7 +33,7 @@ def test_mixed_element(rank, family, cell, degree):
 
     shape = (mesh.geometry.dim,) * rank
     norms = []
-    U_el = element(family, cell.cellname(), degree, shape=shape)
+    U_el = element(family, cell.cellname(), degree, shape=shape, dtype=default_real_type)
     for i in range(3):
         U = functionspace(mesh, U_el)
         u = ufl.TrialFunction(U)
@@ -78,8 +79,10 @@ def test_vector_element():
 @pytest.mark.parametrize("d2", range(1, 4))
 def test_element_product(d1, d2):
     mesh = create_unit_square(MPI.COMM_WORLD, 2, 2)
-    P3 = element("Lagrange", mesh.basix_cell(), d1, shape=(mesh.geometry.dim,))
-    P1 = element("Lagrange", mesh.basix_cell(), d2)
+    P3 = element(
+        "Lagrange", mesh.basix_cell(), d1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    P1 = element("Lagrange", mesh.basix_cell(), d2, dtype=default_real_type)
     TH = mixed_element([P3, P1])
     W = functionspace(mesh, TH)
 
diff --git a/python/test/unit/fem/test_mixed_mesh_dofmap.py b/python/test/unit/fem/test_mixed_mesh_dofmap.py
index 7c528558cc4..b3745feb398 100644
--- a/python/test/unit/fem/test_mixed_mesh_dofmap.py
+++ b/python/test/unit/fem/test_mixed_mesh_dofmap.py
@@ -4,7 +4,6 @@
 
 import basix
 import dolfinx.cpp as _cpp
-from dolfinx import jit
 from dolfinx.cpp.mesh import Mesh_float64, create_geometry, create_topology
 from dolfinx.fem import coordinate_element
 from dolfinx.fem.dofmap import DofMap
@@ -14,22 +13,11 @@
 
 def create_element_dofmap(mesh, cell_types, degree):
     cpp_elements = []
-    dofmaps = []
     for cell_type in cell_types:
-        ufl_e = basix.ufl.element("P", cell_type, degree)
-        form_compiler_options = {"scalar_type": np.float64}
-        (ufcx_element, ufcx_dofmap), module, code = jit.ffcx_jit(
-            mesh.comm, ufl_e, form_compiler_options=form_compiler_options
-        )
-        ffi = module.ffi
-        cpp_elements += [
-            _cpp.fem.FiniteElement_float64(ffi.cast("uintptr_t", ffi.addressof(ufcx_element)))
-        ]
-        dofmaps += [ufcx_dofmap]
-
-    cpp_dofmaps = _cpp.fem.create_dofmaps(
-        mesh.comm, [ffi.cast("uintptr_t", ffi.addressof(dm)) for dm in dofmaps], mesh.topology
-    )
+        ufl_e = basix.ufl.element("P", cell_type, degree, dtype=np.float64)
+        cpp_elements += [_cpp.fem.FiniteElement_float64(ufl_e.basix_element._e, 1, False)]
+
+    cpp_dofmaps = _cpp.fem.create_dofmaps(mesh.comm, mesh.topology, cpp_elements)
 
     return (cpp_elements, cpp_dofmaps)
 
diff --git a/python/test/unit/fem/test_nonlinear_assembler.py b/python/test/unit/fem/test_nonlinear_assembler.py
index d952a2d02fc..1fe195172c7 100644
--- a/python/test/unit/fem/test_nonlinear_assembler.py
+++ b/python/test/unit/fem/test_nonlinear_assembler.py
@@ -15,6 +15,7 @@
 
 import ufl
 from basix.ufl import element, mixed_element
+from dolfinx import default_real_type
 from dolfinx.cpp.la.petsc import scatter_local_vectors
 from dolfinx.fem import (
     Function,
@@ -67,8 +68,8 @@ def test_matrix_assembly_block_nl():
     in the nonlinear setting."""
     mesh = create_unit_square(MPI.COMM_WORLD, 4, 8)
     p0, p1 = 1, 2
-    P0 = element("Lagrange", mesh.basix_cell(), p0)
-    P1 = element("Lagrange", mesh.basix_cell(), p1)
+    P0 = element("Lagrange", mesh.basix_cell(), p0, dtype=default_real_type)
+    P1 = element("Lagrange", mesh.basix_cell(), p1, dtype=default_real_type)
     V0 = functionspace(mesh, P0)
     V1 = functionspace(mesh, P1)
 
@@ -314,7 +315,7 @@ def test_assembly_solve_block_nl():
     matrix approaches and test that solution is the same."""
     mesh = create_unit_square(MPI.COMM_WORLD, 12, 11)
     p = 1
-    P = element("Lagrange", mesh.basix_cell(), p)
+    P = element("Lagrange", mesh.basix_cell(), p, dtype=default_real_type)
     V0 = functionspace(mesh, P)
     V1 = V0.clone()
 
@@ -641,8 +642,10 @@ def nested():
 
     def monolithic():
         """Monolithic"""
-        P2_el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
-        P1_el = element("Lagrange", mesh.basix_cell(), 1)
+        P2_el = element(
+            "Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), dtype=default_real_type
+        )
+        P1_el = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
         TH = mixed_element([P2_el, P1_el])
         W = functionspace(mesh, TH)
         U = Function(W)
diff --git a/python/test/unit/fem/test_petsc_discrete_operators.py b/python/test/unit/fem/test_petsc_discrete_operators.py
index cb40a8c4b1a..83ea90efa17 100644
--- a/python/test/unit/fem/test_petsc_discrete_operators.py
+++ b/python/test/unit/fem/test_petsc_discrete_operators.py
@@ -126,8 +126,10 @@ def test_interpolation_matrix_petsc(cell_type, p, q, from_lagrange):
         mesh = create_unit_cube(comm, 3, 2, 2, ghost_mode=GhostMode.none, cell_type=cell_type)
         lagrange = "Lagrange" if from_lagrange else "DG"
         nedelec = "Nedelec 1st kind H(curl)"
-    v_el = element(lagrange, mesh.basix_cell(), p, shape=(mesh.geometry.dim,))
-    s_el = element(nedelec, mesh.basix_cell(), q)
+    v_el = element(
+        lagrange, mesh.basix_cell(), p, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    s_el = element(nedelec, mesh.basix_cell(), q, dtype=default_real_type)
     if from_lagrange:
         el0 = v_el
         el1 = s_el
diff --git a/python/test/unit/fem/test_symmetric_tensor.py b/python/test/unit/fem/test_symmetric_tensor.py
index cf374c19bbb..588890a3d58 100644
--- a/python/test/unit/fem/test_symmetric_tensor.py
+++ b/python/test/unit/fem/test_symmetric_tensor.py
@@ -12,7 +12,14 @@
 @pytest.mark.parametrize("symmetry", [True, False])
 def test_transpose(degree, symmetry):
     mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
-    e = basix.ufl.element("Lagrange", "triangle", degree, shape=(2, 2), symmetry=symmetry)
+    e = basix.ufl.element(
+        "Lagrange",
+        "triangle",
+        degree,
+        shape=(2, 2),
+        symmetry=symmetry,
+        dtype=dolfinx.default_real_type,
+    )
 
     space = dolfinx.fem.functionspace(mesh, e)
 
@@ -31,8 +38,12 @@ def tensor(x):
         mat = np.array([[0], [1], [2], [1], [3], [4], [2], [4], [5]])
         return np.broadcast_to(mat, (9, x.shape[1]))
 
-    element = basix.ufl.element("DG", mesh.basix_cell(), 0, shape=(3, 3))
-    symm_element = basix.ufl.element("DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=True)
+    element = basix.ufl.element(
+        "DG", mesh.basix_cell(), 0, shape=(3, 3), dtype=dolfinx.default_real_type
+    )
+    symm_element = basix.ufl.element(
+        "DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=True, dtype=dolfinx.default_real_type
+    )
     space = dolfinx.fem.functionspace(mesh, element)
     symm_space = dolfinx.fem.functionspace(mesh, symm_element)
     f = dolfinx.fem.Function(space)
@@ -55,7 +66,9 @@ def test_eval():
     def tensor(x):
         return np.broadcast_to(mat.reshape((9, 1)), (9, x.shape[1]))
 
-    element = basix.ufl.element("DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=True)
+    element = basix.ufl.element(
+        "DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=True, dtype=dolfinx.default_real_type
+    )
     space = dolfinx.fem.functionspace(mesh, element)
     f = dolfinx.fem.Function(space)
 
diff --git a/python/test/unit/fem/test_symmetry.py b/python/test/unit/fem/test_symmetry.py
index 9a167e3ca88..48e5a29650a 100644
--- a/python/test/unit/fem/test_symmetry.py
+++ b/python/test/unit/fem/test_symmetry.py
@@ -135,8 +135,8 @@ def test_mixed_element_form(cell_type, sign, order, dtype):
 
     U_el = mixed_element(
         [
-            element(basix.ElementFamily.P, cell_type.name, order),
-            element(basix.ElementFamily.N1E, cell_type.name, order),
+            element(basix.ElementFamily.P, cell_type.name, order, dtype=dtype),
+            element(basix.ElementFamily.N1E, cell_type.name, order, dtype=dtype),
         ]
     )
 
@@ -164,8 +164,14 @@ def test_mixed_element_vector_element_form(cell_type, sign, order, dtype):
 
     U_el = mixed_element(
         [
-            element(basix.ElementFamily.P, cell_type.name, order, shape=(mesh.geometry.dim,)),
-            element(basix.ElementFamily.N1E, cell_type.name, order),
+            element(
+                basix.ElementFamily.P,
+                cell_type.name,
+                order,
+                shape=(mesh.geometry.dim,),
+                dtype=dtype,
+            ),
+            element(basix.ElementFamily.N1E, cell_type.name, order, dtype=dtype),
         ]
     )
 
diff --git a/python/test/unit/io/test_adios2.py b/python/test/unit/io/test_adios2.py
index 907da78481d..beb26db1af5 100644
--- a/python/test/unit/io/test_adios2.py
+++ b/python/test/unit/io/test_adios2.py
@@ -243,7 +243,7 @@ def vel(x):
 def test_save_vtkx_cell_point(tempdir):
     """Test writing point-wise data"""
     mesh = create_unit_square(MPI.COMM_WORLD, 8, 5)
-    P = element("Discontinuous Lagrange", mesh.basix_cell(), 0)
+    P = element("Discontinuous Lagrange", mesh.basix_cell(), 0, dtype=default_real_type)
 
     V = functionspace(mesh, P)
     u = Function(V)
diff --git a/python/test/unit/io/test_vtk.py b/python/test/unit/io/test_vtk.py
index 52611320a41..2d0c486e9e7 100644
--- a/python/test/unit/io/test_vtk.py
+++ b/python/test/unit/io/test_vtk.py
@@ -99,7 +99,7 @@ def f(x):
         vals[1] = 2 * x[0] * x[0]
         return vals
 
-    e = element("Lagrange", mesh.basix_cell(), 2, shape=(2,))
+    e = element("Lagrange", mesh.basix_cell(), 2, shape=(2,), dtype=default_real_type)
     u = Function(functionspace(mesh, e))
     u.interpolate(f)
     filename = Path(tempdir, "u.pvd")
@@ -149,8 +149,10 @@ def test_save_2d_vector_CG2(tempdir):
 
 def test_save_vtk_mixed(tempdir):
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)
-    P2 = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
-    P1 = element("Lagrange", mesh.basix_cell(), 1)
+    P2 = element(
+        "Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,), dtype=default_real_type
+    )
+    P1 = element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type)
     W = functionspace(mesh, mixed_element([P2, P1]))
     V1 = functionspace(mesh, P1)
     V2 = functionspace(mesh, P2)
@@ -203,8 +205,8 @@ def f(x):
 def test_save_vtk_cell_point(tempdir):
     """Test writing cell-wise and point-wise data"""
     mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)
-    P2 = element("Lagrange", mesh.basix_cell(), 1, shape=(3,))
-    P1 = element("Discontinuous Lagrange", mesh.basix_cell(), 0)
+    P2 = element("Lagrange", mesh.basix_cell(), 1, shape=(3,), dtype=default_real_type)
+    P1 = element("Discontinuous Lagrange", mesh.basix_cell(), 0, dtype=default_real_type)
 
     V2, V1 = functionspace(mesh, P2), functionspace(mesh, P1)
     U2, U1 = Function(V2), Function(V1)
@@ -222,7 +224,7 @@ def test_save_vtk_cell_point(tempdir):
 
 def test_save_1d_tensor(tempdir):
     mesh = create_unit_interval(MPI.COMM_WORLD, 32)
-    e = element("Lagrange", mesh.basix_cell(), 2, shape=(2, 2))
+    e = element("Lagrange", mesh.basix_cell(), 2, shape=(2, 2), dtype=default_real_type)
     u = Function(functionspace(mesh, e))
     u.x.array[:] = 1.0
     filename = Path(tempdir, "u.pvd")
diff --git a/python/test/unit/io/test_xdmf_function.py b/python/test/unit/io/test_xdmf_function.py
index b231bb17a9a..dfcadf4c48f 100644
--- a/python/test/unit/io/test_xdmf_function.py
+++ b/python/test/unit/io/test_xdmf_function.py
@@ -12,6 +12,7 @@
 import pytest
 
 import basix
+from dolfinx import default_real_type
 from dolfinx.fem import Function, functionspace
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import CellType, create_unit_cube, create_unit_interval, create_unit_square
@@ -339,7 +340,11 @@ def gmsh_hex_model(order):
 
     # Write P3 GLL Function (exception expected)
     ufl_e = basix.ufl.element(
-        basix.ElementFamily.P, basix.CellType.tetrahedron, 3, basix.LagrangeVariant.gll_warped
+        basix.ElementFamily.P,
+        basix.CellType.tetrahedron,
+        3,
+        basix.LagrangeVariant.gll_warped,
+        dtype=default_real_type,
     )
     u = Function(functionspace(msh, ufl_e))
     with pytest.raises(RuntimeError):
@@ -350,7 +355,11 @@ def gmsh_hex_model(order):
 
     # Write P3 equispaced Function
     ufl_e = basix.ufl.element(
-        basix.ElementFamily.P, basix.CellType.tetrahedron, 3, basix.LagrangeVariant.equispaced
+        basix.ElementFamily.P,
+        basix.CellType.tetrahedron,
+        3,
+        basix.LagrangeVariant.equispaced,
+        dtype=default_real_type,
     )
     u1 = Function(functionspace(msh, ufl_e))
     u1.interpolate(u)
diff --git a/python/test/unit/la/test_vector_scatter.py b/python/test/unit/la/test_vector_scatter.py
index 729997a84de..b77c1ee7593 100644
--- a/python/test/unit/la/test_vector_scatter.py
+++ b/python/test/unit/la/test_vector_scatter.py
@@ -11,13 +11,17 @@
 import pytest
 
 from basix.ufl import element
-from dolfinx import la
+from dolfinx import default_real_type, la
 from dolfinx.fem import Function, functionspace
 from dolfinx.mesh import create_unit_square
 
 
 @pytest.mark.parametrize(
-    "e", [element("Lagrange", "triangle", 1), element("Lagrange", "triangle", 1, shape=(2,))]
+    "e",
+    [
+        element("Lagrange", "triangle", 1, dtype=default_real_type),
+        element("Lagrange", "triangle", 1, shape=(2,), dtype=default_real_type),
+    ],
 )
 def test_scatter_forward(e):
     mesh = create_unit_square(MPI.COMM_WORLD, 5, 5)
@@ -46,7 +50,11 @@ def test_scatter_forward(e):
 
 
 @pytest.mark.parametrize(
-    "e", [element("Lagrange", "triangle", 1), element("Lagrange", "triangle", 1, shape=(2,))]
+    "e",
+    [
+        element("Lagrange", "triangle", 1, dtype=default_real_type),
+        element("Lagrange", "triangle", 1, shape=(2,), dtype=default_real_type),
+    ],
 )
 def test_scatter_reverse(e):
     comm = MPI.COMM_WORLD