diff --git a/firedrake/cofunction.py b/firedrake/cofunction.py
index 2dca78e961..0547ee0705 100644
--- a/firedrake/cofunction.py
+++ b/firedrake/cofunction.py
@@ -349,6 +349,9 @@ def vector(self):
         :class:`Cofunction`"""
         return vector.Vector(self)
 
+    def nodal_dat(self):
+        return op3.HierarchicalArray(self.function_space().nodal_axes, data=self.dat.data_rw_with_halos)
+
     @property
     def node_set(self):
         r"""A :class:`pyop2.types.set.Set` containing the nodes of this
@@ -356,7 +359,7 @@ def node_set(self):
         :class:`.FunctionSpace`\s) more degrees of freedom are stored
         at each node.
         """
-        return self.function_space().node_set
+        return self.function_space().nodes
 
     def ufl_id(self):
         return self.uid
@@ -386,5 +389,10 @@ def __str__(self):
         else:
             return super(Cofunction, self).__str__()
 
-    def cell_node_map(self):
-        return self.function_space().cell_node_map()
+    @property
+    def cell_node_list(self):
+        return self.function_space().cell_node_list
+
+    @property
+    def owned_cell_node_list(self):
+        return self.function_space().owned_cell_node_list
diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx
index dd575b3dca..328c7c38e9 100644
--- a/firedrake/cython/dmcommon.pyx
+++ b/firedrake/cython/dmcommon.pyx
@@ -1224,8 +1224,8 @@ def create_section(mesh, nodes_per_entity, on_base=False, block_size=1):
     if isinstance(dm, PETSc.DMSwarm) and on_base:
         raise NotImplementedError("Vertex Only Meshes cannot be extruded.")
     variable = mesh.variable_layers
-    extruded = mesh.cell_set._extruded
-    extruded_periodic = mesh.cell_set._extruded_periodic
+    extruded = mesh.extruded
+    extruded_periodic = mesh.extruded_periodic
     on_base_ = on_base
     nodes_per_entity = np.asarray(nodes_per_entity, dtype=IntType)
     if variable:
diff --git a/firedrake/cython/mgimpl.pyx b/firedrake/cython/mgimpl.pyx
index bd7f9d13f8..6cc32ddf6b 100644
--- a/firedrake/cython/mgimpl.pyx
+++ b/firedrake/cython/mgimpl.pyx
@@ -64,8 +64,8 @@ def coarse_to_fine_nodes(Vc, Vf, np.ndarray[PetscInt, ndim=2, mode="c"] coarse_t
         PetscInt fine_layer, fine_layers, coarse_layer, coarse_layers, ratio
         bint extruded
 
-    fine_map = Vf.cell_node_map().values
-    coarse_map = Vc.cell_node_map().values
+    fine_map = Vf.owned_cell_node_list
+    coarse_map = Vc.owned_cell_node_list
 
     fine_cell_per_coarse_cell = coarse_to_fine_cells.shape[1]
     extruded = Vc.extruded
@@ -85,7 +85,7 @@ def coarse_to_fine_nodes(Vc, Vf, np.ndarray[PetscInt, ndim=2, mode="c"] coarse_t
     ndof = fine_per_cell * fine_cell_per_coarse_cell
     if extruded:
         ndof *= ratio
-    coarse_to_fine_map = np.full((Vc.dof_dset.total_size,
+    coarse_to_fine_map = np.full((Vc.node_count,
                                   ndof),
                                  -1,
                                  dtype=IntType)
@@ -124,8 +124,8 @@ def fine_to_coarse_nodes(Vf, Vc, np.ndarray[PetscInt, ndim=2, mode="c"] fine_to_
         PetscInt coarse_per_cell, fine_per_cell, coarse_cell, fine_cells
         bint extruded
 
-    fine_map = Vf.cell_node_map().values
-    coarse_map = Vc.cell_node_map().values
+    fine_map = Vf.owned_cell_node_list
+    coarse_map = Vc.owned_cell_node_list
 
     extruded = Vc.extruded
 
@@ -142,7 +142,7 @@ def fine_to_coarse_nodes(Vf, Vc, np.ndarray[PetscInt, ndim=2, mode="c"] fine_to_
     coarse_per_fine = fine_to_coarse_cells.shape[1]
     coarse_per_cell = coarse_map.shape[1]
     fine_per_cell = fine_map.shape[1]
-    fine_to_coarse_map = np.full((Vf.dof_dset.total_size,
+    fine_to_coarse_map = np.full((Vf.node_count,
                                   coarse_per_fine*coarse_per_cell),
                                  -1,
                                  dtype=IntType)
@@ -255,8 +255,8 @@ def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps):
     fdm = mf.topology_dm
     dim = cdm.getDimension()
     nref = 2 ** dim
-    ncoarse = mc.cell_set.size
-    nfine = mf.cell_set.size
+    ncoarse = mc.cells.owned.size
+    nfine = mf.cells.owned.size
     co2n, _ = get_entity_renumbering(cdm, mc._cell_numbering, "cell")
     _, fn2o = get_entity_renumbering(fdm, mf._cell_numbering, "cell")
     coarse_to_fine = np.full((ncoarse, nref), -1, dtype=PETSc.IntType)
@@ -274,7 +274,7 @@ def coarse_to_fine_cells(mc, mf, clgmaps, flgmaps):
         # Need to permute order of co2n so it maps from non-overlapped
         # cells to new cells (these may have changed order).  Need to
         # map all known cells through.
-        idx = np.arange(mc.cell_set.total_size, dtype=PETSc.IntType)
+        idx = np.arange(mc.cells.size, dtype=PETSc.IntType)
         # LocalToGlobal
         co.apply(idx, result=idx)
         # GlobalToLocal
diff --git a/firedrake/function.py b/firedrake/function.py
index 4fcbfeed33..a4f45e4961 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -381,6 +381,10 @@ def vector(self):
         r"""Return a :class:`.Vector` wrapping the data in this :class:`Function`"""
         return vector.Vector(self)
 
+    def nodal_dat(self):
+        return op3.HierarchicalArray(self.function_space().nodal_axes,
+                                     data=self.dat.data_rw_with_halos)
+
     @PETSc.Log.EventDecorator()
     def interpolate(
         self,
diff --git a/firedrake/functionspacedata.py b/firedrake/functionspacedata.py
index f7bd664bbb..237bd38c64 100644
--- a/firedrake/functionspacedata.py
+++ b/firedrake/functionspacedata.py
@@ -103,7 +103,7 @@ def get_node_set(mesh, key):
     node_classes = mesh.node_classes(nodes_per_entity, real_tensorproduct=real_tensorproduct)
     halo = halo_mod.Halo(mesh.topology_dm, global_numbering, comm=mesh.comm)
     node_set = op2.Set(node_classes, halo=halo, comm=mesh.comm)
-    extruded = mesh.cell_set._extruded
+    extruded = mesh.extruded
 
     assert global_numbering.getStorageSize() == node_set.total_size
     if not extruded and node_set.total_size >= (1 << (IntType.itemsize * 8 - 4)):
@@ -211,7 +211,7 @@ def get_boundary_masks(mesh, key, finat_element):
         with points in the closure of point p.  The basis function
         indices are in the index array, starting at section.getOffset(p).
     """
-    if not mesh.cell_set._extruded:
+    if not mesh.extruded:
         return None
     _, kind = key
     assert kind in {"cell", "interior_facet"}
@@ -453,7 +453,7 @@ def __init__(self, mesh, ufl_element):
         self.node_set = node_set
         self.cell_boundary_masks = get_boundary_masks(mesh, (edofs_key, "cell"), finat_element)
         self.interior_facet_boundary_masks = get_boundary_masks(mesh, (edofs_key, "interior_facet"), finat_element)
-        self.extruded = mesh.cell_set._extruded
+        self.extruded = mesh.extruded
         self.mesh = mesh
         # self.global_numbering = global_numbering
 
diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py
index d2faa0ea08..0e6a96806d 100644
--- a/firedrake/functionspaceimpl.py
+++ b/firedrake/functionspaceimpl.py
@@ -102,6 +102,7 @@ def __init__(self, mesh, element, component=None, cargo=None):
         self.cargo = cargo
         self.comm = mesh.comm
         self._comm = mpi.internal_comm(mesh.comm, self)
+        self.extruded = mesh.extruded
 
     @classmethod
     def create(cls, function_space, mesh):
@@ -562,7 +563,7 @@ def __init__(self, mesh, element, name=None):
                 mesh._shared_data_cache[key] = (axes, block_axes)
 
         self.axes = axes
-        self.block_axes = axes
+        self.block_axes = block_axes
 
     # These properties are overridden in ProxyFunctionSpaces, but are
     # provided by FunctionSpace so that we don't have to special case.
@@ -694,10 +695,44 @@ def local_section(self):
     def _cdim(self):
         return self.value_size
 
+    @utils.cached_property
+    def nodes(self):
+        ax = self.block_axes
+        return op3.Axis([op3.AxisComponent((ax.owned.size, ax.size),
+                         "XXX", rank_equal=False)], "nodes", numbering=None, sf=ax.sf)
+
+    @utils.cached_property
+    def nodal_axes(self):
+        return op3.AxisTree.from_iterable([self.nodes, self.value_size])
+
+    @utils.cached_property
+    def cell_node_dat(self):
+        from firedrake.parloops import pack_pyop3_tensor
+        cells = self.mesh().cells
+        # Pass self.sub(0) to get nodes from the scalar version of this function space
+        packed_axes = pack_pyop3_tensor(self.block_axes, self.sub(0), cells.index(include_ghost_points=True), "cell")
+        return packed_axes.tabulated_offsets
+
+    @utils.cached_property
+    def cell_node_map(self):
+        from pyrsistent import freeze
+        return op3.Map({
+            freeze({self.mesh().topology.name: self.mesh().cells.owned.root.component.label}): [
+                op3.TabulatedMapComponent(self.nodes.label, self.nodes.component.label, self.cell_node_dat)
+            ]
+        })
+
     @utils.cached_property
     def cell_node_list(self):
-        r"""A numpy array mapping mesh cells to function space nodes."""
-        return self._shared_data.entity_node_lists[self.mesh().cell_set]
+        r"""A numpy array mapping mesh cells to function space nodes (includes halo)."""
+        cells = self.mesh().cells
+        return self.cell_node_dat.buffer.data_rw_with_halos.reshape((cells.size, -1))
+
+    @utils.cached_property
+    def owned_cell_node_list(self):
+        r"""A numpy array mapping owned mesh cells to function space nodes."""
+        cells = self.mesh().cells
+        return self.cell_node_list[:cells.owned.size]
 
     @utils.cached_property
     def topological(self):
@@ -771,13 +806,13 @@ def node_count(self):
         this process.  If the :class:`FunctionSpace` has :attr:`FunctionSpace.rank` 0, this
         is equal to the :attr:`FunctionSpace.dof_count`, otherwise the :attr:`FunctionSpace.dof_count` is
         :attr:`dim` times the :attr:`node_count`."""
-        return self.node_set.total_size
+        return self.block_axes.size
 
     @utils.cached_property
     def dof_count(self):
         r"""The number of degrees of freedom (includes halo dofs) of this
         function space on this process. Cf. :attr:`FunctionSpace.node_count` ."""
-        return self.node_count*self.value_size
+        return self.axes.size
 
     def dim(self):
         r"""The global number of degrees of freedom for this function space.
diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 19768a92f3..0e859de5ec 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -787,6 +787,7 @@ def _entity_numbering(self, label):
         for old_component_num, new_component_num in enumerate(renumbering):
             old_pt = old_component_num + self.points._component_offsets[component_index]
             section.setOffset(old_pt, new_component_num)
+            section.setDof(old_pt, 1)
 
         return section
 
diff --git a/firedrake/mg/interface.py b/firedrake/mg/interface.py
index 41b0185bb3..a2bd89168a 100644
--- a/firedrake/mg/interface.py
+++ b/firedrake/mg/interface.py
@@ -1,4 +1,4 @@
-from pyop2 import op2
+import pyop3 as op3
 
 import firedrake
 from firedrake import ufl_expr
@@ -81,13 +81,18 @@ def prolong(coarse, fine):
         # Have to do this, because the node set core size is not right for
         # this expanded stencil
         for d in [coarse, coarse_coords]:
-            d.dat.global_to_local_begin(op2.READ)
-            d.dat.global_to_local_end(op2.READ)
-        op2.par_loop(kernel, next.node_set,
-                     next.dat(op2.WRITE),
-                     coarse.dat(op2.READ, fine_to_coarse),
-                     node_locations.dat(op2.READ),
-                     coarse_coords.dat(op2.READ, fine_to_coarse_coords))
+            d.dat.assemble()
+
+        id_map = utils.owned_node_map(Vf)
+        op3.do_loop(
+            n := Vf.nodes.owned.index(),
+            kernel(
+                next.nodal_dat()[id_map(n)],
+                coarse.nodal_dat()[fine_to_coarse(n)],
+                node_locations.nodal_dat()[id_map(n)],
+                coarse_coords.nodal_dat()[fine_to_coarse_coords(n)],
+            ),
+        )
         coarse = next
         Vc = Vf
     return fine
@@ -125,7 +130,7 @@ def restrict(fine_dual, coarse_dual):
     for j in range(repeat):
         next_level -= 1
         if j == repeat - 1:
-            coarse_dual.dat.eager_zero()
+            coarse_dual.dat.zero()
             next = coarse_dual
         else:
             Vc = firedrake.FunctionSpace(meshes[next_level], element)
@@ -142,14 +147,19 @@ def restrict(fine_dual, coarse_dual):
         # Have to do this, because the node set core size is not right for
         # this expanded stencil
         for d in [coarse_coords]:
-            d.dat.global_to_local_begin(op2.READ)
-            d.dat.global_to_local_end(op2.READ)
+            d.dat.assemble()
         kernel = kernels.restrict_kernel(Vf, Vc)
-        op2.par_loop(kernel, fine_dual.node_set,
-                     next.dat(op2.INC, fine_to_coarse),
-                     fine_dual.dat(op2.READ),
-                     node_locations.dat(op2.READ),
-                     coarse_coords.dat(op2.READ, fine_to_coarse_coords))
+
+        id_map = utils.owned_node_map(Vf)
+        op3.do_loop(
+            n := Vf.nodes.owned.index(),
+            kernel(
+                next.nodal_dat()[fine_to_coarse(n)],
+                fine_dual.nodal_dat()[id_map(n)],
+                node_locations.nodal_dat()[id_map(n)],
+                coarse_coords.nodal_dat()[fine_to_coarse_coords(n)],
+            ),
+        )
         fine_dual = next
         Vf = Vc
     return coarse_dual
@@ -201,7 +211,7 @@ def inject(fine, coarse):
     for j in range(repeat):
         next_level -= 1
         if j == repeat - 1:
-            coarse.dat.eager_zero()
+            coarse.dat.zero()
             next = coarse
             Vc = next.function_space()
         else:
@@ -217,13 +227,18 @@ def inject(fine, coarse):
             # Have to do this, because the node set core size is not right for
             # this expanded stencil
             for d in [fine, fine_coords]:
-                d.dat.global_to_local_begin(op2.READ)
-                d.dat.global_to_local_end(op2.READ)
-            op2.par_loop(kernel, next.node_set,
-                         next.dat(op2.INC),
-                         node_locations.dat(op2.READ),
-                         fine.dat(op2.READ, coarse_node_to_fine_nodes),
-                         fine_coords.dat(op2.READ, coarse_node_to_fine_coords))
+                d.dat.assemble()
+
+            id_map = utils.owned_node_map(Vc)
+            op3.do_loop(
+                n := Vc.nodes.owned.index(),
+                kernel(
+                    next.nodal_dat()[id_map(n)],
+                    node_locations.nodal_dat()[id_map(n)],
+                    fine.nodal_dat()[coarse_node_to_fine_nodes(n)],
+                    fine_coords.nodal_dat()[coarse_node_to_fine_coords(n)],
+                ),
+            )
         else:
             coarse_coords = Vc.mesh().coordinates
             fine_coords = Vf.mesh().coordinates
@@ -232,13 +247,17 @@ def inject(fine, coarse):
             # Have to do this, because the node set core size is not right for
             # this expanded stencil
             for d in [fine, fine_coords]:
-                d.dat.global_to_local_begin(op2.READ)
-                d.dat.global_to_local_end(op2.READ)
-            op2.par_loop(kernel, Vc.mesh().cell_set,
-                         next.dat(op2.INC, next.cell_node_map()),
-                         fine.dat(op2.READ, coarse_cell_to_fine_nodes),
-                         fine_coords.dat(op2.READ, coarse_cell_to_fine_coords),
-                         coarse_coords.dat(op2.READ, coarse_coords.cell_node_map()))
+                d.dat.assemble()
+            op3.do_loop(
+                c := Vc.mesh().cells.owned.index(),
+                kernel(
+                    next.dat[c],
+                    fine.nodal_dat()[coarse_cell_to_fine_nodes(c)],
+                    fine_coords.nodal_dat()[coarse_cell_to_fine_coords(c)],
+                    coarse_coords.nodal_dat()[coarse_coords.function_space().cell_node_map(c)],
+                ),
+            )
+
         fine = next
         Vf = Vc
     return coarse
diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py
index 2399b6a57e..2498d9ccc0 100644
--- a/firedrake/mg/kernels.py
+++ b/firedrake/mg/kernels.py
@@ -1,7 +1,7 @@
 import numpy
 import string
 from fractions import Fraction
-from pyop2 import op2
+import pyop3 as op3
 from firedrake.utils import IntType, as_cstr, complex_mode, ScalarType
 from firedrake.functionspacedata import entity_dofs_key
 from firedrake.functionspaceimpl import FiredrakeDualSpace
@@ -213,10 +213,10 @@ def prolong_kernel(expression):
     levelf = level + Fraction(1, hierarchy.refinements_per_level)
     cache = hierarchy._shared_data_cache["transfer_kernels"]
     coordinates = extract_unique_domain(expression).coordinates
-    if meshc.cell_set._extruded:
+    if meshc.extruded:
         idx = levelf * hierarchy.refinements_per_level
         assert idx == int(idx)
-        assert hierarchy._meshes[int(idx)].cell_set._extruded
+        assert hierarchy._meshes[int(idx)].extruded
     key = (("prolong",)
            + expression.ufl_element().value_shape
            + entity_dofs_key(expression.function_space().finat_element.entity_dofs())
@@ -225,17 +225,12 @@ def prolong_kernel(expression):
         return cache[key]
     except KeyError:
         mesh = extract_unique_domain(coordinates)
-        eval_code = compile_element(expression)
+        evaluate_code = compile_element(expression)
         to_reference_kernel = to_reference_coordinates(coordinates.ufl_element())
         element = create_element(expression.ufl_element())
         coords_element = create_element(coordinates.ufl_element())
 
-        my_kernel = """#include <petsc.h>
-        %(to_reference)s
-        %(evaluate)s
-        __attribute__((noinline)) /* Clang bug */
-        static void pyop2_kernel_prolong(PetscScalar *R, PetscScalar *f, const PetscScalar *X, const PetscScalar *Xc)
-        {
+        kernel = """
             PetscScalar Xref[%(tdim)d];
             int cell = -1;
             int bestcell = -1;
@@ -276,9 +271,7 @@ def prolong_kernel(expression):
                 R[i] = 0;
             }
             pyop2_kernel_evaluate(R, coarsei, Xref);
-        }
-        """ % {"to_reference": str(to_reference_kernel),
-               "evaluate": eval_code,
+        """ % {
                "spacedim": element.cell.get_spatial_dimension(),
                "ncandidate": hierarchy.fine_to_coarse_cells[levelf].shape[1],
                "Rdim": numpy.prod(element.value_shape),
@@ -286,9 +279,27 @@ def prolong_kernel(expression):
                "celldist_l1_c_expr": celldist_l1_c_expr(element.cell, X="Xref"),
                "Xc_cell_inc": coords_element.space_dimension(),
                "coarse_cell_inc": element.space_dimension(),
-               "tdim": mesh.topological_dimension()}
-
-        return cache.setdefault(key, op2.Kernel(my_kernel, name="pyop2_kernel_prolong"))
+               "tdim": mesh.topological_dimension(),
+        }
+        domain = "{ [i]: 0 < i <= 1 }"
+        preambles = [("XXX", "#include <petsc.h>"),
+                     ("YYY", str(to_reference_kernel)),
+                     ("ZZZ", evaluate_code)]
+        varnames = ["R", "f", "X", "Xc"]
+        intents = [op3.WRITE, op3.READ, op3.READ, op3.READ]
+        loopy_kernel = lp.make_kernel(domain,
+            [lp.CInstruction((), kernel, frozenset(varnames[1:]), varnames[:1])],
+            [lp.GlobalArg(varnames[0], ScalarType, None, is_input=False, is_output=True),
+             lp.GlobalArg(varnames[1], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[2], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[3], ScalarType, None, is_input=True, is_output=False),
+            ],
+            name="pyop2_kernel_prolong",
+            preambles=preambles,
+            target=tsfc.parameters.target,
+            lang_version=(2018, 2),
+            )
+        return cache.setdefault(key, op3.Function(loopy_kernel, intents))
 
 
 def restrict_kernel(Vf, Vc):
@@ -313,13 +324,7 @@ def restrict_kernel(Vf, Vc):
         coords_element = create_element(coordinates.ufl_element())
         element = create_element(Vc.ufl_element())
 
-        my_kernel = """#include <petsc.h>
-        %(to_reference)s
-        %(evaluate)s
-
-        __attribute__((noinline)) /* Clang bug */
-        static void pyop2_kernel_restrict(PetscScalar *R, PetscScalar *b, const PetscScalar *X, const PetscScalar *Xc)
-        {
+        kernel = """
             PetscScalar Xref[%(tdim)d];
             int cell = -1;
             int bestcell = -1;
@@ -361,18 +366,35 @@ def restrict_kernel(Vf, Vc):
             const PetscScalar *Ri = R + cell*%(coarse_cell_inc)d;
             pyop2_kernel_evaluate(Ri, b, Xref);
             }
-        }
-        """ % {"to_reference": str(to_reference_kernel),
-               "evaluate": evaluate_code,
+        """ % {
                "ncandidate": hierarchy.fine_to_coarse_cells[levelf].shape[1],
                "inside_cell": inside_check(element.cell, eps=1e-8, X="Xref"),
                "celldist_l1_c_expr": celldist_l1_c_expr(element.cell, X="Xref"),
                "Xc_cell_inc": coords_element.space_dimension(),
                "coarse_cell_inc": element.space_dimension(),
                "spacedim": element.cell.get_spatial_dimension(),
-               "tdim": mesh.topological_dimension()}
-
-        return cache.setdefault(key, op2.Kernel(my_kernel, name="pyop2_kernel_restrict"))
+               "tdim": mesh.topological_dimension(),
+        }
+        domain = "{ [i]: 0 < i <= 1 }"
+        preambles = [("XXX", "#include <petsc.h>"),
+                     ("YYY", str(to_reference_kernel)),
+                     ("ZZZ", evaluate_code)]
+        varnames = ["R", "b", "X", "Xc"]
+        intents = [op3.INC, op3.READ, op3.READ, op3.READ]
+        loopy_kernel = lp.make_kernel(
+            domain,
+            [lp.CInstruction((), kernel, frozenset(varnames[1:]), varnames[:1])],
+            [lp.GlobalArg(varnames[0], ScalarType, None, is_input=True, is_output=True),
+             lp.GlobalArg(varnames[1], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[2], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[3], ScalarType, None, is_input=True, is_output=False),
+            ],
+            name="pyop2_kernel_restrict",
+            preambles=preambles,
+            target=tsfc.parameters.target,
+            lang_version=(2018, 2),
+            )
+        return cache.setdefault(key, op3.Function(loopy_kernel, intents))
 
 
 def inject_kernel(Vf, Vc):
@@ -404,12 +426,6 @@ def inject_kernel(Vf, Vc):
         coords_element = create_element(coordinates.ufl_element())
         Vf_element = create_element(Vf.ufl_element())
         kernel = """
-        %(to_reference)s
-        %(evaluate)s
-
-        __attribute__((noinline)) /* Clang bug */
-        static void pyop2_kernel_inject(PetscScalar *R, const PetscScalar *X, const PetscScalar *f, const PetscScalar *Xf)
-        {
             PetscScalar Xref[%(tdim)d];
             int cell = -1;
             int bestcell = -1;
@@ -449,10 +465,7 @@ def inject_kernel(Vf, Vc):
                 R[i] = 0;
             }
             pyop2_kernel_evaluate(R, fi, Xref);
-        }
         """ % {
-            "to_reference": str(to_reference_kernel),
-            "evaluate": evaluate_code,
             "inside_cell": inside_check(Vc.finat_element.cell, eps=1e-8, X="Xref"),
             "spacedim": Vc.finat_element.cell.get_spatial_dimension(),
             "celldist_l1_c_expr": celldist_l1_c_expr(Vc.finat_element.cell, X="Xref"),
@@ -460,9 +473,28 @@ def inject_kernel(Vf, Vc):
             "ncandidate": ncandidate,
             "Rdim": numpy.prod(Vf_element.value_shape),
             "Xf_cell_inc": coords_element.space_dimension(),
-            "f_cell_inc": Vf_element.space_dimension()
+            "f_cell_inc": Vf_element.space_dimension(),
         }
-        return cache.setdefault(key, (op2.Kernel(kernel, name="pyop2_kernel_inject"), False))
+        domain = "{ [i]: 0 < i <= 1 }"
+        preambles = [("XXX", "#include <petsc.h>"),
+                     ("YYY", str(to_reference_kernel)),
+                     ("ZZZ", evaluate_code)]
+        varnames = ["R", "X", "f", "Xf"]
+        intents = [op3.INC, op3.READ, op3.READ, op3.READ]
+        loopy_kernel = lp.make_kernel(
+            domain,
+            [lp.CInstruction((), kernel, frozenset(varnames[1:]), varnames[:1])],
+            [lp.GlobalArg(varnames[0], ScalarType, None, is_input=True, is_output=True),
+             lp.GlobalArg(varnames[1], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[2], ScalarType, None, is_input=True, is_output=False),
+             lp.GlobalArg(varnames[3], ScalarType, None, is_input=True, is_output=False),
+            ],
+            name="pyop2_kernel_inject",
+            preambles=preambles,
+            target=tsfc.parameters.target,
+            lang_version=(2018, 2),
+            )
+        return cache.setdefault(key, (op3.Function(loopy_kernel, intents), False))
 
 
 class MacroKernelBuilder(firedrake_interface.KernelBuilderBase):
@@ -707,9 +739,12 @@ def name_multiindex(multiindex, name):
         domains, instructions, kernel_data, name=kernel_name,
         target=tsfc.parameters.target, lang_version=(2018, 2))
     kernel = lp.merge([kernel, *subkernels])
-    return op2.Kernel(
-        kernel, name=kernel_name, include_dirs=Ainv.include_dirs,
-        headers=Ainv.headers, events=Ainv.events)
+    kernel = kernel.with_entrypoints({kernel_name})
+    intents = [op3.INC, op3.READ, op3.READ, op3.READ]
+    return op3.Function(kernel, intents)
+    #return op2.Kernel(
+        #kernel, name=kernel_name, include_dirs=Ainv.include_dirs,
+        #headers=Ainv.headers, events=Ainv.events)
 
 
 def _generate_call_insn(name, args, *, iname_prefix=None, **kwargs):
diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py
index a40fad62f7..d686927a19 100644
--- a/firedrake/mg/mesh.py
+++ b/firedrake/mg/mesh.py
@@ -230,7 +230,7 @@ def ExtrudedMeshHierarchy(base_hierarchy, height, base_layer=-1, refinement_rati
     """
     if not isinstance(base_hierarchy, HierarchyBase):
         raise ValueError("Expecting a HierarchyBase, not a %r" % type(base_hierarchy))
-    if any(m.cell_set._extruded for m in base_hierarchy):
+    if any(m.extruded for m in base_hierarchy):
         raise ValueError("Meshes in base hierarchy must not be extruded")
 
     if layers is None:
@@ -287,7 +287,7 @@ def SemiCoarsenedExtrudedHierarchy(base_mesh, height, nref=1, base_layer=-1, ref
     if not isinstance(base_mesh, firedrake.mesh.MeshGeometry):
         raise ValueError(f"Can only extruded a mesh, not a {type(base_mesh)}")
     base_mesh.init()
-    if base_mesh.cell_set._extruded:
+    if base_mesh.extruded:
         raise ValueError("Base mesh must not be extruded")
     if layers is None:
         if base_layer == -1:
@@ -306,7 +306,7 @@ def SemiCoarsenedExtrudedHierarchy(base_mesh, height, nref=1, base_layer=-1, ref
                            gdim=gdim)
               for layer in layers]
     refinements_per_level = 1
-    identity = np.arange(base_mesh.cell_set.size, dtype=IntType).reshape(-1, 1)
+    identity = np.arange(base_mesh.owned_cells.size, dtype=IntType).reshape(-1, 1)
     coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), identity)
                                 for i in range(nref))
     fine_to_coarse_cells = dict((Fraction(i+1, refinements_per_level), identity)
diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py
index f3eb475b50..16fdc6c62c 100644
--- a/firedrake/mg/utils.py
+++ b/firedrake/mg/utils.py
@@ -1,17 +1,52 @@
 import numpy
 from fractions import Fraction
-from pyop2 import op2
+import pyop3 as op3
 from firedrake.utils import IntType
 from firedrake.functionspacedata import entity_dofs_key
 import finat.ufl
 import firedrake
 from firedrake.cython import mgimpl as impl
+from pyrsistent import freeze
+
+
+def coarse_to_fine_cell_map(coarse_mesh, fine_mesh, coarse_to_fine_data):
+    connectivity = {
+        freeze({coarse_mesh.name: coarse_mesh.cell_label}): [
+            op3.TabulatedMapComponent(fine_mesh.name, fine_mesh.cell_label, coarse_to_fine_data)
+        ]
+    }
+    return op3.Map(connectivity)
+
+
+def create_node_map(iterset, toset, arity=1, values=None):
+    axes = op3.AxisTree.from_iterable([iterset, arity])
+    if values is None:
+        values = numpy.arange(iterset.size, dtype=IntType)
+    dat = op3.HierarchicalArray(axes, data=values)
+    return op3.Map({
+        freeze({iterset.label: iterset.owned.component.label}): [
+            op3.TabulatedMapComponent(toset.label, toset.component.label, dat)
+        ]
+    })
+
+
+def owned_node_map(V):
+    """ This should not be necessary
+    """
+    mesh = V.mesh()
+    key = entity_dofs_key(V.finat_element.entity_dofs())
+    cache = mesh._shared_data_cache["owned_node_map"]
+    try:
+        return cache[key]
+    except KeyError:
+        node_map = create_node_map(V.nodes.owned, V.nodes)
+        return cache.setdefault(key, node_map)
 
 
 def fine_node_to_coarse_node_map(Vf, Vc):
     if len(Vf) > 1:
         assert len(Vf) == len(Vc)
-        return op2.MixedMap(fine_node_to_coarse_node_map(f, c) for f, c in zip(Vf, Vc))
+        return op3.MixedMap(fine_node_to_coarse_node_map(f, c) for f, c in zip(Vf, Vc))
     mesh = Vf.mesh()
     assert hasattr(mesh, "_shared_data_cache")
     hierarchyf, levelf = get_level(Vf.mesh())
@@ -41,15 +76,15 @@ def fine_node_to_coarse_node_map(Vf, Vc):
 
         fine_to_coarse = hierarchy.fine_to_coarse_cells[levelf]
         fine_to_coarse_nodes = impl.fine_to_coarse_nodes(Vf, Vc, fine_to_coarse)
-        return cache.setdefault(key, op2.Map(Vf.node_set, Vc.node_set,
-                                             fine_to_coarse_nodes.shape[1],
-                                             values=fine_to_coarse_nodes))
+        return cache.setdefault(key, create_node_map(Vf.nodes, Vc.nodes,
+                                                     arity=fine_to_coarse_nodes.shape[1],
+                                                     values=fine_to_coarse_nodes))
 
 
 def coarse_node_to_fine_node_map(Vc, Vf):
     if len(Vf) > 1:
         assert len(Vf) == len(Vc)
-        return op2.MixedMap(coarse_node_to_fine_node_map(f, c) for f, c in zip(Vf, Vc))
+        return op3.MixedMap(coarse_node_to_fine_node_map(f, c) for f, c in zip(Vf, Vc))
     mesh = Vc.mesh()
     assert hasattr(mesh, "_shared_data_cache")
     hierarchyf, levelf = get_level(Vf.mesh())
@@ -79,15 +114,15 @@ def coarse_node_to_fine_node_map(Vc, Vf):
 
         coarse_to_fine = hierarchy.coarse_to_fine_cells[levelc]
         coarse_to_fine_nodes = impl.coarse_to_fine_nodes(Vc, Vf, coarse_to_fine)
-        return cache.setdefault(key, op2.Map(Vc.node_set, Vf.node_set,
-                                             coarse_to_fine_nodes.shape[1],
-                                             values=coarse_to_fine_nodes))
+        return cache.setdefault(key, create_node_map(Vc.nodes, Vf.nodes,
+                                                     arity=coarse_to_fine_nodes.shape[1],
+                                                     values=coarse_to_fine_nodes))
 
 
 def coarse_cell_to_fine_node_map(Vc, Vf):
     if len(Vf) > 1:
         assert len(Vf) == len(Vc)
-        return op2.MixedMap(coarse_cell_to_fine_node_map(f, c) for f, c in zip(Vf, Vc))
+        return op3.MixedMap(coarse_cell_to_fine_node_map(f, c) for f, c in zip(Vf, Vc))
     mesh = Vc.mesh()
     assert hasattr(mesh, "_shared_data_cache")
     hierarchyf, levelf = get_level(Vf.mesh())
@@ -115,22 +150,27 @@ def coarse_cell_to_fine_node_map(Vc, Vf):
             level_ratio = 1
         coarse_to_fine = hierarchy.coarse_to_fine_cells[levelc]
         _, ncell = coarse_to_fine.shape
-        iterset = Vc.mesh().cell_set
+        iterset = Vc.mesh().owned_cells
         arity = Vf.finat_element.space_dimension() * ncell
-        coarse_to_fine_nodes = numpy.full((iterset.total_size, arity*level_ratio), -1, dtype=IntType)
-        values = Vf.cell_node_map().values[coarse_to_fine, :].reshape(iterset.size, arity)
+        coarse_to_fine_nodes = numpy.full((Vc.mesh().cells.size, arity*level_ratio), -1, dtype=IntType)
+        values = Vf.owned_cell_node_list[coarse_to_fine, :].reshape(iterset.size, arity)
 
         if Vc.extruded:
             off = numpy.tile(Vf.offset, ncell)
-            coarse_to_fine_nodes[:Vc.mesh().cell_set.size, :] = numpy.hstack([values + off*i for i in range(level_ratio)])
+            coarse_to_fine_nodes[:iterset.size, :] = numpy.hstack([values + off*i for i in range(level_ratio)])
         else:
-            coarse_to_fine_nodes[:Vc.mesh().cell_set.size, :] = values
-        offset = Vf.offset
-        if offset is not None:
-            offset = numpy.tile(offset*level_ratio, ncell*level_ratio)
-        return cache.setdefault(key, op2.Map(iterset, Vf.node_set,
-                                             arity=arity*level_ratio, values=coarse_to_fine_nodes,
-                                             offset=offset))
+            coarse_to_fine_nodes[:iterset.size, :] = values
+        if Vf.extruded:
+            offset = numpy.tile(Vf.offset*level_ratio, ncell*level_ratio)
+
+        axes = op3.AxisTree.from_iterable([Vc.nodes, arity*level_ratio])
+        coarse_cell_to_fine_node_dat = op3.HierarchicalArray(axes, data=coarse_to_fine_nodes)
+        coarse_cell_to_fine_node_map = op3.Map({
+            freeze({Vc.mesh().topology.name: Vc.mesh().cell_label}): [
+                op3.TabulatedMapComponent(Vf.nodes.label, Vf.nodes.component.label, coarse_cell_to_fine_node_dat)
+            ]
+        })
+        return cache.setdefault(key, coarse_cell_to_fine_node_map)
 
 
 def physical_node_locations(V):
diff --git a/firedrake/output/vtk_output.py b/firedrake/output/vtk_output.py
index b6f7548216..e6550e6c62 100644
--- a/firedrake/output/vtk_output.py
+++ b/firedrake/output/vtk_output.py
@@ -190,7 +190,7 @@ def get_topology(coordinates):
 
     # Repeat up the column
     num_cells = mesh.cell_set.size
-    if not mesh.cell_set._extruded:
+    if not mesh.extruded:
         cell_layers = 1
         offsets = 0
     else:
diff --git a/firedrake/parloops.py b/firedrake/parloops.py
index b5ab7540d9..e5fbf61c88 100644
--- a/firedrake/parloops.py
+++ b/firedrake/parloops.py
@@ -368,6 +368,60 @@ def pack_pyop3_tensor(tensor: Any, *args, **kwargs):
     raise TypeError(f"No handler defined for {type(tensor).__name__}")
 
 
+@pack_pyop3_tensor.register
+def _(
+    axes: op3.AxisTree,
+    V: WithGeometry,
+    index: op3.LoopIndex,
+    integral_type: str,
+):
+    plex = V.mesh().topology
+
+    if integral_type == "cell":
+        # TODO ideally the FIAT permutation would not need to be known
+        # about by the mesh topology and instead be handled here. This
+        # would probably require an oriented mesh
+        # (see https://github.com/firedrakeproject/firedrake/pull/3332)
+        # indexed = array.getitem(plex._fiat_closure(index), strict=True)
+        pack_indices = _cell_integral_pack_indices(V, index)
+    elif integral_type in {"exterior_facet", "interior_facet"}:
+        pack_indices = _facet_integral_pack_indices(V, index)
+    else:
+        raise NotImplementedError
+
+    indexed = axes.getitem(pack_indices, strict=True)
+
+    if plex.ufl_cell().is_simplex():
+        return indexed
+
+    if plex.ufl_cell() == ufl.hexahedron:
+        perms = _entity_permutations(V)
+        mytree = _orientations(plex, perms, index, integral_type)
+        mytree = _with_shape_indices(V, mytree, integral_type in {"exterior_facet", "interior_facet"})
+
+        indexed = indexed.getitem(mytree, strict=True)
+
+    tensor_axes, ttarget_paths, tindex_exprs = _tensorify_axes(V.finat_element.product)
+    tensor_axes, ttarget_paths, tindex_exprs = _with_shape_axes(V, tensor_axes, ttarget_paths, tindex_exprs, integral_type)
+
+    # This should be cleaned up - basically we need to accumulate the target_paths
+    # and index_exprs along the nodes. This is done inside index_axes in pyop3.
+    from pyop3.itree.tree import _acc_target_paths
+    ttarget_paths = _acc_target_paths(tensor_axes, ttarget_paths)
+    tindex_exprs = _acc_target_paths(tensor_axes, tindex_exprs)
+
+    tensor_axes = op3.IndexedAxisTree(
+        tensor_axes.node_map,
+        indexed.axes.unindexed,
+        target_paths=ttarget_paths,
+        index_exprs=tindex_exprs,
+        layout_exprs={},
+        outer_loops=indexed.axes.outer_loops,
+    )
+    composed_axes = compose_axes(tensor_axes, indexed.axes)
+    return composed_axes
+
+
 @pack_pyop3_tensor.register
 def _(
     array: op3.HierarchicalArray,
diff --git a/firedrake/preconditioners/asm.py b/firedrake/preconditioners/asm.py
index 72caabb5c6..b4d735bc08 100644
--- a/firedrake/preconditioners/asm.py
+++ b/firedrake/preconditioners/asm.py
@@ -143,7 +143,7 @@ class ASMStarPC(ASMPatchPC):
     def get_patches(self, V):
         mesh = V._mesh
         mesh_dm = mesh.topology_dm
-        if mesh.cell_set._extruded:
+        if mesh.extruded:
             warning("applying ASMStarPC on an extruded mesh")
 
         # Obtain the topological entities to use to construct the stars
@@ -286,7 +286,7 @@ class ASMLinesmoothPC(ASMPatchPC):
 
     def get_patches(self, V):
         mesh = V._mesh
-        assert mesh.cell_set._extruded
+        assert mesh.extruded
         dm = mesh.topology_dm
         section = V.dm.getDefaultSection()
         # Obtain the codimensions to loop over from options, if present
@@ -387,7 +387,7 @@ def get_patches(self, V):
         mesh = V.mesh()
         mesh_dm = mesh.topology_dm
         nlayers = mesh.layers
-        if not mesh.cell_set._extruded:
+        if not mesh.extruded:
             return super(ASMExtrudedStarPC, self).get_patches(V)
 
         # Obtain the topological entities to use to construct the stars
diff --git a/firedrake/preconditioners/fdm.py b/firedrake/preconditioners/fdm.py
index 893920d49b..b6425a6af4 100644
--- a/firedrake/preconditioners/fdm.py
+++ b/firedrake/preconditioners/fdm.py
@@ -2138,7 +2138,7 @@ def assemble_coefficients(self, J, fcp):
         if Piola:
             # make DGT functions with the second order coefficient
             # and the Piola tensor for each side of each facet
-            extruded = mesh.cell_set._extruded
+            extruded = mesh.extruded
             dS_int = ufl.dS_h(degree=quad_deg) + ufl.dS_v(degree=quad_deg) if extruded else ufl.dS(degree=quad_deg)
             area = ufl.FacetArea(mesh)
             ifacet_inner = lambda v, u: ((ufl.inner(v('+'), u('+')) + ufl.inner(v('-'), u('-')))/area)*dS_int
@@ -2306,7 +2306,7 @@ def extrude_interior_facet_maps(V):
     facet_to_nodes = facet_node_map.values
     nbase = facet_to_nodes.shape[0]
 
-    if mesh.cell_set._extruded:
+    if mesh.extruded:
         facet_offset = facet_node_map.offset
         local_facet_data_h = numpy.array([5, 4], local_facet_data.dtype)
 
diff --git a/firedrake/preconditioners/patch.py b/firedrake/preconditioners/patch.py
index e73c41b129..e9f883db34 100644
--- a/firedrake/preconditioners/patch.py
+++ b/firedrake/preconditioners/patch.py
@@ -793,7 +793,7 @@ def initialize(self, obj):
         self.ctx = ctx
         self.plex.setAttr("__firedrake_ctx__", weakref.proxy(ctx))
 
-        if mesh.cell_set._extruded:
+        if mesh.extruded:
             raise NotImplementedError("Not implemented on extruded meshes")
 
         if "overlap_type" not in mesh._distribution_parameters:
diff --git a/firedrake/slate/slac/kernel_builder.py b/firedrake/slate/slac/kernel_builder.py
index edf1fabfda..4587f56327 100644
--- a/firedrake/slate/slac/kernel_builder.py
+++ b/firedrake/slate/slac/kernel_builder.py
@@ -191,7 +191,7 @@ def layer_integral_predicates(self, tensor, integral_type):
     def facet_integral_predicates(self, mesh, integral_type, kinfo, subdomain_id):
         self.bag.needs_cell_facets = True
         # Number of recerence cell facets
-        if mesh.cell_set._extruded:
+        if mesh.extruded:
             self.num_facets = mesh._base_mesh.ufl_cell().num_facets()
         else:
             self.num_facets = mesh.ufl_cell().num_facets()
diff --git a/firedrake/slate/static_condensation/hybridization.py b/firedrake/slate/static_condensation/hybridization.py
index 5dd592c44e..5b1c501dcf 100644
--- a/firedrake/slate/static_condensation/hybridization.py
+++ b/firedrake/slate/static_condensation/hybridization.py
@@ -128,7 +128,7 @@ def initialize(self, pc):
         n = ufl.FacetNormal(mesh)
         sigma = TrialFunctions(V_d)[self.vidx]
 
-        if mesh.cell_set._extruded:
+        if mesh.extruded:
             Kform = (gammar('+') * ufl.jump(sigma, n=n) * ufl.dS_h
                      + gammar('+') * ufl.jump(sigma, n=n) * ufl.dS_v)
         else:
@@ -166,7 +166,7 @@ def initialize(self, pc):
             integrand = gammar * ufl.dot(sigma, n)
             measures = []
             trace_subdomains = []
-            if mesh.cell_set._extruded:
+            if mesh.extruded:
                 ds = ufl.ds_v
                 for subdomain in sorted(extruded_neumann_subdomains):
                     measures.append({"top": ufl.ds_t, "bottom": ufl.ds_b}[subdomain])
@@ -192,7 +192,7 @@ def initialize(self, pc):
             # the exterior boundary. Extruded cells will have both
             # horizontal and vertical facets
             trace_subdomains = ["on_boundary"]
-            if mesh.cell_set._extruded:
+            if mesh.extruded:
                 trace_subdomains.extend(["bottom", "top"])
             trace_bcs = [DirichletBC(TraceSpace, Constant(0.0), subdomain) for subdomain in trace_subdomains]
 
diff --git a/tests/conftest.py b/tests/conftest.py
index 232bb46925..5e33c8a95a 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,7 +4,7 @@
 
 
 # TODO pyop3
-collect_ignore_glob = ["vertexonly/*", "extrusion/*", "multigrid/*", "supermesh/*", "external_operators/*"]
+collect_ignore_glob = ["vertexonly/*", "extrusion/*", "supermesh/*", "external_operators/*"]
 
 
 def pytest_configure(config):
diff --git a/tests/multigrid/test_grid_transfer.py b/tests/multigrid/test_grid_transfer.py
index 40b21bb3b2..02134322d2 100644
--- a/tests/multigrid/test_grid_transfer.py
+++ b/tests/multigrid/test_grid_transfer.py
@@ -21,16 +21,15 @@ def space(request, cell):
         return request.param
 
 
-@pytest.fixture(params=[1, 2], scope="module")
-def refinements_per_level(request):
+@pytest.fixture(params=[(1, 1), (2, 1), (1, 2)], ids=["1-1", "2-1", "1-2"], scope="module")
+def nref_refinements_per_level(request):
     return request.param
 
 
 @pytest.fixture(scope="module")
-def hierarchy(cell, refinements_per_level):
+def hierarchy(cell, nref_refinements_per_level):
     if cell == "interval":
         mesh = UnitIntervalMesh(3)
-        return MeshHierarchy(mesh, 2)
     elif cell in {"triangle", "triangle-nonnested", "prism"}:
         mesh = UnitSquareMesh(3, 3, quadrilateral=False)
     elif cell in {"quadrilateral", "hexahedron"}:
@@ -38,7 +37,7 @@ def hierarchy(cell, refinements_per_level):
     elif cell == "tetrahedron":
         mesh = UnitCubeMesh(2, 2, 2)
 
-    nref = {2: 1, 1: 2}[refinements_per_level]
+    nref, refinements_per_level = nref_refinements_per_level
     hierarchy = MeshHierarchy(mesh, nref, refinements_per_level=refinements_per_level)
 
     if cell in {"prism", "hexahedron"}:
@@ -324,3 +323,9 @@ def test_grid_transfer_periodic(periodic_hierarchy, periodic_space):
         run_injection(periodic_hierarchy, vector, periodic_space, degrees, exact=exact_primal_periodic)
     run_prolongation(periodic_hierarchy, vector, periodic_space, degrees, exact=exact_primal_periodic)
     run_restriction(periodic_hierarchy, vector, periodic_space, degrees)
+
+
+if __name__ == "__main__":
+    bmesh = UnitIntervalMesh(3)
+    mh = MeshHierarchy(bmesh, 2, refinements_per_level=1)
+    run_prolongation(mh, False, "CG", [1])
diff --git a/tests/regression/test_fdm.py b/tests/regression/test_fdm.py
index d476081e8f..0f39b3d3fa 100644
--- a/tests/regression/test_fdm.py
+++ b/tests/regression/test_fdm.py
@@ -85,7 +85,7 @@
 def build_riesz_map(V, d):
     beta = Constant(1E-4)
     subs = [(1, 3)]
-    if V.mesh().cell_set._extruded:
+    if V.mesh().extruded:
         subs += ["top"]
 
     x = SpatialCoordinate(V.mesh())
@@ -191,7 +191,7 @@ def test_variable_coefficient(mesh):
     L = inner(v, Constant(1))*dx
 
     subs = ("on_boundary",)
-    if mesh.cell_set._extruded:
+    if mesh.extruded:
         subs += ("top", "bottom")
     bcs = [DirichletBC(V, zero(V.ufl_element().value_shape), sub) for sub in subs]
 
@@ -250,7 +250,7 @@ def test_ipdg_direct_solver(fs):
     alpha = lambda grad_u: dot(dot(A2, grad_u), A1)
     beta = diag(Constant(range(2, ncomp+2)))
 
-    extruded = mesh.cell_set._extruded
+    extruded = mesh.extruded
     subs = (1,)
     if gdim > 1:
         subs += (3,)