diff --git a/.github/workflows/test_suite_main.yml b/.github/workflows/test_suite_main.yml
index 532108b..8329951 100644
--- a/.github/workflows/test_suite_main.yml
+++ b/.github/workflows/test_suite_main.yml
@@ -3,7 +3,8 @@ name: 'Run Animate test suite (main)'
 on:
   # Run test suite whenever the main branch is updated
   push:
-    branches: [main]
+    branches:
+      - main
     paths:
       - '**.py'
       - '**.cxx'
@@ -13,7 +14,6 @@ on:
 
   # Run test suite whenever commits are pushed to an open PR
   pull_request:
-    branches: [main]
     paths:
       - '**.py'
       - '**.cxx'
diff --git a/.github/workflows/test_suite_release.yml b/.github/workflows/test_suite_release.yml
index 2b7e436..af664a9 100644
--- a/.github/workflows/test_suite_release.yml
+++ b/.github/workflows/test_suite_release.yml
@@ -3,13 +3,13 @@ name: 'Run Animate test suite (release)'
 on:
   # Run test suite whenever the release branch is updated
   push:
-    branches: [release]
+    branches:
+      - release
     paths:
       - '.github/workflows/test_suite_release.yml'
 
   # Run test suite whenever commits that change this workflow are pushed to an open PR
   pull_request:
-    branches: [release]
     paths:
       - '.github/workflows/test_suite_release.yml'
 
diff --git a/animate/adapt.py b/animate/adapt.py
index 601164a..e9ed300 100644
--- a/animate/adapt.py
+++ b/animate/adapt.py
@@ -1,4 +1,5 @@
 import abc
+import gc
 import os
 from functools import cached_property
 from shutil import rmtree
@@ -97,10 +98,9 @@ def adapted_mesh(self):
         :rtype: :class:`firedrake.mesh.MeshGeometry`
         """
         self.metric.enforce_spd(restrict_sizes=True, restrict_anisotropy=True)
-        size = self.metric.dat.dataset.layout_vec.getSizes()
-        data = self.metric.dat._data[: size[0]]
+        data = self.metric.dat.data_ro_with_halos
         v = PETSc.Vec().createWithArray(
-            data, size=size, bsize=self.metric.dat.cdim, comm=self.mesh.comm
+            data, size=data.size, bsize=self.metric.dat.cdim, comm=COMM_SELF
         )
         reordered = to_petsc_local_numbering(v, self.metric.function_space())
         v.destroy()
@@ -176,7 +176,7 @@ def adapt(mesh, *metrics, name=None, serialise=None, remove_checkpoints=True):
     """
     nprocs = COMM_WORLD.size
 
-    dim = mesh.topological_dimension()
+    dim = mesh.topological_dimension
     if serialise is None:
         serialise = nprocs > 1 and dim != 3
     elif not serialise and dim != 3:
@@ -193,13 +193,20 @@ def adapt(mesh, *metrics, name=None, serialise=None, remove_checkpoints=True):
         chk_fpath = os.path.join(chk_dir, "adapted_mesh_checkpoint.h5")
         metric_name = "tmp_metric"
         save_checkpoint(chk_fpath, metric, metric_name)
+        # Ensure all processes are finished writing
+        COMM_WORLD.barrier()
 
         if COMM_WORLD.rank == 0:
             metric0 = load_checkpoint(chk_fpath, mesh.name, metric_name, comm=COMM_SELF)
             adaptor0 = MetricBasedAdaptor(metric0._mesh, metric0, name=name)
             with fchk.CheckpointFile(chk_fpath, "w", comm=COMM_SELF) as chk:
                 chk.save_mesh(adaptor0.adapted_mesh)
+        # Ensure rank 0 is finished writing
         COMM_WORLD.barrier()
+        # Garbage collection might be called at different times on different ranks due
+        # to diverging paths, which appears to stall in final cleanup on Python
+        # system exit. Ensure everything is in-sync again at this point
+        gc.collect()
 
         # In parallel, load from the checkpoint
         if not os.path.exists(chk_fpath):
diff --git a/animate/cython/numbering.pyx b/animate/cython/numbering.pyx
index 3abaee2..c80b124 100644
--- a/animate/cython/numbering.pyx
+++ b/animate/cython/numbering.pyx
@@ -26,25 +26,27 @@ cdef extern from "petscis.h" nogil:
 def to_petsc_local_numbering(PETSc.Vec vec, V):
     """
     Reorder a PETSc Vec corresponding to a Firedrake Function w.r.t.
-    the PETSc natural numbering.
+    the PETSc natural numbering, i.e. the numbering consistent with
+    that of the DMPlex topological points.
 
-    :arg vec: the PETSc Vec to reorder; must be a global vector
+    :arg vec: the PETSc Vec to reorder; must be a local vector, i.e.
+              a sequential vector that includes all (owned and halo) DoFs
     :arg V: the FunctionSpace of the Function which the Vec comes from
     :ret out: a copy of the Vec, ordered with the PETSc natural numbering
     """
-    cdef int dim, idx, start, end, p, d, k
+    cdef int dim, idx, lsize, p, d, k
     cdef PetscInt dof, off
     cdef PETSc.Vec out
     cdef PETSc.Section section
     cdef np.ndarray varray, oarray
 
-    section = V.dm.getGlobalSection()
+    section = V.dm.getLocalSection()
     out = vec.duplicate()
     varray = vec.array_r
     oarray = out.array
     dim = V.value_size
     idx = 0
-    start, end = vec.getOwnershipRange()
+    lsize = vec.getSize()
     for p in range(*section.getChart()):
         CHKERR(PetscSectionGetDof(section.sec, p, &dof))
         if dof > 0:
@@ -53,7 +55,12 @@ def to_petsc_local_numbering(PETSc.Vec vec, V):
             off *= dim
             for d in range(dof):
                 for k in range(dim):
-                    oarray[idx] = varray[off + dim * d + k - start]
+                    oarray[idx] = varray[off + dim * d + k]
                     idx += 1
-    assert idx == (end - start)
+    if idx != lsize:
+        raise ValueError(
+           f"Number of local section entries not the same as vector size"
+           f"({idx} vs. {lsize}). Need to provide local vector including halo DoFs."
+        )
+        
     return out
diff --git a/animate/interpolation.py b/animate/interpolation.py
index 75423ef..a735939 100644
--- a/animate/interpolation.py
+++ b/animate/interpolation.py
@@ -324,7 +324,7 @@ def clement_interpolant(source, target_space=None, boundary=False):
     if rank not in (0, 1, 2):
         raise ValueError(f"Rank-{rank + 1} tensors are not supported.")
     mesh = Vs.mesh()
-    dim = mesh.topological_dimension()
+    dim = mesh.topological_dimension
 
     # Process target space
     Vt = target_space
diff --git a/animate/math.py b/animate/math.py
index 403723a..4d492b1 100644
--- a/animate/math.py
+++ b/animate/math.py
@@ -73,7 +73,7 @@ def construct_basis(vector, normalise=True):
         if not isinstance(vector, ufl.core.expr.Expr):
             raise TypeError(f"Expected UFL Expr, not '{type(vector)}'.")
         as_vector = ufl.as_vector
-        dim = ufl.domain.extract_unique_domain(vector).topological_dimension()
+        dim = ufl.domain.extract_unique_domain(vector).topological_dimension
 
     if dim not in (2, 3):
         raise ValueError(f"Dimension {dim} not supported.")
diff --git a/animate/metric.py b/animate/metric.py
index e30d08d..b449dbd 100644
--- a/animate/metric.py
+++ b/animate/metric.py
@@ -89,7 +89,7 @@ def __init__(self, function_space, *args, **kwargs):
         # Check that we have an appropriate tensor P1 function
         fs = self.function_space()
         mesh = fs.mesh()
-        tdim = mesh.topological_dimension()
+        tdim = mesh.topological_dimension
         if tdim not in (2, 3):
             raise ValueError(f"Riemannian metric should be 2D or 3D, not {tdim}D.")
         if isinstance(fs.dof_count, Iterable):
@@ -303,7 +303,7 @@ def _set_plex_coordinates(self):
         consistent.
         """
         entity_dofs = np.zeros(self._tdim + 1, dtype=np.int32)
-        entity_dofs[0] = self._mesh.geometric_dimension()
+        entity_dofs[0] = self._mesh.geometric_dimension
         coord_section = self._mesh.create_section(entity_dofs)[0]
         # NOTE: section doesn't have any fields, but PETSc assumes it to have one
         coord_dm = self._plex.getCoordinateDM()
@@ -349,6 +349,14 @@ def compute_hessian(self, field, method="mixed_L2", **kwargs):
         for the gradient recovery. The target space for the Hessian recovery is
         inherited from the metric itself.
         """
+        if len(field.ufl_shape) > 0:
+            value_error = (
+                "RiemannianMetric.compute_hessian only accepts scalar fields. To "
+                "recover Hessians of higher rank fields, call the method on separate "
+                "RiemannianMetrics for each component and combine them with "
+                "RiemannianMetric.combine."
+            )
+            raise ValueError(value_error)
         if method == "L2":
             gradient = recover_gradient_l2(
                 field, target_space=kwargs.get("target_space")
@@ -518,7 +526,7 @@ def interp(f):
         if not np.isclose(_a_max, 1.0) and _a_max < 1.0:
             raise ValueError(f"Encountered a_max value smaller than unity: {_a_max}.")
 
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         boundary_tag = self._variable_parameters.get("dm_plex_metric_boundary_tag")
         if boundary_tag is None:
             node_set = self.function_space().node_set
@@ -621,25 +629,28 @@ def intersect(self, *metrics):
                     "Cannot intersect metrics with different function spaces."
                 )
 
-        # Intersect the metrics recursively one at a time
         if len(metrics) == 0:
-            pass
-        elif len(metrics) == 1:
-            v1 = self._create_from_array(self.dat.data_with_halos)
-            v2 = self._create_from_array(metrics[0].dat.data_ro_with_halos)
-            vout = self._create_from_array(np.zeros_like(self.dat.data_with_halos))
-
-            # Compute the intersection on the PETSc level
-            self._plex.metricIntersection2(v1, v2, vout)
-
-            # Assign to the output of the intersection
-            size = np.shape(self.dat.data_with_halos)
-            self.dat.data_with_halos[:] = np.reshape(vout.array, size)
-            v2.destroy()
-            v1.destroy()
-            vout.destroy()
-        else:
+            return self
+
+        # Intersect the metrics recursively, starting with metrics[0]
+        v1 = self._create_from_array(self.dat.data_with_halos)
+        v2 = self._create_from_array(metrics[0].dat.data_ro_with_halos)
+        vout = self._create_from_array(np.zeros_like(self.dat.data_with_halos))
+
+        # Compute the intersection on the PETSc level
+        self._plex.metricIntersection2(v1, v2, vout)
+
+        # Assign to the output of the intersection
+        size = np.shape(self.dat.data_with_halos)
+        self.dat.data_with_halos[:] = np.reshape(vout.array, size)
+        v2.destroy()
+        v1.destroy()
+        vout.destroy()
+
+        # Intersect with the remaining metrics
+        if len(metrics) > 1:
             self.intersect(*metrics[1:])
+
         return self
 
     @PETSc.Log.EventDecorator()
@@ -734,7 +745,7 @@ def compute_eigendecomposition(self, reorder=False):
         mesh = V_ten.mesh()
         fe = (V_ten.ufl_element().family(), V_ten.ufl_element().degree())
         V_vec = firedrake.VectorFunctionSpace(mesh, *fe)
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         evectors, evalues = firedrake.Function(V_ten), firedrake.Function(V_vec)
         if reorder:
             name = "get_reordered_eigendecomposition"
@@ -788,7 +799,7 @@ def assemble_eigendecomposition(self, evectors, evalues):
                 "Mismatching finite element space degrees:"
                 f" {fe_ten.degree()} vs. {fe_vec.degree()}."
             )
-        dim = V_ten.mesh().topological_dimension()
+        dim = V_ten.mesh().topological_dimension
         op2.par_loop(
             get_metric_kernel("set_eigendecomposition", dim),
             V_ten.node_set,
@@ -848,7 +859,7 @@ def density_and_quotients(self, reorder=False):
         fs_ten = self.function_space()
         mesh = fs_ten.mesh()
         fe = (fs_ten.ufl_element().family(), fs_ten.ufl_element().degree())
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         evectors, evalues = self.compute_eigendecomposition(reorder=reorder)
 
         # Extract density and quotients
@@ -886,7 +897,7 @@ def compute_isotropic_metric(
         mesh = ufl.domain.extract_unique_domain(error_indicator)
         if mesh != self.function_space().mesh():
             raise ValueError("Cannot use an error indicator from a different mesh.")
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
 
         # Interpolate P0 indicators into P1 space
         if interpolant == "Clement":
@@ -980,7 +991,7 @@ def compute_anisotropic_dwr_metric(
         mesh = ufl.domain.extract_unique_domain(error_indicator)
         if mesh != self.function_space().mesh():
             raise ValueError("Cannot use an error indicator from a different mesh.")
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         if convergence_rate < 1.0:
             raise ValueError(
                 f"Convergence rate must be at least one, not {convergence_rate}."
@@ -1133,7 +1144,7 @@ def determine_metric_complexity(H_interior, H_boundary, target, p, **kwargs):
     :returns: unique solution of algebraic problem
     :rtype: :class:`float`
     """
-    d = H_interior.function_space().mesh().topological_dimension()
+    d = H_interior.function_space().mesh().topological_dimension
     if d not in (2, 3):
         raise ValueError(f"Spatial dimension {d} not supported.")
     if np.isinf(p):
@@ -1183,7 +1194,7 @@ def intersect_on_boundary(*metrics, boundary_tag="on_boundary"):
     n = len(metrics)
     assert n > 0, "Nothing to combine"
     fs = metrics[0].function_space()
-    dim = fs.mesh().topological_dimension()
+    dim = fs.mesh().topological_dimension
     if dim not in (2, 3):
         raise ValueError(
             f"Spatial dimension {dim} not supported. Must be either 2 or 3."
diff --git a/animate/quality.py b/animate/quality.py
index c05402d..2cd6e29 100644
--- a/animate/quality.py
+++ b/animate/quality.py
@@ -8,10 +8,10 @@
 import ufl
 from firedrake.__future__ import interpolate
 from firedrake.petsc import PETSc
+from petsctools import get_petsc_dirs
 from pyop2 import op2
-from pyop2.utils import get_petsc_dir
 
-petsc_dirs = get_petsc_dir()
+petsc_dirs = get_petsc_dirs()
 include_dir = ["%s/include/eigen3" % petsc_dirs[-1]]
 
 __all__ = ["QualityMeasure"]
@@ -60,7 +60,7 @@ def __init__(self, mesh, metric=None, python=False):
         self.mesh = mesh
         self.metric = metric
         self.python = python
-        self.dim = mesh.topological_dimension()
+        self.dim = mesh.topological_dimension
         self.coords = mesh.coordinates
         self.P0 = firedrake.FunctionSpace(mesh, "DG", 0)
         src_dir = os.path.join(os.path.dirname(__file__), "cxx")
diff --git a/animate/recovery.py b/animate/recovery.py
index 5b1e457..e5ca1b5 100644
--- a/animate/recovery.py
+++ b/animate/recovery.py
@@ -132,7 +132,7 @@ def recover_boundary_hessian(f, method="Clement", target_space=None, **kwargs):
     """
 
     mesh = ufl.domain.extract_unique_domain(f)
-    d = mesh.topological_dimension()
+    d = mesh.topological_dimension
     assert d in (2, 3)
 
     # Apply Gram-Schmidt to get tangent vectors
diff --git a/test/test_checkpoint.py b/test/test_checkpoint.py
index f687f35..9f4e110 100644
--- a/test/test_checkpoint.py
+++ b/test/test_checkpoint.py
@@ -7,7 +7,6 @@
 import ufl
 from firedrake.function import Function
 from firedrake.functionspace import FunctionSpace, TensorFunctionSpace
-from firedrake.mesh import _generate_default_mesh_topology_name
 from test_setup import uniform_mesh
 
 from animate.checkpointing import get_checkpoint_dir, load_checkpoint, save_checkpoint
@@ -57,8 +56,8 @@ def test_save(self, filename="test_save.h5", metric=None):
         metric = metric or self.metric
         save_checkpoint(fpath, metric)
         self.assertTrue(os.path.exists(fpath))
-        mesh_name = metric._mesh.name
-        topology_name = _generate_default_mesh_topology_name(mesh_name)
+        mesh = ufl.domain.extract_unique_domain(metric)
+        topology_name = mesh.topology_dm.name
         with h5py.File(fpath, "r") as h5:
             self.assertTrue("topologies" in h5)
             self.assertTrue(topology_name in h5["topologies"].keys())
diff --git a/test/test_interpolation.py b/test/test_interpolation.py
index 67ad6d6..a02166d 100644
--- a/test/test_interpolation.py
+++ b/test/test_interpolation.py
@@ -66,7 +66,7 @@ def get_space(self, rank, family, degree):
         elif rank == 1:
             return VectorFunctionSpace(self.mesh, family, degree)
         else:
-            shape = tuple(rank * [self.mesh.topological_dimension()])
+            shape = tuple(rank * [self.mesh.topological_dimension])
             return TensorFunctionSpace(self.mesh, family, degree, shape=shape)
 
     def analytic(self, rank):
diff --git a/test/test_metric.py b/test/test_metric.py
index 9faeaa3..e0b3bf7 100644
--- a/test/test_metric.py
+++ b/test/test_metric.py
@@ -194,7 +194,33 @@ class TestHessianMetric(MetricTestCase):
     Unit tests for the :meth:`compute_hessian` method of :class:`RiemannianMetric`.
     """
 
-    def test_bowl(self, dim=2, places=7):
+    def test_vector_rank_error(self):
+        mesh = uniform_mesh(2)
+        P1_vec = VectorFunctionSpace(mesh, "CG", 1)
+        P1_ten = TensorFunctionSpace(mesh, "CG", 1)
+        with self.assertRaises(ValueError) as cm:
+            RiemannianMetric(P1_ten).compute_hessian(Function(P1_vec))
+        msg = (
+            "RiemannianMetric.compute_hessian only accepts scalar fields. To "
+            "recover Hessians of higher rank fields, call the method on separate "
+            "RiemannianMetrics for each component and combine them with "
+            "RiemannianMetric.combine."
+        )
+        self.assertEqual(str(cm.exception), msg)
+
+    def test_tensor_rank_error(self):
+        P1_ten = TensorFunctionSpace(uniform_mesh(2), "CG", 1)
+        with self.assertRaises(ValueError) as cm:
+            RiemannianMetric(P1_ten).compute_hessian(Function(P1_ten))
+        msg = (
+            "RiemannianMetric.compute_hessian only accepts scalar fields. To "
+            "recover Hessians of higher rank fields, call the method on separate "
+            "RiemannianMetrics for each component and combine them with "
+            "RiemannianMetric.combine."
+        )
+        self.assertEqual(str(cm.exception), msg)
+
+    def test_bowl(self, dim=2, places=6):
         mesh = uniform_mesh(dim, 4, recentre=True)
         P1_ten = TensorFunctionSpace(mesh, "CG", 1)
         metric = RiemannianMetric(P1_ten).compute_hessian(bowl(*mesh.coordinates))
@@ -282,16 +308,37 @@ def test_multiple_intersect(self, dim):
         mesh = uniform_mesh(dim, 1)
         P1_ten = TensorFunctionSpace(mesh, "CG", 1)
 
-        metric1 = uniform_metric(P1_ten, 100.0)
-        metric2 = uniform_metric(P1_ten, 40.0)
+        metric1 = uniform_metric(P1_ten, 40.0)
+        metric2 = uniform_metric(P1_ten, 100.0)
         metric3 = uniform_metric(P1_ten, 20.0)
-        expected = metric1
+        expected = metric2
 
         metric = RiemannianMetric(P1_ten)
         metric.assign(metric1)
         metric.intersect(metric2, metric3)
         self.assertAlmostMatching(metric, expected)
 
+    @parameterized.expand([[2], [3]])
+    def test_multiple_variable_intersect(self, dim):
+        mesh = uniform_mesh(dim, 1)
+        P1_ten = TensorFunctionSpace(mesh, "CG", 1)
+
+        expected = RiemannianMetric(P1_ten)
+        expected.interpolate(ufl.Identity(dim))
+        metrics = []
+        for i in range(P1_ten.node_count):
+            metric = RiemannianMetric(P1_ten)
+            metric.interpolate(ufl.Identity(dim))
+            # change to (i+1)I at node i only
+            metric.dat.data[i] *= float(i+1)
+            expected.dat.data[i] *= float(i+1)
+            metrics.append(metric)
+
+        metric = RiemannianMetric(P1_ten)
+        metric.assign(metrics[0])
+        metric.intersect(*metrics[1:])
+        self.assertAlmostMatching(metric, expected)
+
 
 class TestNormalisation(MetricTestCase):
     """
diff --git a/test/test_numbering.py b/test/test_numbering.py
index 9d859c3..130411a 100644
--- a/test/test_numbering.py
+++ b/test/test_numbering.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pytest
 from animate.cython.numbering import to_petsc_local_numbering
+from firedrake.petsc import PETSc
 
 
 @pytest.fixture
@@ -35,13 +36,47 @@ def test_is_permutation(function_space):
     """Test that to_petsc_local_numbering is indeed a permutation."""
     f = fd.Function(function_space)
 
+    comm = function_space.comm
+    # unique float 0<=x<1 for each rank:
+    rank_fraction = comm.rank / (comm.size + 1)
+
+    # Initialise owned and in particular halo DoFs, so we can check they
+    # have been updated later on
+    f.assign(-1)
+
     # Fill the vector with arbitrary values
+    # f.dat.vec below is a global Vec, so owned entries only
+    # halos values should be updated automatically coming out of the context
     with f.dat.vec as vec:
-        vec.array[:] = np.arange(vec.size)
+        vec.array[:] = np.arange(vec.sizes[0]) + rank_fraction
+
+    # Check that this is the case:
+    owned = f.dat.data_ro.flatten()
+    halos = f.dat.data_ro_with_halos[owned.size:].flatten()
+    np.testing.assert_equal(owned, np.arange(owned.size) + rank_fraction)
+    if halos.size > 0:
+        # all halo enties should be set and come from different rank
+        assert all(halos >= 0.0)
+        assert all(np.modf(halos)[0] != rank_fraction)
 
     # Verify that reordering according to the PETSc numbering is a permutation
-    with f.dat.vec_ro as vec:
-        reordered_vec = to_petsc_local_numbering(vec, function_space)
-        assert reordered_vec is not None
-        assert reordered_vec.size == vec.size
-        assert sorted(reordered_vec.array) == sorted(vec.array)
+    data = f.dat.data_ro_with_halos
+    lvec = PETSc.Vec().createWithArray(
+        data, size=data.size, bsize=f.dat.cdim, comm=fd.COMM_SELF
+    )
+    reordered_vec = to_petsc_local_numbering(lvec, function_space)
+    assert reordered_vec is not None
+    assert reordered_vec.size == lvec.size
+    assert sorted(reordered_vec.array) == sorted(lvec.array)
+    lvec.destroy()
+    reordered_vec.destroy()
+
+
+@pytest.mark.parallel(nprocs=2)
+def test_is_permutation_np2(function_space):
+    test_is_permutation(function_space)
+
+
+@pytest.mark.parallel(nprocs=3)
+def test_is_permutation_np3(function_space):
+    test_is_permutation(function_space)
diff --git a/test/test_quality.py b/test/test_quality.py
index 703b94b..133163c 100644
--- a/test/test_quality.py
+++ b/test/test_quality.py
@@ -37,7 +37,7 @@ class TestQuality(unittest.TestCase):
     """
 
     def quality(self, name, mesh, **kwargs):
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         if name == "metric":
             P1_ten = TensorFunctionSpace(mesh, "CG", 1)
             M = Function(P1_ten).interpolate(ufl.Identity(dim))
diff --git a/test/test_recovery.py b/test/test_recovery.py
index 39cd1a6..5c7bf61 100644
--- a/test/test_recovery.py
+++ b/test/test_recovery.py
@@ -15,6 +15,7 @@
 
 from animate.math import construct_basis
 from animate.metric import RiemannianMetric
+from animate.recovery import recover_gradient_l2
 
 # ---------------------------
 # standard tests for pytest
@@ -47,7 +48,7 @@ def test_clement_function_error(self):
         self.assertEqual(str(cm.exception), msg)
 
     def test_clement_space_error(self):
-        f = self.get_func_ones("RT", 1)
+        f = self.get_func_ones("DG", 0)
         with self.assertRaises(ValueError) as cm:
             self.metric.compute_hessian(f, method="Clement")
         msg = (
@@ -69,14 +70,14 @@ def test_clement_degree_error(self):
     def test_l2_projection_function_error(self):
         f = bowl(*ufl.SpatialCoordinate(self.mesh))
         with self.assertRaises(ValueError) as cm:
-            self.metric.compute_hessian(f, method="L2")
+            recover_gradient_l2(f)
         msg = "If a target space is not provided then the input must be a Function."
         self.assertEqual(str(cm.exception), msg)
 
     def test_l2_projection_rank_error(self):
         f = Function(TensorFunctionSpace(self.mesh, "CG", 1))
         with self.assertRaises(ValueError) as cm:
-            self.metric.compute_hessian(f, method="L2")
+            recover_gradient_l2(f)
         msg = (
             "L2 projection can only be used to compute gradients of scalar or vector"
             " Functions, not Functions of rank 2."
@@ -120,7 +121,7 @@ def metric(mesh):
     @staticmethod
     def relative_error(approx, ignore_boundary=False):
         mesh = approx.function_space().mesh()
-        dim = mesh.topological_dimension()
+        dim = mesh.topological_dimension
         P1_ten = TensorFunctionSpace(mesh, "CG", 1)
         identity = Function(P1_ten).interpolate(ufl.Identity(dim))
 
diff --git a/test/test_setup.py b/test/test_setup.py
index 29364fe..42589ac 100644
--- a/test/test_setup.py
+++ b/test/test_setup.py
@@ -49,7 +49,7 @@ def uniform_metric(mesh, a=100.0, metric_parameters=None):
     else:
         function_space = mesh
         mesh = function_space.mesh()
-    dim = mesh.topological_dimension()
+    dim = mesh.topological_dimension
     metric = RiemannianMetric(function_space)
     metric.interpolate(a * ufl.Identity(dim))
     metric.set_parameters(metric_parameters)