diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 311de9d4..5b3a080b 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -180,6 +180,17 @@ jobs:
         run: |
           . venv/bin/activate
           python3 -m pip install --break-system-packages -e './fiat-repo'
+    
+      - uses: actions/checkout@v5
+        with:
+          repository: firedrakeproject/ufl
+          path: ufl-repo
+          ref: refs/heads/indiamai/hash_directional_sobolev
+        
+      - name: Install checked out UFL
+        run: |
+          . venv/bin/activate
+          python3 -m pip install --break-system-packages -e './ufl-repo'
 
       - name: Run tests
         run: |
diff --git a/fuse/cells.py b/fuse/cells.py
index b8a5e903..c3e51e70 100644
--- a/fuse/cells.py
+++ b/fuse/cells.py
@@ -10,11 +10,13 @@
 from matplotlib.patches import FancyArrowPatch
 from mpl_toolkits.mplot3d import proj3d
 from sympy.combinatorics.named_groups import SymmetricGroup
-from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple, orientation_value
+from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple, orientation_value, _SYMBOLS, as_tuple
 from FIAT.reference_element import Simplex, TensorProductCell as FiatTensorProductCell, Hypercube
 from FIAT.quadrature_schemes import create_quadrature
 from ufl.cell import Cell, TensorProductCell
 from functools import cache
+from itertools import product
+from collections import defaultdict
 
 
 class Arrow3D(FancyArrowPatch):
@@ -855,6 +857,34 @@ def quadrature(self, degree):
         pts, wts = Q.get_points(), Q.get_weights()
         return pts, wts
 
+    def volume(self):
+        vertices = np.asarray(self.ordered_vertex_coords())
+        if self.get_spatial_dimension() == 0:
+            return 1
+        elif self.get_spatial_dimension() == 1:
+            return abs(vertices[1] - vertices[0])[0]
+        elif self.get_spatial_dimension() == 2:
+            x = vertices[:, 0]
+            y = vertices[:, 1]
+            return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+        elif self.get_spatial_dimension() == 3:
+            vertices = np.asarray(vertices)
+            V = 0.0
+            for face in self.d_entities(2):
+                pts = np.array([self.get_node(v, return_coords=True) for v in face.ordered_vertices()])
+                c = pts.mean(axis=0)
+
+                n = np.zeros(3)
+                for i in range(len(pts)):
+                    v0 = pts[i]
+                    v1 = pts[(i + 1) % len(pts)]
+                    n += np.cross(v0, v1)
+                V += np.dot(c, n)
+            V = abs(V) / 3.0
+            return V
+        else:
+            raise NotImplementedError("Dimension not accounted for")
+
     def cartesian_to_barycentric(self, pts):
         verts = np.array(self.ordered_vertex_coords())
         v_0 = self.ordered_vertex_coords()[0]
@@ -987,41 +1017,49 @@ def _from_dict(o_dict):
 class TensorProductPoint():
     id_iter = itertools.count()
 
-    def __init__(self, A, B):
+    def __init__(self, *factors):
         self.id = next(self.id_iter)
-        self.A = A
-        self.B = B
+        self.A = factors[0]
+        self.B = factors[1]
+        self.factors = factors
         self.dimension = self.A.dimension + self.B.dimension
         self.flat = False
         self.fiat_elem = None
         self.group = self.compute_cell_group()
         self.entities = {}
-        for d in [(a_d, b_d) for a_d in range(self.A.dimension + 1) for b_d in range(self.B.dimension + 1)]:
-            if d == (self.A.dimension, self.B.dimension):
+
+        degree_tuples = list(product(*(range(f.dimension + 1) for f in factors)))
+        for d in degree_tuples:
+            if d == tuple(f.dimension for f in factors):
                 self.entities[d] = [self]
             else:
-                self.entities[d] = [TensorProductPoint(e_a, e_b) for e_a in self.A.d_entities(d[0], True) for e_b in self.B.d_entities(d[1], True)]
+                self.entities[d] = [TensorProductPoint(*entities) for entities in product(*(f.d_entities(degree, True) for f, degree in zip(factors, d)))]
 
     def ordered_vertices(self):
-        return self.A.ordered_vertices() + self.B.ordered_vertices()
+        return self.entities[0]
+
+    def ordered_vertex_coords(self):
+        return [sum(verts, ()) for verts in product(*(f.vertices(return_coords=True) for f in self.factors))]
 
     def component_orientations(self):
         from fuse.utils import orientation_value
         self.component_os_to_os = {}
         for dim in self.to_fiat().get_topology():
             self.component_os_to_os[dim] = {}
-            a_ent = self.A.d_entities(dim[0])[0]
-            b_ent = self.B.d_entities(dim[1])[0]
-            verts = [(v_a, v_b) for v_a in a_ent.vertices() for v_b in b_ent.vertices()]
+            ents = [f.d_entities(d)[0] for f, d in zip(self.factors, dim)]
+            verts = list(product(*(e.vertices() for e in ents)))
             ident = [i for i in range(len(verts))]
-            group = [(g_a, g_b) for g_a in a_ent.group.members() for g_b in b_ent.group.members()]
-            for g_a, g_b in group:
-                new_verts = [(v_a, v_b) for v_a in g_a.permute(a_ent.vertices()) for v_b in g_b.permute(b_ent.vertices())]
+            group = list(product(*(e.group.members() for e in ents)))
+            for gs in group:
+                new_verts = list(product(*(g.permute(f.vertices()) for g, f in zip(gs, ents))))
                 perm = [verts.index(v) for v in new_verts]
                 o_val = orientation_value(ident, perm)
                 if sum(dim) == self.dimension and self.group.group_rep_numbering is not None:
                     o_val = self.group.group_rep_numbering[o_val]
-                self.component_os_to_os[dim][(g_a.numeric_rep(), g_b.numeric_rep())] = o_val
+                # if o_val in [m.numeric_rep() for m in self.group.members()]:
+                self.component_os_to_os[dim][tuple(g.numeric_rep() for g in gs)] = o_val
+                # else:
+                #     breakpoint()
         return self.component_os_to_os
 
     def compute_cell_group(self):
@@ -1029,10 +1067,12 @@ def compute_cell_group(self):
         Systematically work out the symmetry group of the tensor product cell.
         """
         verts = self.vertices()
-        group = [(g_a, g_b) for g_a in self.A.group.members() for g_b in self.B.group.members()]
+        group = list(product(*(f.group.members() for f in self.factors)))
+        # group = [(g_a, g_b) for g_a in self.A.group.members() for g_b in self.B.group.members()]
         perms = []
-        for g_a, g_b in group:
-            new_verts = [(v_a, v_b) for v_a in g_a.permute(self.A.vertices()) for v_b in g_b.permute(self.B.vertices())]
+        for gs in group:
+            new_verts = list(product(*(g.permute(f.vertices()) for g, f in zip(gs, self.factors))))
+            # new_verts = [(v_a, v_b) for v_a in g_a.permute(self.A.vertices()) for v_b in g_b.permute(self.B.vertices())]
             perm = [verts.index(v) for v in new_verts]
             perms += [fuse_groups.Permutation(perm)]
 
@@ -1041,6 +1081,7 @@ def compute_cell_group(self):
 
     def get_starter_ids(self):
         # this doesn't actually make sense - remove when confirmed all changes to eliminate min ids from triple is done
+        raise NotImplementedError
         a_starts = self.A.get_starter_ids()
         b_starts = self.B.get_starter_ids()
         ids = []
@@ -1055,7 +1096,7 @@ def get_sub_entities(self):
         return self.to_fiat().sub_entities
 
     def dim(self):
-        return (self.A.dimension, self.B.dimension)
+        return tuple(f.dimension for f in self.factors)
 
     def d_entities(self, d, get_class=True):
         if isinstance(d, tuple):
@@ -1063,18 +1104,20 @@ def d_entities(self, d, get_class=True):
                 return self.entities[d]
             return [e.id for e in self.entities[d]]
         raise NotImplementedError("not sure this is right")
-        return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class)
+        # return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class)
 
     def vertices(self, get_class=True, return_coords=False):
         # TODO maybe refactor with get_node
         if return_coords:
-            a_verts = self.A.vertices(return_coords=return_coords)
-            b_verts = self.B.vertices(return_coords=return_coords)
-            return [a + b for a in a_verts for b in b_verts]
-        return [(a, b) for a in self.A.vertices() for b in self.B.vertices()]
+            # a_verts = self.A.vertices(return_coords=return_coords)
+            # b_verts = self.B.vertices(return_coords=return_coords)
+            # return [a + b for a in a_verts for b in b_verts]
+            return [sum(verts, ()) for verts in product(*(f.vertices(return_coords=True) for f in self.factors))]
+        # return [(a, b) for a in self.A.vertices() for b in self.B.vertices()]
+        return list(product(*(f.vertices() for f in self.factors)))
 
     def to_ufl(self, name=None):
-        return TensorProductCell(self.A.to_ufl(), self.B.to_ufl())
+        return TensorProductCell(*[f.to_ufl() for f in self.factors])
 
     def to_fiat(self, name=None):
         if self.fiat_elem is None:
@@ -1082,23 +1125,24 @@ def to_fiat(self, name=None):
         return self.fiat_elem
 
     def flatten(self):
-        assert self.A.equivalent(self.B)
-        return FlattenedPoint(self.A, self.B)
+        assert all(self.factors[0].equivalent(f) for f in self.factors)
+        return FlattenedPoint(*self.factors)
 
 
 class FlattenedPoint(Point, TensorProductPoint):
     d_entities_by_total_d = Point.d_entities
 
-    def __init__(self, A, B):
-        self.A = A
-        self.B = B
-        self.dimension = self.A.dimension + self.B.dimension
+    def __init__(self, *factors):
+        self.A = factors[0]
+        self.B = factors[1]
+        self.factors = factors
+        self.dimension = sum(f.dimension for f in factors)
         self.flat = True
         fuse_edges = self.construct_fuse_rep()
         super().__init__(self.dimension, fuse_edges)
 
     def to_ufl(self, name=None):
-        return CellComplexToUFL(self, "quadrilateral")
+        return CellComplexToUFL(self, name=name)
 
     def to_fiat(self, name=None):
         # TODO this should check if it actually is a hypercube
@@ -1112,68 +1156,118 @@ def d_entities(self, d, get_class=True):
             return self.all_subpoints[d]
         return self.d_entities_by_total_d(d, get_class)
 
+    def tensor_attachment_expr(self, axis, factor_edge, parent_mask, child_mask):
+        """
+        Build the tensor-product attachment as a tuple of SymPy expressions.
+
+        Parameters
+        ----------
+        axis:
+            The factor in which the parent cell is being restricted to a facet.
+
+        factor_edge:
+            The Fuse Edge from the factor parent entity to the factor child entity.
+            Its `.attachment` is expected to be a SymPy expression or tuple of
+            SymPy expressions.
+
+        parent_mask:
+            Dimension tuple of the parent product entity.
+
+        child_mask:
+            Dimension tuple of the child product entity.
+
+        Example
+        -------
+        For parent mask (1, 1), child mask (0, 1), axis 0:
+
+            parent coords: (x, y)
+            attachment might be: (0, y) or (1, y)
+
+        For parent mask (1, 1, 1), child mask (1, 0, 1), axis 1:
+
+            parent coords: (x, y, z)
+            attachment might be: (x, 0, z) or (x, 1, z)
+        """
+        child_dim = sum(child_mask)
+        child_syms = _SYMBOLS[:child_dim]
+        result = tuple()
+        child_offset = 0
+
+        for i, (pdim, cdim) in enumerate(zip(parent_mask, child_mask)):
+            if i == axis:
+                local_expr = as_tuple(factor_edge.attachment)
+                # Substitute the child coordinates belonging to this factor.
+                local_child_syms = child_syms[child_offset:child_offset + cdim]
+                local_child_symbols = _SYMBOLS[:cdim]
+                subs = {old: new for old, new in zip(local_child_symbols, local_child_syms)}
+                mapped = tuple(sp.sympify(expr).subs(subs) for expr in local_expr)
+                for comp in mapped:
+                    result += comp
+                child_offset += cdim
+            else:
+                # Identity map on unchanged tensor factors.
+                result += tuple(child_syms[child_offset:child_offset + cdim])
+                child_offset += cdim
+        return result
+
     def construct_fuse_rep(self):
-        sub_cells = [self.A, self.B]
-        dims = (self.A.dimension, self.B.dimension)
-        if sum(dims) > 2:
-            raise NotImplementedError("Flattening 3D tensor products not yet implemented")
-
-        self.all_subpoints = {(a_d, b_d): [] for a_d in range(self.A.dimension + 1) for b_d in range(self.B.dimension + 1)}
-        points = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells}
-        attachments = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells}
-
-        for d in range(max(dims) + 1):
-            for cell in sub_cells:
-                if d <= cell.dimension:
-                    sub_ent = cell.d_entities(d, get_class=True)
-                    points[cell][d].extend(sub_ent)
-                    for s in sub_ent:
-                        attachments[cell][d].extend(s.connections)
-
-        # prod_points = list(itertools.product(*reversed([points[cell][0] for cell in sub_cells])))
-        prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells]))
-        # temp = prod_points[1]
-        # prod_points[1] = prod_points[2]
-        # prod_points[2] = temp
-        point_cls = [Point(0) for i in range(len(prod_points))]
-        self.all_subpoints[(0, 0)] = point_cls
-        edges = []
-
-        # generate edges of tensor product result
-        for a in prod_points:
-            for b in prod_points:
-                # of all combinations of point, take those where at least one changes and at least one is the same
-                if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))):
-                    # ensure if they change, that edge exists in the existing topology
-                    if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i]._topology[1].values()) for i in range(len(sub_cells))]):
-                        edges.append((a, b))
-        # hasse level 1
-        edge_cls1 = {e: None for e in edges}
-        for i in range(len(sub_cells)):
-            for (a, b) in edges:
-                a_idx = prod_points.index(a)
-                b_idx = prod_points.index(b)
-                if a[i] != b[i]:
-                    a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0]
-                    b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0]
-                    edge_cls1[(a, b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o),
-                                                  Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])
-                    edge_tuple_dim = tuple(int((sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i]._topology[1].values())) for i in range(len(sub_cells)))
-                    self.all_subpoints[edge_tuple_dim] += [edge_cls1[(a, b)]]
-        edge_cls2 = []
-        # hasse level 2
-        for i in range(len(sub_cells)):
-            for (a, b) in edges:
-                if a[i] == b[i]:
-                    x = sp.Symbol("x")
-                    a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0]
-                    if i == 0:
-                        attach = (x,) + a_edge.attachment
+        """
+        Construct a Fuse Point for the tensor product of two or three Fuse Point objects.
+        """
+        if len(self.factors) not in (2, 3):
+            raise NotImplementedError("Only 2- and 3-factor tensor products are supported.")
+        top_dim = sum(f.dimension for f in self.factors)
+        # Cache all subentities of each factor by dimension.
+        factor_entities = [{d: tuple(f.d_entities(d, get_class=True)) for d in range(f.dimension + 1)}
+                           for f in self.factors]
+        masks_by_total_dim = defaultdict(list)
+        for mask in product(*(range(f.dimension + 1) for f in self.factors)):
+            masks_by_total_dim[sum(mask)].append(mask)
+
+        product_points = {}
+        all_subpoints = {mask: []
+                         for mask in product(*(range(f.dimension + 1) for f in self.factors))}
+
+        def codim_one_facets(product_entity, mask):
+            """
+            Yield (child_product_entity, axis, factor_edge) for each codim-1 facet.
+            product_entity is a tuple of factor subentities.
+            mask is the corresponding tuple of factor dimensions.
+            """
+            for axis, dim in enumerate(mask):
+                if dim == 0:
+                    continue
+                factor_parent = product_entity[axis]
+                for factor_edge in factor_parent.connections:
+                    child_factor_entity = factor_edge.point
+
+                    child_entity = list(product_entity)
+                    child_entity[axis] = child_factor_entity
+                    child_entity = tuple(child_entity)
+
+                    yield child_entity, axis, factor_edge
+
+        top_level_edges = []
+        for total_dim in range(top_dim + 1):
+            for mask in masks_by_total_dim[total_dim]:
+                for prod_ent in product(*(factor_entities[i][d] for i, d in enumerate(mask))):
+                    if total_dim == 0:
+                        product_point = Point(0)
                     else:
-                        attach = a_edge.attachment + (x,)
-                    edge_cls2.append(Edge(edge_cls1[(a, b)], attach, a_edge.o))
-        self.all_subpoints[(1, 1)] = [self]
-        return edge_cls2
+                        boundary = []
+                        for child_ent, axis, factor_edge in codim_one_facets(prod_ent, mask):
+                            child_point = product_points[child_ent]
+                            child_mask = tuple(e.dimension for e in child_ent)
+                            attach = self.tensor_attachment_expr(axis, factor_edge, mask, child_mask)
+                            boundary.append(Edge(child_point, attach, factor_edge.o))
+                        product_point = Point(total_dim, boundary)
+
+                        if prod_ent == tuple(self.factors):
+                            top_level_edges = boundary
+                    product_points[prod_ent] = product_point
+                    all_subpoints[mask].append(product_point)
+        self.all_subpoints = all_subpoints
+        return top_level_edges
 
     def flatten(self):
         return self
@@ -1204,7 +1298,7 @@ def __init__(self, cell, name=None, renumber=False):
         #     breakpoint()
 
     def cellname(self):
-        return self.name
+        return "FUSE_" + self.name
 
     def construct_subelement(self, dimension, e_id=0, o=None):
         """Constructs the reference element of a cell
@@ -1237,15 +1331,14 @@ def __new__(cls, cell, name=None, *args, **kwargs):
 
     def __init__(self, cell, name=None):
         self.fe_cell = cell
-        self.sub_cells = [cell.A.to_fiat(), cell.B.to_fiat()]
+        fiat_factors = [f.to_fiat() for f in cell.factors]
         if name is None:
-            name = " * ".join([s.name for s in self.sub_cells])
+            name = " * ".join([s.name for s in fiat_factors])
         self.name = name
-# , sub_entities=self.fe_cell.get_sub_entities()
-        super(CellComplexToFiatTensorProduct, self).__init__(cell.A.to_fiat(), cell.B.to_fiat())
+        super(CellComplexToFiatTensorProduct, self).__init__(*fiat_factors)
 
     def cellname(self):
-        return self.name
+        return "FUSE_" + self.name
 
     def construct_subelement(self, dimension):
         """Constructs the reference element of a cell
@@ -1273,10 +1366,11 @@ class CellComplexToFiatHypercube(Hypercube):
 
     def __init__(self, cell, product):
         self.fe_cell = cell
+        self.name = product.name
         super(CellComplexToFiatHypercube, self).__init__(product.get_spatial_dimension(), product)
 
     def cellname(self):
-        return self.name
+        return "FUSE_" + self.name
 
     def construct_subelement(self, dimension):
         """Constructs the reference element of a cell
@@ -1372,10 +1466,11 @@ def constructCellComplex(name):
         # return ufc_tetrahedron().to_ufl(name)
         return make_tetrahedron().to_ufl(name)
     elif name == "hexahedron":
-        import warnings
-        warnings.warn("Hexahedron unimplemented in Fuse")
-        import ufl
-        return ufl.Cell(name)
+        # import warnings
+        # warnings.warn("Hexahedron unimplemented in Fuse")
+        # import ufl
+        # return ufl.Cell(name)
+        return TensorProductPoint(line(), line(), line()).flatten().to_ufl(name)
     elif "*" in name:
         components = [constructCellComplex(c.strip()).cell_complex for c in name.split("*")]
         return TensorProductPoint(*components).to_ufl(name)
diff --git a/fuse/dof.py b/fuse/dof.py
index 376d9c73..3897a843 100644
--- a/fuse/dof.py
+++ b/fuse/dof.py
@@ -3,6 +3,7 @@
 from fuse.traces import TrH1
 import numpy as np
 import sympy as sp
+import numbers
 
 
 class Pairing():
@@ -35,7 +36,7 @@ def __call__(self, kernel, v, cell):
         return v(*kernel.pt)
 
     def tabulate(self):
-        return 1
+        return np.eye(self.entity.dim())
 
     def add_entity(self, entity):
         res = DeltaPairing()
@@ -80,8 +81,8 @@ def tabulate(self):
         if self.orientation:
             new_bvs = np.array(self.entity.orient(self.orientation).basis_vectors())
             basis_change = np.matmul(np.linalg.inv(new_bvs), bvs)
-            return basis_change
-        return np.eye(bvs.shape[0])
+            return (1/self.entity.volume())*basis_change
+        return (1/self.entity.volume())*np.eye(bvs.shape[0])
 
     def add_entity(self, entity):
         res = L2Pairing()
@@ -194,7 +195,7 @@ def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape):
             comps = [[tuple()] for pt in Qpts]
         else:
             comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts]
-        if isinstance(self.pt, tuple) or isinstance(self.pt, int):
+        if isinstance(self.pt, tuple) or isinstance(self.pt, numbers.Number):
             return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps
         if not immersed:
             return Qpts, np.array([wt*np.matmul(self.pt, basis_change) for wt in Qwts]).astype(np.float64), comps
@@ -432,13 +433,7 @@ def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()):
     def to_quadrature(self, arg_degree, value_shape):
         Qpts, Qwts = self.cell_defined_on.quadrature(self.kernel.degree(arg_degree))
         Qwts = Qwts.reshape(Qwts.shape + (1,))
-        dim = self.cell_defined_on.get_spatial_dimension()
-        if dim > 0:
-            bvs = np.array(self.cell_defined_on.basis_vectors())
-            new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors())
-            basis_change = np.matmul(np.linalg.inv(new_bvs), bvs)
-        else:
-            basis_change = np.eye(dim)
+        basis_change = self.pairing.tabulate()
 
         if self.immersed and (isinstance(self.kernel, VectorKernel) or isinstance(self.kernel, BarycentricPolynomialKernel) or isinstance(self.kernel, PolynomialKernel)):
             def immersed(pt):
diff --git a/fuse/enriched.py b/fuse/enriched.py
index 6cedf03d..d417ce0d 100644
--- a/fuse/enriched.py
+++ b/fuse/enriched.py
@@ -10,7 +10,7 @@ class EnrichedElement(ElementTriple):
     In general, FUSE element triples should be represented nodally,
     however this may not be possible for all constructions.
 
-    In particulr, we need to preserve tensor product structure.
+    In particular, we need to preserve tensor product structure.
     """
 
     def __init__(self, A, B, flat=False, symmetric=True, matrices=True):
@@ -87,3 +87,9 @@ def generate(self):
     def to_ufl(self):
         ufl_sub_elements = [e.to_ufl() for e in self.sub_elements]
         return finat.ufl.EnrichedElement(*ufl_sub_elements, triple=self)
+
+    def flatten(self):
+        return EnrichedElement(self.A.flatten(), self.B.flatten(), flat=True, symmetric=self.symmetric, matrices=self.matrices)
+
+    def unflatten(self):
+        return EnrichedElement(self.A.unflatten(), self.B.unflatten(), flat=False, symmetric=self.symmetric, matrices=self.matrices)
diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py
index c7f8d369..8e5fa900 100644
--- a/fuse/spaces/polynomial_spaces.py
+++ b/fuse/spaces/polynomial_spaces.py
@@ -133,6 +133,9 @@ def __hash__(self):
     def restrict(self, mindegree, maxdegree):
         return PolynomialSpace(maxdegree, contains=-1, mindegree=mindegree, set_shape=self.set_shape)
 
+    def to_vector(self):
+        return PolynomialSpace(self.maxdegree, self.contains, self.mindegree, set_shape=True)
+
     def _to_dict(self):
         return {"set_shape": self.set_shape, "min": self.mindegree, "contains": self.contains, "max": self.maxdegree}
 
@@ -231,6 +234,9 @@ def __add__(self, x):
         s.extend([x])
         return ConstructedPolynomialSpace(w, s)
 
+    def to_vector(self):
+        return ConstructedPolynomialSpace(self.weights, [space.to_vector() for space in self.spaces])
+
     def _to_dict(self):
         super_dict = super(ConstructedPolynomialSpace, self)._to_dict()
         super_dict["spaces"] = self.spaces
diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py
index 44275061..c1bed2a6 100644
--- a/fuse/tensor_products.py
+++ b/fuse/tensor_products.py
@@ -4,18 +4,21 @@
 from fuse.enriched import EnrichedElement
 import numpy as np
 from finat.ufl import TensorProductElement, FuseElement, HDivElement, HCurlElement
+from itertools import product
+from functools import reduce
+from collections import defaultdict
 
 
-def tensor_product(A, B, matrices=True):
-    if not (isinstance(A, ElementTriple) and isinstance(B, ElementTriple)):
-        raise ValueError("Both components of Tensor Product need to be a Fuse Triple.")
-    return TensorProductTriple(A, B, matrices=matrices)
+def tensor_product(*factors, matrices=True):
+    if not all(isinstance(f, ElementTriple) for f in factors):
+        raise ValueError("All components of Tensor Product need to be a Fuse Triple.")
+    return TensorProductTriple(*factors, matrices=matrices)
 
 
-def symmetric_tensor_product(A, B, matrices=True):
-    if not (isinstance(A, ElementTriple) and isinstance(B, ElementTriple)):
-        raise ValueError("Both components of Tensor Product need to be a Fuse Triple.")
-    return TensorProductTriple(A, B, matrices=matrices, symmetric=True)
+def symmetric_tensor_product(*factors, matrices=True):
+    if not all(isinstance(f, ElementTriple) for f in factors):
+        raise ValueError("All components of Tensor Product need to be a Fuse Triple.")
+    return TensorProductTriple(*factors, matrices=matrices, symmetric=True)
 
 
 def flatten_dictionary(tensor_dict):
@@ -34,15 +37,14 @@ def flatten_dictionary(tensor_dict):
 
 class TensorProductTriple(ElementTriple):
 
-    def __init__(self, A, B, flat=False, symmetric=True, matrices=True):
-        self.A = A
-        self.B = B
+    def __init__(self, *factors, flat=False, symmetric=True, matrices=True):
+        self.factors = factors
         self.spaces = []
-        for (a, b) in zip(self.A.spaces, self.B.spaces):
-            self.spaces.append(a if a >= b else b)
+        for i in range(len(self.factors[0].spaces)):
+            self.spaces.append(max(f.spaces[i] for f in self.factors))
 
-        self.DOFGenerator = [A.DOFGenerator, B.DOFGenerator]
-        self.cell = TensorProductPoint(A.cell, B.cell)
+        self.DOFGenerator = [f.DOFGenerator for f in self.factors]
+        self.cell = TensorProductPoint(*[f.cell for f in factors])
         self.symmetric = symmetric
         self.flat = flat
         if self.flat:
@@ -59,51 +61,49 @@ def __init__(self, A, B, flat=False, symmetric=True, matrices=True):
 
     @property
     def sub_elements(self):
-        return [self.A, self.B]
+        return self.factors
 
     def __repr__(self):
-        return "TensorProd(%s, %s)" % (repr(self.A), repr(self.B))
+        return f"TensorProd({','.join(['{}' for f in self.factors])})".format(*(repr(f) for f in self.factors))
+
+    def _entity_associations(self, dofs, overall=True):
+        return self.entity_assocs, None, None
 
     def setup_matrices(self):
-        if self.A.cell.dimension > 1 or self.B.cell.dimension > 1:
-            raise NotImplementedError("Combining of matrices not implemented in 3D")
         if self.cell.flat and not self.symmetric:
             raise NotImplementedError("Matrices for flattened cells that are not symmetric not supported")
-        self.A.to_ufl()
-        self.B.to_ufl()
+        for f in self.factors:
+            f.to_ufl()
         oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(self.generate(), tensor=True)
         if self.flat:
             cell = self.unflat_cell
         else:
             cell = self.cell
         top = cell.to_fiat().get_topology()
-        for dim in top.keys():
-            total_dim = sum(dim) if self.flat else dim
-            a_ents = self.A.cell.get_topology()[dim[0]].keys()
-            b_ents = self.B.cell.get_topology()[dim[1]].keys()
-            ents = [(a, b) for a in a_ents for b in b_ents]
-            comp_os = cell.component_orientations()
-            for e, (a, b) in enumerate(ents):
-                ent_dofs = self.entity_dofs[total_dim][self.ent_mapping[dim][(a, b)]]
-                if len(ent_dofs) >= 1:
-                    sub_mat = oriented_mats_by_entity[dim][e]
-                    a_mat = self.A.matrices[dim[0]][a]
-                    a_ent_ids = self.A.entity_ids[dim[0]][a]
-                    b_mat = self.B.matrices[dim[1]][b]
-                    b_ent_ids = self.B.entity_ids[dim[1]][b]
-
-                    os = [(o_a, o_b) for o_a in a_mat.keys() for o_b in b_mat.keys()]
-                    for o in os:
-                        a_sub_mat = a_mat[o[0]][np.ix_(a_ent_ids, a_ent_ids)]
-                        b_sub_mat = b_mat[o[1]][np.ix_(b_ent_ids, b_ent_ids)]
-                        if self.mat_transformer is not None:
-                            o_classes = (self.A.cell.group.get_member_by_val(o[0]), self.B.cell.group.get_member_by_val(o[1]))
-                            combined_sub_mat = self.mat_transformer(a_sub_mat, b_sub_mat, o_classes)
-                        else:
-                            combined_sub_mat = np.kron(a_sub_mat, b_sub_mat)
-                        new_o = comp_os[dim][o]
-                        sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat)
-                        # sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.eye(np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat).shape[0])
+        if len(self.factors) == 2:
+            for dim in top.keys():
+                total_dim = sum(dim) if self.flat else dim
+                f_ents = [f.cell.get_topology()[d].keys() for f, d in zip(self.factors, dim)]
+                ents = list(product(*(f_ents)))
+                comp_os = cell.component_orientations()
+                for e, sub_ents in enumerate(ents):
+                    ent_dofs = self.entity_dofs[total_dim][self.ent_mapping[dim][sub_ents]]
+                    if len(ent_dofs) >= 1:
+                        sub_mat = oriented_mats_by_entity[dim][e]
+                        mats = [f.matrices[d][ent] for f, d, ent in zip(self.factors, dim, sub_ents)]
+                        ent_ids = [f.entity_dofs[d][ent] for f, d, ent in zip(self.factors, dim, sub_ents)]
+                        os = list(product(*([mat.keys() for mat in mats])))
+                        for o in os:
+                            sub_mats = [mat[o_f][np.ix_(ent_id, ent_id)] for mat, o_f, ent_id in zip(mats, o, ent_ids)]
+                            if self.mat_transformer is not None:
+                                o_classes = (f.cell.group.get_member_by_val(o_f) for f, o_f in zip(self.factors, o))
+                                combined_sub_mat = self.mat_transformer(*sub_mats, o_classes)
+                            else:
+                                combined_sub_mat = reduce(lambda acc, x: np.kron(acc, x), sub_mats)
+                            new_o = comp_os[dim][o]
+                            if new_o in sub_mat.keys():
+                                sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat)
+                            # sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.eye(np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat).shape[0])
 
         if self.cell.flat:
             oriented_mats_by_entity = flatten_dictionary(oriented_mats_by_entity)
@@ -112,39 +112,41 @@ def setup_matrices(self):
         self.reversed_matrices = self.reverse_dof_perms(self.matrices)
 
     def generate(self):
-        a_dofs = self.A.generate()
-        b_dofs = self.B.generate()
-        a_ent_assocs, _, _ = self.A._entity_associations(a_dofs, overall=False)
-        b_ent_assocs, _, _ = self.B._entity_associations(b_dofs, overall=False)
+        dofs = [f.generate() for f in self.factors]
+        ent_assocs = [f._entity_associations(dofs_f, overall=False)[0] for f, dofs_f in zip(self.factors, dofs)]
         if self.flat:
             top = self.unflat_cell.to_fiat().get_topology()
         else:
             top = self.cell.to_fiat().get_topology()
-        self.entity_dofs = {}
+        self.entity_dofs = defaultdict(dict)
         self.ent_mapping = {}
+        self.entity_assocs = defaultdict(dict)
+        self.dof_ids = {}
         dofs = []
-        ent_counter = {}
+        ent_counter = defaultdict(lambda: 0)
         dof_counter = 0
         for dim in top.keys():
             total_dim = sum(dim) if self.flat else dim
-            ents_A = a_ent_assocs[dim[0]].keys()
-            ents_B = b_ent_assocs[dim[1]].keys()
-            if total_dim not in self.entity_dofs.keys():
-                self.entity_dofs[total_dim] = {}
-                ent_counter[total_dim] = 0
+            ents = [ent_assoc[d].keys() for ent_assoc, d in zip(ent_assocs, dim)]
+            # if total_dim not in self.entity_dofs.keys():
+            #     self.entity_dofs[total_dim] = {}
+            #     self.entity_assocs[total_dim] = {}
             self.ent_mapping[dim] = {}
             ent_list = []
-            for i, ent in enumerate([(a_e, b_e) for a_e in ents_A for b_e in ents_B]):
+            for i, ent in enumerate(list(product(*ents))):
                 self.ent_mapping[dim][ent] = i + ent_counter[total_dim] if self.flat else ent
-                self.entity_dofs[total_dim][self.ent_mapping[dim][ent]] = tuple()
+                self.entity_dofs[total_dim][self.ent_mapping[dim][ent]] = []
                 ent_list += [ent]
-            for a_e, b_e in ent_list:
-                a_dofs = [d for dofs in a_ent_assocs[dim[0]][a_e].values() for d in dofs]
-                b_dofs = [d for dofs in b_ent_assocs[dim[1]][b_e].values() for d in dofs]
-                new_dofs = [(a, b) for a in a_dofs for b in b_dofs]
+            for es in ent_list:
+                e_dofs = [[d for dofs in ent_assoc[d][e].values() for d in dofs] for ent_assoc, d, e in zip(ent_assocs, dim, es)]
+                new_dofs = list(product(*e_dofs))
                 dofs += new_dofs
-                self.entity_dofs[total_dim][self.ent_mapping[dim][(a_e, b_e)]] = [i + dof_counter for i in range(len(new_dofs))]
-                dof_counter += len(new_dofs)
+                dof_gens = "(" + "*".join([",".join(list(ent_assoc[d][e].keys())) for ent_assoc, d, e in zip(ent_assocs, dim, es)]) + ")"
+                self.entity_assocs[total_dim][self.ent_mapping[dim][es]] = {dof_gens: new_dofs}
+                self.entity_dofs[total_dim][self.ent_mapping[dim][es]] += [i + dof_counter for i in range(len(new_dofs))]
+                for d in new_dofs:
+                    self.dof_ids[d] = dof_counter
+                    dof_counter += 1
                 ent_counter[total_dim] += 1
 
         return dofs
@@ -163,10 +165,10 @@ def __add__(self, other):
         return EnrichedElement(self, other, symmetric=self.symmetric and other.symmetric, matrices=self.apply_matrices or other.apply_matrices)
 
     def flatten(self):
-        return TensorProductTriple(self.A, self.B, flat=True, symmetric=self.symmetric, matrices=self.apply_matrices)
+        return TensorProductTriple(*self.factors, flat=True, symmetric=self.symmetric, matrices=self.apply_matrices)
 
     def unflatten(self):
-        return TensorProductTriple(self.A, self.B, flat=False, symmetric=self.symmetric, matrices=self.apply_matrices)
+        return TensorProductTriple(*self.factors, flat=False, symmetric=self.symmetric, matrices=self.apply_matrices)
 
 
 def compute_matrix_transform(trace, cell, o):
@@ -188,13 +190,17 @@ def compute_matrix_transform(trace, cell, o):
 class HDiv(TensorProductTriple):
 
     def __init__(self, tensor_element):
+        self.base_element = tensor_element
         self.gem_transformer, self.mat_transformer = self.select_fuse_hdiv_transformer(tensor_element)
         self.trace = TrHDiv
-        super(HDiv, self).__init__(tensor_element.A, tensor_element.B, tensor_element.flat, tensor_element.symmetric, tensor_element.matrices)
-        self.spaces[1] = TrHDiv
+        super(HDiv, self).__init__(*tensor_element.factors, flat=tensor_element.flat, symmetric=tensor_element.symmetric, matrices=tensor_element.matrices)
+        # self.spaces = (self.spaces[0], TrHDiv, self.spaces[2])
 
     def to_ufl(self):
-        return HDivElement(super(HDiv, self).to_ufl(), self.gem_transformer)
+        return HDivElement(super(HDiv, self).to_ufl(), transform=self.gem_transformer)
+
+    def repr(self):
+        return "HDiv(" + super(HDiv, self).repr() + ")"
 
     def select_fuse_hdiv_transformer(self, element):
         # Assume: something x interval
@@ -211,14 +217,15 @@ def select_fuse_hdiv_transformer(self, element):
             # Make the scalar value the right hand rule normal on the
             # y-aligned edges.
             cell = element.sub_elements[1].cell
-            bv = cell.basis_vectors()[0]
-            return lambda v: [gem.Product(gem.Literal(bv[0]), v), gem.Zero()], lambda m_a, m_b, o: np.kron(transform(cell, o[1]) @ m_a, m_b)
+            bv = cell.basis_vectors()[0][0]
+            return lambda v: [gem.Product(gem.Literal(bv), v), gem.Zero()], lambda m_a, m_b, o: np.kron(transform(cell, o[1]) @ m_a, m_b)
         elif ks == (1, 0):
             # Make the scalar value the upward-pointing normal on the
             # x-aligned edges.
             cell = element.sub_elements[0].cell
             bv = cell.basis_vectors()[0][0]
-            return lambda v: [gem.Zero(), gem.Product(gem.Literal(bv), v)], lambda m_a, m_b, o: np.kron(m_a, transform(cell, o[0]) @ m_b)
+            return lambda v: [gem.Zero(), v], lambda m_a, m_b, o: np.kron(m_a, transform(cell, o[0]) @ m_b)
+            # return lambda v: [gem.Zero(), gem.Product(gem.Literal(bv), v)], lambda m_a, m_b, o: np.kron(m_a, transform(cell, o[0]) @ m_b)
         # elif ks == (2, 0):
         #     # Same for 3D, so z-plane.
         #     return lambda v: [gem.Zero(), gem.Zero(), v]
@@ -241,17 +248,26 @@ def select_fuse_hdiv_transformer(self, element):
             raise NotImplementedError("Unexpected original mapping!")
             assert False, "Unexpected form degree combination!"
 
+    def flatten(self):
+        return HDiv(self.base_element.flatten())
+
+    def unflatten(self):
+        return HDiv(self.base_element.unflatten())
+
 
 class HCurl(TensorProductTriple):
 
     def __init__(self, tensor_element):
         self.gem_transformer, self.mat_transformer = self.select_fuse_hcurl_transformer(tensor_element)
         self.trace = TrHCurl
-        super(HDiv, self).__init__(tensor_element.A, tensor_element.B, tensor_element.flat, tensor_element.symmetric, tensor_element.matrices)
-        self.spaces[1] = TrHCurl
+        super(HDiv, self).__init__(*tensor_element.factors, tensor_element.flat, tensor_element.symmetric, tensor_element.matrices)
+        self.spaces = (self.spaces[0], TrHCurl, self.spaces[2])
 
     def to_ufl(self):
-        return HCurlElement(super(HDiv, self).to_ufl(), self.gem_transformer)
+        return HCurlElement(super(HCurl, self).to_ufl(), self.gem_transformer)
+
+    def repr(self):
+        return "HCurl(" + super(HCurl, self).repr() + ")"
 
     def select_fuse_hcurl_transformer(element):
         import gem
@@ -292,3 +308,9 @@ def select_fuse_hcurl_transformer(element):
         else:
             raise NotImplementedError("Unexpected original mapping!")
             assert False, "Unexpected original mapping!"
+
+    def flatten(self):
+        return HCurl(self.base_element.flatten())
+
+    def unflatten(self):
+        return HCurl(self.base_element.unflatten())
diff --git a/fuse/triples.py b/fuse/triples.py
index 2ac7e04d..38450063 100644
--- a/fuse/triples.py
+++ b/fuse/triples.py
@@ -7,6 +7,8 @@
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import CiarletElement
 from FIAT.reference_element import ufc_cell
+from functools import cache
+from itertools import product
 import matplotlib as mpl
 mpl.use('Agg')
 import matplotlib.pyplot as plt
@@ -15,7 +17,6 @@
 import warnings
 import numpy as np
 import scipy
-from functools import cache
 
 
 def compute_form_degree(cell, spaces):
@@ -76,11 +77,11 @@ def setup_ids_and_nodes(self):
         value_shape = self.get_value_shape()
         top = self.ref_el.get_topology()
         min_ids = self.cell.get_starter_ids()
-        entity_ids = {}
+        entity_dofs = {}
         nodes = []
 
         for dim in sorted(top):
-            entity_ids[dim] = {i: [] for i in top[dim]}
+            entity_dofs[dim] = {i: [] for i in top[dim]}
 
         self.dof_id_to_fiat_id = {}
         entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])]
@@ -90,17 +91,17 @@ def setup_ids_and_nodes(self):
             for i in range(len(dofs)):
                 if entity[1] == dofs[i].cell_defined_on.id - min_ids[dim]:
                     self.dof_id_to_fiat_id[dofs[i].id] = counter
-                    entity_ids[dim][dofs[i].cell_defined_on.id - min_ids[dim]].append(counter)
+                    entity_dofs[dim][dofs[i].cell_defined_on.id - min_ids[dim]].append(counter)
                     nodes.append(dofs[i].convert_to_fiat(self.ref_el, degree, value_shape))
                     counter += 1
         self.nodes = nodes
         # for i in range(4):
-        #     entity_ids[2][i] = [entity_ids[2][i][-1]] + entity_ids[2][i][:-1]
-        return entity_ids, nodes
+        #     entity_dofs[2][i] = [entity_dofs[2][i][-1]] + entity_dofs[2][i][:-1]
+        return entity_dofs, nodes
 
     def setup_matrices(self):
-        # self.matrices_by_entity = self.make_entity_dense_matrices(self.ref_el, self.entity_ids, self.nodes, self.poly_set)
-        matrices, entity_perms, pure_perm = self.make_dof_perms(self.ref_el, self.entity_ids, self.nodes, self.poly_set)
+        # self.matrices_by_entity = self.make_entity_dense_matrices(self.ref_el, self.entity_dofs, self.nodes, self.poly_set)
+        matrices, entity_perms, pure_perm = self.make_dof_perms(self.ref_el, self.entity_dofs, self.nodes, self.poly_set)
         reversed_matrices = self.reverse_dof_perms(matrices)
         if self.perm:
             self.pure_perm = pure_perm
@@ -159,7 +160,7 @@ def get_dof_info(self, dof, tikz=True):
         return center, colours[tikz][dof.cell_defined_on.dimension]
 
     def get_value_shape(self):
-        # TODO Shape should be specificed somewhere else probably
+        # TODO Shape should be specified somewhere else probably
         if self.spaces[0].set_shape:
             return (self.cell.get_spatial_dimension(),)
         else:
@@ -170,7 +171,7 @@ def to_ufl(self):
             # set up for eventual conversion to FIAT if not already done
             self.ref_el = self.cell.to_fiat()
             self.poly_set = self.spaces[0].to_ON_polynomial_set(self.ref_el)
-            self.entity_ids, self.nodes = self.setup_ids_and_nodes()
+            self.entity_dofs, self.nodes = self.setup_ids_and_nodes()
             self.matrices, self.reversed_matrices = self.setup_matrices()
         return FuseElement(self)
 
@@ -181,11 +182,11 @@ def to_fiat(self):
         form_degree = compute_form_degree(self.cell, self.spaces)
         degree = self.spaces[0].degree()
         # sanity check that the dofs span the space
-        original_V, original_basis = self.compute_dense_matrix(self.ref_el, self.entity_ids, self.nodes, self.poly_set)
+        original_V, original_basis = self.compute_dense_matrix(self.ref_el, self.entity_dofs, self.nodes, self.poly_set)
         if self.pure_perm:
-            dual = DualSet(self.nodes, self.ref_el, self.entity_ids, self.entity_perms)
+            dual = DualSet(self.nodes, self.ref_el, self.entity_dofs, self.entity_perms)
         else:
-            dual = DualSet(self.nodes, self.ref_el, self.entity_ids)
+            dual = DualSet(self.nodes, self.ref_el, self.entity_dofs)
         return CiarletElement(self.poly_set, dual, degree, form_degree)
 
     def to_tikz(self, show=True, scale=3):
@@ -277,8 +278,8 @@ def plot(self, filename="temp.png"):
         else:
             raise ValueError("Plotting not supported in this dimension")
 
-    def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set):
-        dual = DualSet(nodes, ref_el, entity_ids)
+    def compute_dense_matrix(self, ref_el, entity_dofs, nodes, poly_set):
+        dual = DualSet(nodes, ref_el, entity_dofs)
 
         old_coeffs = poly_set.get_coeffs()
         dualmat = dual.to_riesz(poly_set)
@@ -297,7 +298,7 @@ def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set):
                 raise np.linalg.LinAlgError("Singular Vandermonde matrix")
         return A, new_coeffs_flat
 
-    def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set):
+    def make_entity_dense_matrices(self, ref_el, entity_dofs, nodes, poly_set):
         raise NotImplementedError("This should be deprecated")
         degree = self.spaces[0].degree()
         min_ids = self.cell.get_starter_ids()
@@ -314,7 +315,7 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set):
             # dof_ids = [self.dof_id_to_fiat_id[d.id] for d in self.generate() if d.cell_defined_on == e]
             dof_ids = [d.id for d in self.generate() if d.cell_defined_on == e]
             # res_dict[dim][e_id][0] = np.eye(len(dof_ids))
-            original_V, original_basis = self.compute_dense_matrix(ref_el, entity_ids, nodes, poly_set)
+            original_V, original_basis = self.compute_dense_matrix(ref_el, entity_dofs, nodes, poly_set)
 
             for g in self.cell.group.members():
                 permuted_e, permuted_g = self.cell.permute_entities(g, dim)[e_id]
@@ -329,7 +330,7 @@ def make_mat(perm_g):
                             # , entity_o=perm_g
                             new_nodes = [d(g, entity_o=perm_g).convert_to_fiat(ref_el, degree, self.get_value_shape()) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()]
                             # new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()]
-                            transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set)
+                            transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_dofs, new_nodes, poly_set)
                             return np.matmul(transformed_basis, original_V.T)
                         temp = make_mat(permuted_g)
                         # if dim == 1 and e_id == 0:
@@ -346,7 +347,7 @@ def make_mat(perm_g):
                         res_dict[dim][e_id][val] = temp[np.ix_(dof_ids, dof_ids)]
         return res_dict
 
-    def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set):
+    def make_overall_dense_matrices(self, ref_el, entity_dofs, nodes, poly_set):
         raise NotImplementedError("this function should be unnecessary")
         min_ids = self.cell.get_starter_ids()
         dim = self.cell.dim()
@@ -354,14 +355,14 @@ def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set):
         e_id = e.id - min_ids[dim]
         res_dict = {dim: {e_id: {}}}
         degree = self.spaces[0].degree()
-        original_V, original_basis = self.compute_dense_matrix(ref_el, entity_ids, nodes, poly_set)
+        original_V, original_basis = self.compute_dense_matrix(ref_el, entity_dofs, nodes, poly_set)
         for g in self.cell.group.members():
             val = g.numeric_rep()
             if g.perm.is_Identity:
                 res_dict[dim][e_id][val] = np.eye(len(nodes))
             else:
                 new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()]
-                transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set)
+                transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_dofs, new_nodes, poly_set)
                 res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T)
         return res_dict
 
@@ -411,7 +412,7 @@ def _initialise_entity_dicts(self, dofs, tensor=False):
         flat_by_entity = {}
         cell = self.cell
         if tensor:
-            dims = [(a_d, b_d) for a_d in range(cell.A.dimension + 1) for b_d in range(cell.B.dimension + 1)]
+            dims = list(product(*(range(f.dimension + 1) for f in cell.factors)))
         else:
             dims = [i for i in range(cell.dimension + 1)]
         for dim in dims:
@@ -431,14 +432,14 @@ def _initialise_entity_dicts(self, dofs, tensor=False):
                     flat_by_entity[dim][e_id][val] = []
         return oriented_mats_by_entity, flat_by_entity
 
-    def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set):
+    def make_dof_perms(self, ref_el, entity_dofs, nodes, poly_set):
         dofs = self.generate()
         min_ids = self.cell.get_starter_ids()
         entity_associations, pure_perm, sub_pure_perm = self._entity_associations(dofs)
         # if pure_perm is False:
         #    #TODO think about where this call goes
         #    return self.matrices_by_entity, None, pure_perm
-        #    return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), None, pure_perm
+        #    return self.make_overall_dense_matrices(ref_el, entity_dofs, nodes, poly_set), None, pure_perm
 
         oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(dofs)
         # for each entity, look up generation on that entity and permute the
@@ -564,7 +565,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set):
         oriented_mats_overall = oriented_mats_by_entity[dim][0]
         if pure_perm and sub_pure_perm:
             for val, mat in oriented_mats_overall.items():
-                cell_dofs = entity_ids[dim][0]
+                cell_dofs = entity_dofs[dim][0]
                 flat_by_entity[dim][e_id][val] = perm_matrix_to_perm_array(mat[np.ix_(cell_dofs, cell_dofs)])
             return oriented_mats_by_entity, flat_by_entity, True
 
@@ -694,21 +695,21 @@ def generate(self, cell, space, id_counter):
         self.dof_ids = [dof.id for dof in self.ls]
         return self.ls
 
-    def make_entity_ids(self):
+    def make_entity_dofs(self):
         dofs = self.ls
-        entity_ids = {}
+        entity_dofs = {}
         min_ids = dofs[0].cell.get_starter_ids()
 
         top = dofs[0].cell.get_topology()
 
         for dim in sorted(top):
-            entity_ids[dim] = {i: [] for i in top[dim]}
+            entity_dofs[dim] = {i: [] for i in top[dim]}
 
         for i in range(len(dofs)):
             entity = dofs[i].cell_defined_on
             dim = entity.dim()
-            entity_ids[dim][entity.id - min_ids[dim]].append(i)
-        return entity_ids
+            entity_dofs[dim][entity.id - min_ids[dim]].append(i)
+        return entity_dofs
 
     def __repr__(self):
         repr_str = "DOFGen("
diff --git a/fuse/utils.py b/fuse/utils.py
index 73725dea..3d2fb4f9 100644
--- a/fuse/utils.py
+++ b/fuse/utils.py
@@ -2,6 +2,8 @@
 import sympy as sp
 import math
 
+_SYMBOLS = tuple(sp.Symbol(s) for s in ("x", "y", "z"))
+
 
 def fold_reduce(func_list, *prev):
     """
@@ -49,7 +51,7 @@ def tabulate_sympy(expr, pts):
     # returns: evaluation of expr at pts
     res = np.zeros((pts.shape[0],) + (expr.shape[-1],))
     i = 0
-    syms = ["x", "y", "z"]
+    syms = _SYMBOLS
     for pt in pts:
         if not hasattr(pt, "__iter__"):
             pt = (pt,)
@@ -105,3 +107,11 @@ def orientation_value(identity_arg, perm_arg):
         identity.remove(perm[i])
         val += loc * math.factorial(len(perm) - i - 1)
     return val
+
+
+def as_tuple(expr):
+    if isinstance(expr, tuple):
+        return expr
+    if isinstance(expr, list):
+        return tuple(expr)
+    return (expr,)
diff --git a/pyproject.toml b/pyproject.toml
index cc59f90d..db00e3aa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -53,4 +53,4 @@ testpaths = [
 [tool.coverage.run]
 include=[
     "fuse/*",
-]
\ No newline at end of file
+]
diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py
index 85d00b3f..ccf7006c 100644
--- a/test/test_2d_examples_docs.py
+++ b/test/test_2d_examples_docs.py
@@ -39,7 +39,7 @@ def construct_dg0_integral(edge=None):
     return dg0
 
 
-def construct_dg1_integral(cell=None):
+def construct_dg1_integral():
     edge = Point(1, [Point(0), Point(0)], vertex_num=2)
     x = sp.Symbol("x")
     xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))]
@@ -375,6 +375,7 @@ def test_nd_example():
     for dof in ned.generate():
         assert [np.allclose(1, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 1
         assert [np.allclose(0, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 2
+    ned.to_fiat()
 
 
 def construct_rt(tri=None):
diff --git a/test/test_3d_examples_docs.py b/test/test_3d_examples_docs.py
index b7ce26ed..ec8e593f 100644
--- a/test/test_3d_examples_docs.py
+++ b/test/test_3d_examples_docs.py
@@ -680,11 +680,11 @@ def test_tet_nd():
 def construct_V(elem):
     nodes = elem.nodes
     ref_el = elem.ref_el
-    entity_ids = elem.entity_ids
+    entity_dofs = elem.entity_dofs
     poly_set = elem.poly_set
     from FIAT.dual_set import DualSet
 
-    dual = DualSet(nodes, ref_el, entity_ids)
+    dual = DualSet(nodes, ref_el, entity_dofs)
 
     old_coeffs = poly_set.get_coeffs()
     dualmat = dual.to_riesz(poly_set)
diff --git a/test/test_cells.py b/test/test_cells.py
index d31ad091..35f5a9fc 100644
--- a/test/test_cells.py
+++ b/test/test_cells.py
@@ -1,13 +1,13 @@
 from fuse import *
 from firedrake import *
-from fuse.cells import ufc_triangle
+from fuse.cells import ufc_triangle, ufc_tetrahedron
 import pytest
 import numpy as np
 from FIAT.reference_element import default_simplex
 from test_convert_to_fiat import helmholtz_solve
 
 
-@pytest.fixture(scope='module', params=[0, 1, 2])
+@pytest.fixture(scope='module', params=[0, 1, 2, 3])
 def C(request):
     dim = request.param
     if dim == 0:
@@ -16,6 +16,8 @@ def C(request):
         return Point(1, [Point(0), Point(0)], vertex_num=2)
     elif dim == 2:
         return polygon(3)
+    elif dim == 3:
+        return make_tetrahedron()
 
 
 def test_vertices(C):
@@ -246,7 +248,7 @@ def test_new_connectivity(cell):
 def test_compare_tris():
     fuse_tet = polygon(3)
     ufc_tet = ufc_triangle()
-    fiat_tet = ufc_simplex(2)
+    fiat_tet = default_simplex(2)
 
     print(fiat_tet.get_topology())
     print(fuse_tet.get_topology())
@@ -277,7 +279,7 @@ def test_compare_tets():
     # perm = tet.group.get_member([1, 2, 0, 3])
     fuse_tet = tet
     ufc_tet = ufc_tetrahedron()
-    fiat_tet = ufc_simplex(3)
+    fiat_tet = default_simplex(3)
     # breakpoint()
     print(fiat_tet.get_topology())
     print(fuse_tet.get_topology())
diff --git a/test/test_construction.py b/test/test_construction.py
index 4415236e..9250a74c 100644
--- a/test/test_construction.py
+++ b/test/test_construction.py
@@ -121,28 +121,32 @@ def test_polynomial_poisson_solve(deg):
     assert np.allclose(res, 0)
 
 
-# def test_plane():
-#     from fuse import make_tetrahedron
-#     cell = make_tetrahedron()
-#     verts = cell.ordered_vertex_coords()
-#     res = check_below_plane(verts[1], verts[2], verts[3], (verts[1] + verts[2] + verts[3])/3)
-#     print(res)
-
-
-# def test_check_line():
-#     from fuse import polygon
-#     cell = polygon(3)
-#     verts = np.array(sorted(cell.ordered_vertex_coords()))
-#     midpoint = (verts[1] + verts[2])/2
-#     midpoint1 = (verts[0] + verts[2])/2
-#     assert check_below_line(verts[0], midpoint, (0, 0)) == 0
-#     assert check_on_line(verts[0], midpoint, (0, 0))
-#     assert check_on_line(verts[1], verts[2], midpoint)
-#     assert not check_on_line(verts[1], verts[2], midpoint1)
-
-#     assert check_below_line(verts[0], midpoint, (-0.5, 0)) == -1
-#     assert check_below_line(verts[0], midpoint, (0, -0.5)) == 1
-
-#     assert check_below_line(verts[1], midpoint1, verts[0]) == 1
-
-# test_construction3d(1,3, 2)
+def test_ned3():
+    nd3_pt = periodic_table(1, 3, 1, 3)
+    pt_gen = nd3_pt.DOFGenerator[1].x[0].triple.DOFGenerator[0].g1.members()
+    nd3_pt.to_fiat()
+    from test_3d_examples_docs import construct_tet_ned_2nd_kind_3
+    nd3_mn = construct_tet_ned_2nd_kind_3()
+    mn_gen = nd3_mn.DOFGenerator[1].x[0].triple.DOFGenerator[0].g1.members()
+    nd3_mn.to_fiat()
+
+    print([nd3_pt.dofs[i].id for i in range(24, 30)])
+    print([nd3_mn.dofs[i].id for i in range(24, 30)])
+    # for i in range(24, 30):
+    #     print(i)
+    #     print(pt_gen[i - 24], nd3_pt.dofs[i])
+    #     print(mn_gen[i - 24], nd3_mn.dofs[i])
+    
+    def permute_face(elem, o):
+        dof_ids = [d.id for d in elem.dofs]
+        transform_mat = elem.matrices[2][0][o.numeric_rep()]
+        transformed = np.matmul(transform_mat, dof_ids)
+        return transformed[24:30]
+    for o in mn_gen:
+        print(o)
+        print(o.numeric_rep())
+    for o in pt_gen:
+        print(o)
+        print(o.numeric_rep())
+        # print(permute_face(nd3_pt, o))
+        # print(permute_face(nd3_mn, o))
\ No newline at end of file
diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py
index 9121d19c..c1bb0287 100644
--- a/test/test_convert_to_fiat.py
+++ b/test/test_convert_to_fiat.py
@@ -365,10 +365,10 @@ def test_entity_perms(elem_gen, cell):
 
 @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1),
                                                     (create_dg1, "DG", 1),
-                                                    pytest.param(construct_dg0_integral, "DG", 0, marks=pytest.mark.xfail(reason='Passes locally, fails in CI, probably same as dg2')),
+                                                    (construct_dg0_integral, "DG", 0),
                                                     (construct_dg1_integral, "DG", 1),
                                                     (construct_dg2_integral, "DG", 2),
-                                                    pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')),
+                                                    (create_dg2, "DG", 2),
                                                     (create_cg2, "CG", 2)
                                                     ])
 def test_1d(elem_gen, elem_code, deg):
@@ -566,7 +566,7 @@ def test_poisson_analytic(params, elem_gen):
 
 
 @pytest.mark.parametrize(['elem_gen'],
-                         [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='Issue with cell/mesh'))])
+                         [(create_cg1_quad_tensor,), (create_cg1_quad,)])
 def test_quad(elem_gen):
     elem = elem_gen()
     r = 0
@@ -899,10 +899,10 @@ def test_basis_funcs_gen(form_num):
     for v in basis_funcs[:1]:
         print(v)
         vec = as_tensor(sp.lambdify(symbols, v)(x_m[0], x_m[1], x_m[2])[:, 0])
-        min_id1 = min([v for e in elem.entity_ids[2].values() for v in e])
-        max_id1 = max([v for e in elem.entity_ids[2].values() for v in e]) + 1
-        min_id2 = min([v for e in elem2.entity_ids[2].values() for v in e])
-        max_id2 = max([v for e in elem2.entity_ids[2].values() for v in e]) + 1
+        min_id1 = min([v for e in elem.entity_dofs[2].values() for v in e])
+        max_id1 = max([v for e in elem.entity_dofs[2].values() for v in e]) + 1
+        min_id2 = min([v for e in elem2.entity_dofs[2].values() for v in e])
+        max_id2 = max([v for e in elem2.entity_dofs[2].values() for v in e]) + 1
 
         res = assemble(interpolate(vec, V)).dat.data
         res2 = assemble(interpolate(vec, V2)).dat.data
diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py
index 317010f1..e6270a48 100644
--- a/test/test_tensor_prod.py
+++ b/test/test_tensor_prod.py
@@ -3,7 +3,8 @@
 from fuse import *
 from firedrake import *
 from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg0_integral, construct_dg1_integral
-from test_convert_to_fiat import create_cg2, create_dg0
+from test_convert_to_fiat import create_cg2, create_dg0, helmholtz_solve as helmholtz_solve2
+from fuse.tensor_products import HDiv as HDiv_fuse
 # from test_convert_to_fiat import create_cg1
 
 
@@ -22,6 +23,10 @@ def create_cg3_interval(cell=None):
                                                 DOFGenerator(interior, S2, S1)])
     return cg
 
+def rt1_quad():
+    cg1 = construct_cg1()
+    dg0 = construct_dg0_integral()
+    return HDiv_fuse(tensor_product(cg1, dg0).flatten()) + HDiv_fuse(tensor_product(dg0, cg1).flatten())
 
 def helmholtz_solve(mesh, V):
     u = TrialFunction(V)
@@ -103,6 +108,72 @@ def test_helmholtz(elem_gen, elem_code, deg, conv_rate):
     print("convergence order:", conv)
     assert (np.array(conv) > conv_rate).all()
 
+def project_expr(mesh, U, expr):
+    x = SpatialCoordinate(mesh)
+    f = assemble(project(expr(x), U))
+    out = Function(U)
+    u = TrialFunction(U)
+    v = TestFunction(U)
+    a = inner(u, v)*dx
+    L = inner(f, v)*dx
+    solve(a == L, out)
+    res = sqrt(assemble(dot(out - func, out - func) * dx))
+    return res
+
+@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(rt1_quad, "RT", 1, 0.8)])
+def test_project_vec_quad(elem_gen, elem_code, deg, conv_rate):
+    vals = range(3, 6)
+    function = lambda x, i: cos((3/4)*pi*x[i])
+    expr = lambda x: as_vector([function(x, 0), function(x, 1)])
+    res = []
+    res_ufc = []
+    for r in vals:
+        mesh_fuse = UnitSquareMesh(2**r, 2**r, use_fuse=True)
+        U = FunctionSpace(mesh_fuse, elem_gen().to_ufl())
+        res += [project_expr(mesh_fuse, U, expr)]
+
+    print("l2 error norms:", res)
+    res = np.array(res)
+    conv = np.log2(res[:-1] / res[1:])
+    print("convergence order:", conv)
+
+    assert (np.array(conv) > conv_rate).all()
+
+
+@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.8),
+                                                                         (create_cg2, "CG", 2, 3.8),
+                                                                         (create_cg3_interval, "CG", 3, 4.8)])
+def test_helmholtz_3d(elem_gen, elem_code, deg, conv_rate):
+    vals = range(2, 4)
+    res_ufc = []
+    res_fuse = []
+    for r in vals:
+        m = UnitSquareMesh(2**r, 2**r, quadrilateral=True, use_fuse=True)
+        mesh_fuse = ExtrudedMesh(m, 2**r)
+
+        A = elem_gen()
+        B = elem_gen()
+        C = elem_gen()
+        elem = tensor_product(tensor_product(A, B).flatten(), C)
+
+        U1 = FunctionSpace(mesh_fuse, elem.to_ufl())
+        res_fuse += [helmholtz_solve2(U1, mesh_fuse)]
+
+        m = UnitSquareMesh(2**r, 2**r, quadrilateral=True)
+        mesh_ufc = ExtrudedMesh(m, 2**r)
+        U2 = FunctionSpace(mesh_ufc, elem_code, deg)
+        res_ufc += [helmholtz_solve2(U2, mesh_ufc)]
+    print("l2 error norms:", res_ufc)
+    res_ufc = np.array(res_ufc)
+    conv_ufc = np.log2(res_ufc[:-1] / res_ufc[1:])
+    print("convergence order:", conv_ufc)
+    print("l2 error norms:", res_fuse)
+    res_fuse = np.array(res_fuse)
+    conv_fuse = np.log2(res_fuse[:-1] / res_fuse[1:])
+    print("convergence order:", conv_fuse)
+    # assert (np.array(conv_fuse) > conv_rate).all()
+    # assert (np.array(conv_ufc) > conv_rate).all()
+
 
 def test_on_quad_mesh():
     quadrilateral = True
@@ -166,6 +237,70 @@ def test_quad_mesh_helmholtz(elem_gen, elem_code, deg, conv_rate):
     assert (np.array(conv) > conv_rate).all()
 
 
+@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.7),
+                                                                         (create_cg2, "CG", 2, 3.8),
+                                                                         (create_cg3_interval, "CG", 3, 4.8)])
+def test_ext_mesh_helmholtz_3d(elem_gen, elem_code, deg, conv_rate):
+    vals = range(2, 4)
+    res_fuse = []
+    res_fire = []
+    for r in vals:
+        mesh_fuse = ExtrudedMesh(UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=True, use_fuse=True), 2**r)
+        A = elem_gen()
+        B = elem_gen()
+        C = elem_gen()
+        elem = symmetric_tensor_product(A, B, C, matrices=False)
+        U = FunctionSpace(mesh_fuse, elem.to_ufl())
+        res_fuse += [helmholtz_solve2(U, mesh_fuse)]
+
+        mesh_ufc = ExtrudedMesh(UnitSquareMesh(2 ** r, 2 ** r), 2**r)
+        U = FunctionSpace(mesh_ufc, elem_code, deg)
+        res_fire += [helmholtz_solve2(U, mesh_ufc)]
+    print("Fuse l2 error norms:", res_fuse)
+    res = np.array(res_fuse)
+    conv = np.log2(res[:-1] / res[1:])
+    print("Fuse convergence order:", conv)
+
+    print("FIAT l2 error norms:", res_fire)
+    res = np.array(res_fire)
+    conv = np.log2(res[:-1] / res[1:])
+    print("Fiat convergence order:", conv)
+    assert (np.array(conv) > conv_rate).all()
+    assert (np.array(conv) > conv_rate).all()
+
+
+@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.7),
+                                                                         (create_cg2, "CG", 2, 3.8),
+                                                                         (create_cg3_interval, "CG", 3, 4.8)])
+def test_quad_mesh_helmholtz_3d(elem_gen, elem_code, deg, conv_rate):
+    vals = range(2, 4)
+    res_fuse = []
+    res_fire = []
+    for r in vals:
+        mesh_fuse = UnitCubeMesh(2 ** r, 2 ** r, 2 ** r, hexahedral=True, use_fuse=True)
+        A = elem_gen()
+        B = elem_gen()
+        C = elem_gen()
+        elem = symmetric_tensor_product(A, B, C).flatten()
+        U = FunctionSpace(mesh_fuse, elem.to_ufl())
+        res_fuse += [helmholtz_solve2(U, mesh_fuse)]
+
+        mesh_ufc = UnitCubeMesh(2 ** r, 2 ** r, 2 ** r, hexahedral=True)
+        U = FunctionSpace(mesh_ufc, elem_code, deg)
+        res_fire += [helmholtz_solve2(U, mesh_ufc)]
+    print("Fuse l2 error norms:", res_fuse)
+    res = np.array(res_fuse)
+    conv = np.log2(res[:-1] / res[1:])
+    print("Fuse convergence order:", conv)
+
+    print("FIAT l2 error norms:", res_fire)
+    res = np.array(res_fire)
+    conv = np.log2(res[:-1] / res[1:])
+    print("Fiat convergence order:", conv)
+    assert (np.array(conv) > conv_rate).all()
+    assert (np.array(conv) > conv_rate).all()
+
+
 @pytest.mark.parametrize(["A", "B", "res"], [(Point(0), line(), False),
                                              (line(), line(), True),
                                              (polygon(3), line(), False),])
@@ -179,6 +314,20 @@ def test_flattening(A, B, res):
         cell.construct_fuse_rep()
 
 
+@pytest.mark.parametrize(["A", "B", "C"], [(line(), line(), line()),])
+def test_creation(A, B, C):
+    tensor_cell_2d = TensorProductPoint(A, B)
+    tensor_cell_2d.to_ufl()
+    tensor_cell_2d.to_fiat()
+    flat_tensor_cell_2d = tensor_cell_2d.flatten()
+    print(flat_tensor_cell_2d)
+    tensor_cell_3d = TensorProductPoint(A, B, C)
+    tensor_cell_3d.to_ufl()
+    tensor_cell_3d.to_fiat()
+    flat_tensor_cell_3d = tensor_cell_3d.flatten()
+    print(flat_tensor_cell_3d)
+
+
 def test_cg1_dg0():
     A = construct_cg1()
     B = construct_dg1_integral()
@@ -186,12 +335,12 @@ def test_cg1_dg0():
     ba = tensor_product(B, A).flatten()
     combined = ab + ba
     combined.symmetric = True
-    mesh1 = UnitSquareMesh(2, 2, quadrilateral=True)
+    mesh1 = UnitSquareMesh(2, 2, quadrilateral=True, use_fuse=True)
     V = FunctionSpace(mesh1, combined.to_ufl())
     ab = tensor_product(A, B)
     ba = tensor_product(B, A)
     combined = ab + ba
-    m = UnitIntervalMesh(2)
+    m = UnitIntervalMesh(2, use_fuse=True)
     mesh2 = ExtrudedMesh(m, 2)
     V2 = FunctionSpace(mesh2, combined.to_ufl())
     # CG_1 = FiniteElement("CG", "interval", 1)
@@ -236,7 +385,7 @@ def test_cg1_dg0():
 
 
 def test_trace_galerkin_projection():
-    mesh = UnitSquareMesh(10, 10, quadrilateral=True)
+    mesh = UnitSquareMesh(10, 10, quadrilateral=True, use_fuse=True)
 
     x, y = SpatialCoordinate(mesh)
     A = construct_cg1()
@@ -245,7 +394,7 @@ def test_trace_galerkin_projection():
     elem2 = tensor_product(B, A).flatten()
 
     # Define the Trace Space
-    T = FunctionSpace(mesh, elem.to_ufl() + elem2.to_ufl())
+    T = FunctionSpace(mesh, (elem + elem2).to_ufl())
 
     # Define trial and test functions
     lambdar = TrialFunction(T)
@@ -253,7 +402,7 @@ def test_trace_galerkin_projection():
 
     # Define right hand side function
 
-    V = FunctionSpace(mesh, "CG", 1)
+    V = FunctionSpace(mesh, tensor_product(A, construct_cg1()).flatten().to_ufl())
     f = Function(V)
     f.interpolate(cos(x*pi*2)*cos(y*pi*2))
 
@@ -274,57 +423,50 @@ def test_trace_galerkin_projection():
 
 
 def test_hdiv():
-    from fuse.tensor_products import HDiv
+    # from fuse.tensor_products import HDiv
     np.set_printoptions(linewidth=90, precision=4, suppress=True)
     m = UnitIntervalMesh(2)
     mesh = ExtrudedMesh(m, 2)
-    # CG_1 = FiniteElement("CG", "interval", 1)
-    # DG_0 = FiniteElement("DG", "interval", 0)
-    # cg1 = construct_cg1()
-    # dg0 = construct_dg0_integral()
-    # p1p0 = HDiv(tensor_product(cg1, dg0))
-    # P1P0 = TensorProductElement(CG_1, DG_0)
-    # RT_horiz = HDivElement(p1p0.to_ufl(), transform=hdiv_transform(p1p0))
-    # RT_horiz = p1p0.to_ufl()
-    # RT_horiz = HDivElement(P1P0)
-    # p0p1 = HDiv(tensor_product(dg0, cg1))
-    # P0P1 = TensorProductElement(DG_0, CG_1)
-    # RT_vert = p0p1.to_ufl()
-    # RT_vert = HDivElement(P0P1)
-    # elt = RT_horiz
-    # + RT_vert
-    # + RT_vert
-    # mesh = UnitSquareMesh(1, 1, quadrilateral=True)
-    A = construct_cg1()
-    B = construct_dg0_integral()
-    non_sym1 = tensor_product(A, B).flatten()
-    print(non_sym1.matrices[1])
-    # .flatten()
-    non_sym2 = tensor_product(B, A).flatten()
-    combined = non_sym1 + non_sym2
-    combined = combined
-    combined.symmetric = True
-    elt = HDiv(combined).to_ufl()
+    CG_1 = FiniteElement("CG", "interval", 1)
+    DG_0 = FiniteElement("DG", "interval", 0)
+    cg1 = construct_cg1()
+    dg0 = construct_dg0_integral()
+    p1p0 = HDiv_fuse(tensor_product(cg1, dg0)) + HDiv_fuse(tensor_product(dg0, cg1))
+    P1P0 = TensorProductElement(CG_1, DG_0)
+    RT_horiz = HDivElement(P1P0)
+    P0P1 = TensorProductElement(DG_0, CG_1)
+
+    RT_vert = HDivElement(P0P1)
+    elt = RT_horiz + RT_vert
+    elt2 = p1p0.to_ufl()
+    m = UnitIntervalMesh(2, use_fuse=True)
+    mesh2 = ExtrudedMesh(m, 2)
+    # mesh2 = UnitSquareMesh(2, 2, quadrilateral=True, use_fuse=True)
+    # A = construct_cg1()
+    # B = construct_dg0_integral()
+    # non_sym1 = tensor_product(A, B)
+    # # .flatten()
+    # non_sym2 = tensor_product(B, A)
+    # combined = HDiv(non_sym1) + HDiv(non_sym2)
+    # combined = combined
+    # combined.symmetric = True
+    # elt = combined.flatten().to_ufl()
     V = FunctionSpace(mesh, elt)
-    u = TrialFunction(V)
-    v = TestFunction(V)
-    f = Function(V)
-    x, y = SpatialCoordinate(mesh)
-    # f_vec = as_vector(((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2), (1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2)))
-    f_vec = as_vector((2, 3))
-    f = project(f_vec, V)
-    a = (inner(grad(u), grad(v)) + inner(u, v)) * dx
-    L = inner(f, v) * dx
-    u = Function(V)
-    solve(a == L, u)
+    V2 = FunctionSpace(mesh2, elt2)
+    for V, mesh in zip([V, V2], (mesh, mesh2)):
+        u = TrialFunction(V)
+        v = TestFunction(V)
+        f = Function(V)
+        x, y = SpatialCoordinate(mesh)
+        # f_vec = as_vector(((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2), (1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2)))
+        f_vec = as_vector((2, 3))
+        f = project(f_vec, V)
+        a = (inner(grad(u), grad(v)) + inner(u, v)) * dx
+        L = inner(f, v) * dx
+        u = Function(V)
+        solve(a == L, u)
+        print(u.dat.data)
     breakpoint()
-    # f.interpolate(cos(x*pi*2)*cos(y*pi*2))
-    # V.finat_element.basis_evaluation(1, [(0, 0)])
-    # tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)])
-    # for ent, arr in tabulation.items():
-    #     print(ent)
-    #     for comp in arr:
-    #         print(comp[0], comp[1])
 
 
 def test_transforms():
@@ -351,7 +493,6 @@ def test_transforms():
     print(HDiv(tensor_product(cg1, rev_dg0))(v))
     print(HDiv(tensor_product(dg0, rev_cg1))(v))
     print(HDiv(tensor_product(rev_cg1, dg0))(v))
-    breakpoint()
 
 
 def test_sum_fac():
@@ -384,20 +525,18 @@ def test_sum_fac():
         assert (kernel_vanilla.flop_count / kernel_spectral.flop_count) > 3
 
 
-@pytest.mark.xfail(reason="3D tensor products not implemented")
+# @pytest.mark.xfail(reason="3D tensor products not implemented")
 def test_sum_fac_3d():
     # In 2d we have O(N_q^3N_i^6) -> O(p^9)
     # Sum factorisation gains 2 factors so we expect O(p^7)
     # For CG3 p = 3 so it should be 9x faster - seems that it is faster than this in regular firedrake
-    mesh = ExtrudedMesh(UnitSquareMesh(10, 10, quadrilateral=True), 10)
+    mesh = ExtrudedMesh(UnitSquareMesh(10, 10, use_fuse=True), 10)
     A = create_cg3_interval()
     B = create_cg3_interval()
     C = create_cg3_interval()
-    elem = tensor_product(tensor_product(A, B).flatten(), C)
-    mesh2 = UnitSquareMesh(10, 10, quadrilateral=True)
-    C = create_cg3_interval()
-    D = create_cg3_interval()
-    elem2 = symmetric_tensor_product(C, D).flatten()
+    elem = tensor_product(A, B, C)
+    mesh2 = UnitCubeMesh(10, 10, 10, hexahedral=True, use_fuse=True)
+    elem2 = symmetric_tensor_product(A, B, C).flatten()
     V = FunctionSpace(mesh, elem.to_ufl())
     V1 = FunctionSpace(mesh, "CG", 3)
     V2 = FunctionSpace(mesh2, elem2.to_ufl())