diff --git a/tests/test_create_fiat_element.py b/tests/test_create_fiat_element.py
index 0924edb3..b54d9782 100644
--- a/tests/test_create_fiat_element.py
+++ b/tests/test_create_fiat_element.py
@@ -4,6 +4,7 @@
 from FIAT.discontinuous_lagrange import HigherOrderDiscontinuousLagrange as FIAT_DiscontinuousLagrange
 
 import ufl
+import finat.ufl
 from tsfc.finatinterface import create_element as _create_element
 
 
@@ -40,7 +41,7 @@ def triangle_names(request):
 
 @pytest.fixture
 def ufl_element(triangle_names):
-    return ufl.FiniteElement(triangle_names, ufl.triangle, 2)
+    return finat.ufl.FiniteElement(triangle_names, ufl.triangle, 2)
 
 
 def test_triangle_basic(ufl_element):
@@ -58,16 +59,16 @@ def tensor_name(request):
                 ids=lambda x: x.cellname(),
                 scope="module")
 def ufl_A(request, tensor_name):
-    return ufl.FiniteElement(tensor_name, request.param, 1)
+    return finat.ufl.FiniteElement(tensor_name, request.param, 1)
 
 
 @pytest.fixture
 def ufl_B(tensor_name):
-    return ufl.FiniteElement(tensor_name, ufl.interval, 1)
+    return finat.ufl.FiniteElement(tensor_name, ufl.interval, 1)
 
 
 def test_tensor_prod_simple(ufl_A, ufl_B):
-    tensor_ufl = ufl.TensorProductElement(ufl_A, ufl_B)
+    tensor_ufl = finat.ufl.TensorProductElement(ufl_A, ufl_B)
 
     tensor = create_element(tensor_ufl)
     A = create_element(ufl_A)
@@ -84,7 +85,7 @@ def test_tensor_prod_simple(ufl_A, ufl_B):
                           ('DP', FIAT.GaussLegendre),
                           ('DP L2', FIAT.GaussLegendre)])
 def test_interval_variant_default(family, expected_cls):
-    ufl_element = ufl.FiniteElement(family, ufl.interval, 3)
+    ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3)
     assert isinstance(create_element(ufl_element), expected_cls)
 
 
@@ -96,42 +97,42 @@ def test_interval_variant_default(family, expected_cls):
                           ('DP L2', 'equispaced', FIAT_DiscontinuousLagrange),
                           ('DP L2', 'spectral', FIAT.GaussLegendre)])
 def test_interval_variant(family, variant, expected_cls):
-    ufl_element = ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
+    ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
     assert isinstance(create_element(ufl_element), expected_cls)
 
 
 def test_triangle_variant_spectral_fail():
-    ufl_element = ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral')
+    ufl_element = finat.ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral')
     with pytest.raises(ValueError):
         create_element(ufl_element)
 
 
 def test_triangle_variant_spectral_fail_l2():
-    ufl_element = ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
+    ufl_element = finat.ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
     with pytest.raises(ValueError):
         create_element(ufl_element)
 
 
 def test_quadrilateral_variant_spectral_q():
-    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
     assert isinstance(element.element.A, FIAT.GaussLobattoLegendre)
     assert isinstance(element.element.B, FIAT.GaussLobattoLegendre)
 
 
 def test_quadrilateral_variant_spectral_dq():
-    element = create_element(ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral'))
     assert isinstance(element.element.A, FIAT.GaussLegendre)
     assert isinstance(element.element.B, FIAT.GaussLegendre)
 
 
 def test_quadrilateral_variant_spectral_dq_l2():
-    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
     assert isinstance(element.element.A, FIAT.GaussLegendre)
     assert isinstance(element.element.B, FIAT.GaussLegendre)
 
 
 def test_quadrilateral_variant_spectral_rtcf():
-    element = create_element(ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='spectral'))
     assert isinstance(element.element._elements[0].A, FIAT.GaussLobattoLegendre)
     assert isinstance(element.element._elements[0].B, FIAT.GaussLegendre)
     assert isinstance(element.element._elements[1].A, FIAT.GaussLegendre)
diff --git a/tests/test_create_finat_element.py b/tests/test_create_finat_element.py
index 2f937a86..93001065 100644
--- a/tests/test_create_finat_element.py
+++ b/tests/test_create_finat_element.py
@@ -1,6 +1,7 @@
 import pytest
 
 import ufl
+import finat.ufl
 import finat
 from tsfc.finatinterface import create_element, supported_elements
 
@@ -18,7 +19,7 @@ def triangle_names(request):
 
 @pytest.fixture
 def ufl_element(triangle_names):
-    return ufl.FiniteElement(triangle_names, ufl.triangle, 2)
+    return finat.ufl.FiniteElement(triangle_names, ufl.triangle, 2)
 
 
 def test_triangle_basic(ufl_element):
@@ -28,7 +29,7 @@ def test_triangle_basic(ufl_element):
 
 @pytest.fixture
 def ufl_vector_element(triangle_names):
-    return ufl.VectorElement(triangle_names, ufl.triangle, 2)
+    return finat.ufl.VectorElement(triangle_names, ufl.triangle, 2)
 
 
 def test_triangle_vector(ufl_element, ufl_vector_element):
@@ -48,16 +49,16 @@ def tensor_name(request):
                         ufl.quadrilateral],
                 ids=lambda x: x.cellname())
 def ufl_A(request, tensor_name):
-    return ufl.FiniteElement(tensor_name, request.param, 1)
+    return finat.ufl.FiniteElement(tensor_name, request.param, 1)
 
 
 @pytest.fixture
 def ufl_B(tensor_name):
-    return ufl.FiniteElement(tensor_name, ufl.interval, 1)
+    return finat.ufl.FiniteElement(tensor_name, ufl.interval, 1)
 
 
 def test_tensor_prod_simple(ufl_A, ufl_B):
-    tensor_ufl = ufl.TensorProductElement(ufl_A, ufl_B)
+    tensor_ufl = finat.ufl.TensorProductElement(ufl_A, ufl_B)
 
     tensor = create_element(tensor_ufl)
     A = create_element(ufl_A)
@@ -73,7 +74,7 @@ def test_tensor_prod_simple(ufl_A, ufl_B):
                           ('DP', finat.GaussLegendre),
                           ('DP L2', finat.GaussLegendre)])
 def test_interval_variant_default(family, expected_cls):
-    ufl_element = ufl.FiniteElement(family, ufl.interval, 3)
+    ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3)
     assert isinstance(create_element(ufl_element), expected_cls)
 
 
@@ -85,36 +86,36 @@ def test_interval_variant_default(family, expected_cls):
                           ('DP L2', 'equispaced', finat.DiscontinuousLagrange),
                           ('DP L2', 'spectral', finat.GaussLegendre)])
 def test_interval_variant(family, variant, expected_cls):
-    ufl_element = ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
+    ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
     assert isinstance(create_element(ufl_element), expected_cls)
 
 
 def test_triangle_variant_spectral_fail():
-    ufl_element = ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral')
+    ufl_element = finat.ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral')
     with pytest.raises(ValueError):
         create_element(ufl_element)
 
 
 def test_triangle_variant_spectral_fail_l2():
-    ufl_element = ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
+    ufl_element = finat.ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
     with pytest.raises(ValueError):
         create_element(ufl_element)
 
 
 def test_quadrilateral_variant_spectral_q():
-    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
     assert isinstance(element.product.factors[0], finat.GaussLobattoLegendre)
     assert isinstance(element.product.factors[1], finat.GaussLobattoLegendre)
 
 
 def test_quadrilateral_variant_spectral_dq():
-    element = create_element(ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral'))
     assert isinstance(element.product.factors[0], finat.GaussLegendre)
     assert isinstance(element.product.factors[1], finat.GaussLegendre)
 
 
 def test_quadrilateral_variant_spectral_dq_l2():
-    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
+    element = create_element(finat.ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
     assert isinstance(element.product.factors[0], finat.GaussLegendre)
     assert isinstance(element.product.factors[1], finat.GaussLegendre)
 
diff --git a/tests/test_dual_evaluation.py b/tests/test_dual_evaluation.py
index 461226c2..b4f6e977 100644
--- a/tests/test_dual_evaluation.py
+++ b/tests/test_dual_evaluation.py
@@ -1,12 +1,13 @@
 import pytest
 import ufl
+import finat.ufl
 from tsfc.finatinterface import create_element
 from tsfc import compile_expression_dual_evaluation
 
 
 def test_ufl_only_simple():
-    mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
+    mesh = ufl.Mesh(finat.ufl.VectorElement("P", ufl.triangle, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("P", ufl.triangle, 2))
     v = ufl.Coefficient(V)
     expr = ufl.inner(v, v)
     W = V
@@ -16,8 +17,8 @@ def test_ufl_only_simple():
 
 
 def test_ufl_only_spatialcoordinate():
-    mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
+    mesh = ufl.Mesh(finat.ufl.VectorElement("P", ufl.triangle, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("P", ufl.triangle, 2))
     x, y = ufl.SpatialCoordinate(mesh)
     expr = x*y - y**2 + x
     W = V
@@ -27,30 +28,30 @@ def test_ufl_only_spatialcoordinate():
 
 
 def test_ufl_only_from_contravariant_piola():
-    mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
+    mesh = ufl.Mesh(finat.ufl.VectorElement("P", ufl.triangle, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("RT", ufl.triangle, 1))
     v = ufl.Coefficient(V)
     expr = ufl.inner(v, v)
-    W = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
+    W = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("P", ufl.triangle, 2))
     to_element = create_element(W.ufl_element())
     kernel = compile_expression_dual_evaluation(expr, to_element, W.ufl_element())
     assert kernel.needs_external_coords is True
 
 
 def test_ufl_only_to_contravariant_piola():
-    mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
+    mesh = ufl.Mesh(finat.ufl.VectorElement("P", ufl.triangle, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("P", ufl.triangle, 2))
     v = ufl.Coefficient(V)
     expr = ufl.as_vector([v, v])
-    W = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
+    W = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("RT", ufl.triangle, 1))
     to_element = create_element(W.ufl_element())
     kernel = compile_expression_dual_evaluation(expr, to_element, W.ufl_element())
     assert kernel.needs_external_coords is True
 
 
 def test_ufl_only_shape_mismatch():
-    mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
+    mesh = ufl.Mesh(finat.ufl.VectorElement("P", ufl.triangle, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement("RT", ufl.triangle, 1))
     v = ufl.Coefficient(V)
     expr = ufl.inner(v, v)
     assert expr.ufl_shape == ()
diff --git a/tests/test_estimated_degree.py b/tests/test_estimated_degree.py
index c3c80f02..e4842f2e 100644
--- a/tests/test_estimated_degree.py
+++ b/tests/test_estimated_degree.py
@@ -3,6 +3,7 @@
 import pytest
 
 import ufl
+import finat.ufl
 from tsfc import compile_form
 from tsfc.logging import logger
 
@@ -14,8 +15,8 @@ def emit(self, record):
 
 def test_estimated_degree():
     cell = ufl.tetrahedron
-    mesh = ufl.Mesh(ufl.VectorElement('P', cell, 1))
-    V = ufl.FunctionSpace(mesh, ufl.FiniteElement('P', cell, 1))
+    mesh = ufl.Mesh(finat.ufl.VectorElement('P', cell, 1))
+    V = ufl.FunctionSpace(mesh, finat.ufl.FiniteElement('P', cell, 1))
     f = ufl.Coefficient(V)
     u = ufl.TrialFunction(V)
     v = ufl.TestFunction(V)
diff --git a/tests/test_firedrake_972.py b/tests/test_firedrake_972.py
index 0c0cab0e..f9ec8d09 100644
--- a/tests/test_firedrake_972.py
+++ b/tests/test_firedrake_972.py
@@ -1,8 +1,9 @@
 import numpy
 import pytest
 
-from ufl import (Mesh, FunctionSpace, VectorElement, TensorElement,
+from ufl import (Mesh, FunctionSpace,
                  Coefficient, TestFunction, interval, indices, dx)
+from finat.ufl import VectorElement, TensorElement
 from ufl.classes import IndexSum, Product, MultiIndex
 
 from tsfc import compile_form
diff --git a/tests/test_gem_failure.py b/tests/test_gem_failure.py
index 8bf31338..0a79a39a 100644
--- a/tests/test_gem_failure.py
+++ b/tests/test_gem_failure.py
@@ -1,5 +1,6 @@
-from ufl import (triangle, tetrahedron, FiniteElement,
+from ufl import (triangle, tetrahedron, FunctionSpace, Mesh,
                  TrialFunction, TestFunction, inner, grad, dx, dS)
+from finat.ufl import FiniteElement, VectorElement
 from tsfc import compile_form
 from FIAT.hdiv_trace import TraceError
 import pytest
@@ -12,8 +13,10 @@ def test_cell_error(cell, degree):
     cell triggers `gem.Failure` to raise the TraceError exception.
     """
     trace_element = FiniteElement("HDiv Trace", cell, degree)
-    lambdar = TrialFunction(trace_element)
-    gammar = TestFunction(trace_element)
+    domain = Mesh(VectorElement("Lagrange", cell, 1))
+    space = FunctionSpace(domain, trace_element)
+    lambdar = TrialFunction(space)
+    gammar = TestFunction(space)
 
     with pytest.raises(TraceError):
         compile_form(lambdar * gammar * dx)
@@ -27,8 +30,10 @@ def test_gradient_error(cell, degree):
     exception.
     """
     trace_element = FiniteElement("HDiv Trace", cell, degree)
-    lambdar = TrialFunction(trace_element)
-    gammar = TestFunction(trace_element)
+    domain = Mesh(VectorElement("Lagrange", cell, 1))
+    space = FunctionSpace(domain, trace_element)
+    lambdar = TrialFunction(space)
+    gammar = TestFunction(space)
 
     with pytest.raises(TraceError):
         compile_form(inner(grad(lambdar('+')), grad(gammar('+'))) * dS)
diff --git a/tests/test_idempotency.py b/tests/test_idempotency.py
index 9ecb0f37..cf554f46 100644
--- a/tests/test_idempotency.py
+++ b/tests/test_idempotency.py
@@ -1,4 +1,5 @@
 import ufl
+import finat.ufl
 from tsfc import compile_form
 import loopy
 import pytest
@@ -21,13 +22,13 @@ def coord_degree(request):
 
 @pytest.fixture
 def mesh(cell, coord_degree):
-    c = ufl.VectorElement("CG", cell, coord_degree)
+    c = finat.ufl.VectorElement("CG", cell, coord_degree)
     return ufl.Mesh(c)
 
 
-@pytest.fixture(params=[ufl.FiniteElement,
-                        ufl.VectorElement,
-                        ufl.TensorElement],
+@pytest.fixture(params=[finat.ufl.FiniteElement,
+                        finat.ufl.VectorElement,
+                        finat.ufl.TensorElement],
                 ids=["FE", "VE", "TE"])
 def V(request, mesh):
     return ufl.FunctionSpace(mesh, request.param("CG", mesh.ufl_cell(), 2))
diff --git a/tests/test_impero_loopy_flop_counts.py b/tests/test_impero_loopy_flop_counts.py
index da357cf2..76bdad96 100644
--- a/tests/test_impero_loopy_flop_counts.py
+++ b/tests/test_impero_loopy_flop_counts.py
@@ -5,10 +5,11 @@
 import numpy
 import loopy
 from tsfc import compile_form
-from ufl import (FiniteElement, FunctionSpace, Mesh, TestFunction,
-                 TrialFunction, VectorElement, dx, grad, inner,
+from ufl import (FunctionSpace, Mesh, TestFunction,
+                 TrialFunction, dx, grad, inner,
                  interval, triangle, quadrilateral,
                  TensorProductCell)
+from finat.ufl import FiniteElement, VectorElement
 from tsfc.parameters import target
 
 
diff --git a/tests/test_interpolation_factorisation.py b/tests/test_interpolation_factorisation.py
index b588b34d..f2404a84 100644
--- a/tests/test_interpolation_factorisation.py
+++ b/tests/test_interpolation_factorisation.py
@@ -2,9 +2,9 @@
 import numpy
 import pytest
 
-from ufl import (Mesh, FunctionSpace, FiniteElement, VectorElement,
-                 TensorElement, Coefficient,
+from ufl import (Mesh, FunctionSpace, Coefficient,
                  interval, quadrilateral, hexahedron)
+from finat.ufl import FiniteElement, VectorElement, TensorElement
 
 from tsfc import compile_expression_dual_evaluation
 from tsfc.finatinterface import create_element
@@ -54,11 +54,11 @@ def test_sum_factorisation_scalar_tensor(mesh, element):
     source = element(degree - 1)
     target = element(degree)
     tensor_flops = flop_count(mesh, source, target)
-    expect = numpy.prod(target.value_shape())
+    expect = numpy.prod(target.value_shape)
     if isinstance(target, FiniteElement):
         scalar_flops = tensor_flops
     else:
-        target = target.sub_elements()[0]
-        source = source.sub_elements()[0]
+        target = target.sub_elements[0]
+        source = source.sub_elements[0]
         scalar_flops = flop_count(mesh, source, target)
     assert numpy.allclose(tensor_flops / scalar_flops, expect, rtol=1e-2)
diff --git a/tests/test_sum_factorisation.py b/tests/test_sum_factorisation.py
index cb6a992e..3e785c5b 100644
--- a/tests/test_sum_factorisation.py
+++ b/tests/test_sum_factorisation.py
@@ -1,11 +1,11 @@
 import numpy
 import pytest
 
-from ufl import (Mesh, FunctionSpace, FiniteElement, VectorElement,
-                 TestFunction, TrialFunction, TensorProductCell,
-                 EnrichedElement, HCurlElement, HDivElement,
-                 TensorProductElement, dx, action, interval, triangle,
+from ufl import (Mesh, FunctionSpace, TestFunction, TrialFunction,
+                 TensorProductCell, dx, action, interval, triangle,
                  quadrilateral, curl, dot, div, grad)
+from finat.ufl import (FiniteElement, VectorElement, EnrichedElement,
+                       TensorProductElement, HCurlElement, HDivElement)
 
 from tsfc import compile_form
 
diff --git a/tests/test_tensor.py b/tests/test_tensor.py
index 39485abb..9d09a467 100644
--- a/tests/test_tensor.py
+++ b/tests/test_tensor.py
@@ -1,9 +1,10 @@
 import numpy
 import pytest
 
-from ufl import (Mesh, FunctionSpace, FiniteElement, VectorElement,
+from ufl import (Mesh, FunctionSpace,
                  Coefficient, TestFunction, TrialFunction, dx, div,
                  inner, interval, triangle, tetrahedron, dot, grad)
+from finat.ufl import FiniteElement, VectorElement
 
 from tsfc import compile_form
 
diff --git a/tests/test_tsfc_182.py b/tests/test_tsfc_182.py
index 47735630..556a6baf 100644
--- a/tests/test_tsfc_182.py
+++ b/tests/test_tsfc_182.py
@@ -1,7 +1,7 @@
 import pytest
 
-from ufl import (Coefficient, FiniteElement, MixedElement,
-                 TestFunction, VectorElement, dx, inner, tetrahedron)
+from ufl import Coefficient, TestFunction, dx, inner, tetrahedron, Mesh, FunctionSpace
+from finat.ufl import FiniteElement, MixedElement, VectorElement
 
 from tsfc import compile_form
 
@@ -19,9 +19,11 @@ def test_delta_elimination(mode):
                                   dim=6, quad_scheme=scheme)
 
     element_chi_lambda = MixedElement(element_eps_p, element_lambda)
+    domain = Mesh(VectorElement("Lagrange", tetrahedron, 1))
+    space = FunctionSpace(domain, element_chi_lambda)
 
-    chi_lambda = Coefficient(element_chi_lambda)
-    delta_chi_lambda = TestFunction(element_chi_lambda)
+    chi_lambda = Coefficient(space)
+    delta_chi_lambda = TestFunction(space)
 
     L = inner(delta_chi_lambda, chi_lambda) * dx(degree=degree, scheme=scheme)
     kernel, = compile_form(L, parameters={'mode': mode})
diff --git a/tests/test_tsfc_204.py b/tests/test_tsfc_204.py
index e41e9c0c..f6108c7c 100644
--- a/tests/test_tsfc_204.py
+++ b/tests/test_tsfc_204.py
@@ -1,7 +1,8 @@
 from tsfc import compile_form
-from ufl import (BrokenElement, Coefficient, FacetNormal, FiniteElement,
-                 FunctionSpace, Mesh, MixedElement, VectorElement, as_matrix,
+from ufl import (Coefficient, FacetNormal,
+                 FunctionSpace, Mesh, as_matrix,
                  dot, dS, ds, dx, facet, grad, inner, outer, split, triangle)
+from finat.ufl import BrokenElement, FiniteElement, MixedElement, VectorElement
 
 
 def test_physically_mapped_facet():
diff --git a/tests/test_tsfc_274.py b/tests/test_tsfc_274.py
index 00f7a96a..453d8746 100644
--- a/tests/test_tsfc_274.py
+++ b/tests/test_tsfc_274.py
@@ -3,7 +3,8 @@
 from finat.point_set import PointSet
 from gem.interpreter import evaluate
 from tsfc.finatinterface import create_element
-from ufl import FiniteElement, RestrictedElement, quadrilateral
+from ufl import quadrilateral
+from finat.ufl import FiniteElement, RestrictedElement
 
 
 def test_issue_274():
diff --git a/tests/test_underintegration.py b/tests/test_underintegration.py
index 05372261..24221e05 100644
--- a/tests/test_underintegration.py
+++ b/tests/test_underintegration.py
@@ -3,9 +3,9 @@
 import numpy
 import pytest
 
-from ufl import (Mesh, FunctionSpace, FiniteElement, VectorElement,
-                 TestFunction, TrialFunction, TensorProductCell, dx,
+from ufl import (Mesh, FunctionSpace, TestFunction, TrialFunction, TensorProductCell, dx,
                  action, interval, quadrilateral, dot, grad)
+from finat.ufl import FiniteElement, VectorElement
 
 from FIAT import ufc_cell
 from FIAT.quadrature import GaussLobattoLegendreQuadratureLineRule, GaussLegendreQuadratureLineRule
diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py
index fb4e9cef..285d0eaa 100644
--- a/tsfc/finatinterface.py
+++ b/tsfc/finatinterface.py
@@ -25,6 +25,7 @@
 import FIAT
 import finat
 import ufl
+import finat.ufl
 
 
 __all__ = ("as_fiat_cell", "create_base_element",
@@ -110,9 +111,9 @@ def convert(element, **kwargs):
 
 
 # Base finite elements first
-@convert.register(ufl.FiniteElement)
+@convert.register(finat.ufl.FiniteElement)
 def convert_finiteelement(element, **kwargs):
-    cell = as_fiat_cell(element.cell())
+    cell = as_fiat_cell(element.cell)
     if element.family() == "Quadrature":
         degree = element.degree()
         scheme = element.quadrature_scheme()
@@ -121,14 +122,14 @@ def convert_finiteelement(element, **kwargs):
 
         return finat.make_quadrature_element(cell, degree, scheme), set()
     lmbda = supported_elements[element.family()]
-    if element.family() == "Real" and element.cell().cellname() in {"quadrilateral", "hexahedron"}:
+    if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}:
         lmbda = None
-        element = ufl.FiniteElement("DQ", element.cell(), 0)
+        element = finat.ufl.FiniteElement("DQ", element.cell, 0)
     if lmbda is None:
-        if element.cell().cellname() == "quadrilateral":
+        if element.cell.cellname() == "quadrilateral":
             # Handle quadrilateral short names like RTCF and RTCE.
             element = element.reconstruct(cell=quadrilateral_tpc)
-        elif element.cell().cellname() == "hexahedron":
+        elif element.cell.cellname() == "hexahedron":
             # Handle hexahedron short names like NCF and NCE.
             element = element.reconstruct(cell=hexahedron_tpc)
         else:
@@ -138,7 +139,7 @@ def convert_finiteelement(element, **kwargs):
         return finat.FlattenedDimensions(finat_elem), deps
 
     kind = element.variant()
-    is_interval = element.cell().cellname() == 'interval'
+    is_interval = element.cell.cellname() == 'interval'
     if kind is None:
         kind = 'spectral' if is_interval else 'equispaced'  # default variant
 
@@ -164,7 +165,7 @@ def convert_finiteelement(element, **kwargs):
             deps = {"shift_axes", "restriction"}
             return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction), deps
         else:
-            raise ValueError("Variant %r not supported on %s" % (kind, element.cell()))
+            raise ValueError("Variant %r not supported on %s" % (kind, element.cell))
     elif element.family() in {"Raviart-Thomas", "Nedelec 1st kind H(curl)",
                               "Brezzi-Douglas-Marini", "Nedelec 2nd kind H(curl)"}:
         lmbda = partial(lmbda, variant=element.variant())
@@ -188,96 +189,96 @@ def convert_finiteelement(element, **kwargs):
             deps = {"shift_axes", "restriction"}
             return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction, continuous=False), deps
         else:
-            raise ValueError("Variant %r not supported on %s" % (kind, element.cell()))
+            raise ValueError("Variant %r not supported on %s" % (kind, element.cell))
     elif element.family() == ["DPC", "DPC L2"]:
-        if element.cell().geometric_dimension() == 2:
+        if element.cell.geometric_dimension() == 2:
             element = element.reconstruct(cell=ufl.cell.hypercube(2))
-        elif element.cell().geometric_dimension() == 3:
+        elif element.cell.geometric_dimension() == 3:
             element = element.reconstruct(cell=ufl.cell.hypercube(3))
     elif element.family() == "S":
-        if element.cell().geometric_dimension() == 2:
+        if element.cell.geometric_dimension() == 2:
             element = element.reconstruct(cell=ufl.cell.hypercube(2))
-        elif element.cell().geometric_dimension() == 3:
+        elif element.cell.geometric_dimension() == 3:
             element = element.reconstruct(cell=ufl.cell.hypercube(3))
 
     return lmbda(cell, element.degree()), set()
 
 
 # Element modifiers and compound element types
-@convert.register(ufl.BrokenElement)
+@convert.register(finat.ufl.BrokenElement)
 def convert_brokenelement(element, **kwargs):
     finat_elem, deps = _create_element(element._element, **kwargs)
     return finat.DiscontinuousElement(finat_elem), deps
 
 
-@convert.register(ufl.EnrichedElement)
+@convert.register(finat.ufl.EnrichedElement)
 def convert_enrichedelement(element, **kwargs):
     elements, deps = zip(*[_create_element(elem, **kwargs)
                            for elem in element._elements])
     return finat.EnrichedElement(elements), set.union(*deps)
 
 
-@convert.register(ufl.NodalEnrichedElement)
+@convert.register(finat.ufl.NodalEnrichedElement)
 def convert_nodalenrichedelement(element, **kwargs):
     elements, deps = zip(*[_create_element(elem, **kwargs)
                            for elem in element._elements])
     return finat.NodalEnrichedElement(elements), set.union(*deps)
 
 
-@convert.register(ufl.MixedElement)
+@convert.register(finat.ufl.MixedElement)
 def convert_mixedelement(element, **kwargs):
     elements, deps = zip(*[_create_element(elem, **kwargs)
-                           for elem in element.sub_elements()])
+                           for elem in element.sub_elements])
     return finat.MixedElement(elements), set.union(*deps)
 
 
-@convert.register(ufl.VectorElement)
-@convert.register(ufl.TensorElement)
+@convert.register(finat.ufl.VectorElement)
+@convert.register(finat.ufl.TensorElement)
 def convert_tensorelement(element, **kwargs):
-    inner_elem, deps = _create_element(element.sub_elements()[0], **kwargs)
-    shape = element.reference_value_shape()
+    inner_elem, deps = _create_element(element.sub_elements[0], **kwargs)
+    shape = element.reference_value_shape
     shape = shape[:len(shape) - len(inner_elem.value_shape)]
     shape_innermost = kwargs["shape_innermost"]
     return (finat.TensorFiniteElement(inner_elem, shape, not shape_innermost),
             deps | {"shape_innermost"})
 
 
-@convert.register(ufl.TensorProductElement)
+@convert.register(finat.ufl.TensorProductElement)
 def convert_tensorproductelement(element, **kwargs):
-    cell = element.cell()
+    cell = element.cell
     if type(cell) is not ufl.TensorProductCell:
         raise ValueError("TensorProductElement not on TensorProductCell?")
     shift_axes = kwargs["shift_axes"]
     dim_offset = 0
     elements = []
     deps = set()
-    for elem in element.sub_elements():
+    for elem in element.sub_elements:
         kwargs["shift_axes"] = shift_axes + dim_offset
-        dim_offset += elem.cell().topological_dimension()
+        dim_offset += elem.cell.topological_dimension()
         finat_elem, ds = _create_element(elem, **kwargs)
         elements.append(finat_elem)
         deps.update(ds)
     return finat.TensorProductElement(elements), deps
 
 
-@convert.register(ufl.HDivElement)
+@convert.register(finat.ufl.HDivElement)
 def convert_hdivelement(element, **kwargs):
     finat_elem, deps = _create_element(element._element, **kwargs)
     return finat.HDivElement(finat_elem), deps
 
 
-@convert.register(ufl.HCurlElement)
+@convert.register(finat.ufl.HCurlElement)
 def convert_hcurlelement(element, **kwargs):
     finat_elem, deps = _create_element(element._element, **kwargs)
     return finat.HCurlElement(finat_elem), deps
 
 
-@convert.register(ufl.WithMapping)
+@convert.register(finat.ufl.WithMapping)
 def convert_withmapping(element, **kwargs):
     return _create_element(element.wrapee, **kwargs)
 
 
-@convert.register(ufl.RestrictedElement)
+@convert.register(finat.ufl.RestrictedElement)
 def convert_restrictedelement(element, **kwargs):
     finat_elem, deps = _create_element(element._element, **kwargs)
     return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps
@@ -323,7 +324,7 @@ def _create_element(ufl_element, **kwargs):
             return finat_element, set(param for param, value in key)
 
     # Convert if cache miss
-    if ufl_element.cell() is None:
+    if ufl_element.cell is None:
         raise ValueError("Don't know how to build element when cell is not given")
 
     finat_element, deps = convert(ufl_element, **kwargs)
diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py
index db7cf53c..67f7dac8 100644
--- a/tsfc/kernel_interface/common.py
+++ b/tsfc/kernel_interface/common.py
@@ -467,7 +467,7 @@ def prepare_coefficient(coefficient, name, interior_facet=False):
 
     if coefficient.ufl_element().family() == 'Real':
         # Constant
-        value_size = coefficient.ufl_element().value_size()
+        value_size = coefficient.ufl_element().value_size
         expression = gem.reshape(gem.Variable(name, (value_size,)),
                                  coefficient.ufl_shape)
         return expression
diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py
index f538acfb..e51fd239 100644
--- a/tsfc/kernel_interface/firedrake_loopy.py
+++ b/tsfc/kernel_interface/firedrake_loopy.py
@@ -2,7 +2,8 @@
 from collections import namedtuple, OrderedDict
 from functools import partial
 
-from ufl import Coefficient, MixedElement as ufl_MixedElement, FunctionSpace, FiniteElement
+from ufl import Coefficient, FunctionSpace
+from finat.ufl import MixedElement as ufl_MixedElement, FiniteElement
 
 import gem
 from gem.flop_count import count_flops
@@ -321,7 +322,7 @@ def set_coefficients(self, integral_data, form_data):
                         k += 1
                     else:
                         self.coefficient_split[coefficient] = []
-                        for j, element in enumerate(coefficient.ufl_element().sub_elements()):
+                        for j, element in enumerate(coefficient.ufl_element().sub_elements):
                             c = Coefficient(FunctionSpace(coefficient.ufl_domain(), element))
                             self.coefficient_split[coefficient].append(c)
                             self._coefficient(c, f"w_{k}")
diff --git a/tsfc/ufl_utils.py b/tsfc/ufl_utils.py
index 21f94d37..0cec4dbf 100644
--- a/tsfc/ufl_utils.py
+++ b/tsfc/ufl_utils.py
@@ -22,9 +22,10 @@
 from ufl.geometry import Jacobian, JacobianDeterminant, JacobianInverse
 from ufl.classes import (Abs, Argument, CellOrientation, Coefficient,
                          ComponentTensor, Expr, FloatValue, Division,
-                         MixedElement, MultiIndex, Product,
+                         MultiIndex, Product,
                          ScalarValue, Sqrt, Zero, CellVolume, FacetArea)
 from ufl.domain import extract_unique_domain
+from finat.ufl import MixedElement
 from ufl.utils.sorting import sorted_by_count
 
 from gem.node import MemoizerArg
@@ -202,7 +203,7 @@ def modified_terminal(self, o):
             # Apply terminal modifiers onto the subcoefficient
             component = construct_modified_terminal(mt, subcoeff)
             # Collect components of the subcoefficient
-            for alpha in numpy.ndindex(subcoeff.ufl_element().reference_value_shape()):
+            for alpha in numpy.ndindex(subcoeff.ufl_element().reference_value_shape):
                 # New modified terminal: component[alpha + beta]
                 components.append(component[alpha + beta])
         # Repack derivative indices to shape
@@ -363,7 +364,7 @@ def apply_mapping(expression, element, domain):
          advertise the pullback to apply.
     :arg domain: Optional domain to provide in case expression does
          not contain a domain (used for constructing geometric quantities).
-    :returns: A new UFL expression with shape element.reference_value_shape()
+    :returns: A new UFL expression with shape element.reference_value_shape
     :raises NotImplementedError: If we don't know how to apply the
         inverse of the pullback.
     :raises ValueError: If we get shape mismatches.
@@ -411,8 +412,8 @@ def apply_mapping(expression, element, domain):
         mesh = domain
     if domain is not None and mesh != domain:
         raise NotImplementedError("Multiple domains not supported")
-    if expression.ufl_shape != element.value_shape():
-        raise ValueError(f"Mismatching shapes, got {expression.ufl_shape}, expected {element.value_shape()}")
+    if expression.ufl_shape != element.value_shape:
+        raise ValueError(f"Mismatching shapes, got {expression.ufl_shape}, expected {element.value_shape}")
     mapping = element.mapping().lower()
     if mapping == "identity":
         rexpression = expression
@@ -447,11 +448,11 @@ def apply_mapping(expression, element, domain):
         # We're going to apply the inverse of the physical to
         # reference space mapping.
         fcm = element.flattened_sub_element_mapping()
-        sub_elem = element.sub_elements()[0]
+        sub_elem = element.sub_elements[0]
         shape = expression.ufl_shape
         flat = ufl.as_vector([expression[i] for i in numpy.ndindex(shape)])
-        vs = sub_elem.value_shape()
-        rvs = sub_elem.reference_value_shape()
+        vs = sub_elem.value_shape
+        rvs = sub_elem.reference_value_shape
         seen = set()
         rpieces = []
         gm = int(numpy.prod(vs, dtype=int))
@@ -469,11 +470,11 @@ def apply_mapping(expression, element, domain):
             # Concatenate with the other pieces
             rpieces.extend([piece[idx] for idx in numpy.ndindex(rvs)])
         # And reshape
-        rexpression = as_tensor(numpy.asarray(rpieces).reshape(element.reference_value_shape()))
+        rexpression = as_tensor(numpy.asarray(rpieces).reshape(element.reference_value_shape))
     else:
-        raise NotImplementedError(f"Don't know how to handle mapping type {mapping} for expression of rank {element.value_shape()}")
-    if rexpression.ufl_shape != element.reference_value_shape():
-        raise ValueError(f"Mismatching reference shapes, got {rexpression.ufl_shape} expected {element.reference_value_shape()}")
+        raise NotImplementedError(f"Don't know how to handle mapping type {mapping} for expression of rank {element.value_shape}")
+    if rexpression.ufl_shape != element.reference_value_shape:
+        raise ValueError(f"Mismatching reference shapes, got {rexpression.ufl_shape} expected {element.reference_value_shape}")
     return rexpression