diff --git a/finat/element_factory.py b/finat/element_factory.py
index 64053bbc7..85830e2e5 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -149,13 +149,14 @@ def convert(element, **kwargs):
 @convert.register(finat.ufl.FiniteElement)
 def convert_finiteelement(element, **kwargs):
     cell = as_fiat_cell(element.cell)
-    if element.family() == "Quadrature":
+    if element.family() in {"Quadrature", "Boundary Quadrature"}:
         degree = element.degree()
-        scheme = element.quadrature_scheme()
+        scheme = element.quadrature_scheme() or "default"
         if degree is None or scheme is None:
             raise ValueError("Quadrature scheme and degree must be specified!")
 
-        return finat.make_quadrature_element(cell, degree, scheme), set()
+        codim = 1 if element.family() == "Boundary Quadrature" else 0
+        return finat.make_quadrature_element(cell, degree, scheme, codim), set()
     lmbda = supported_elements[element.family()]
     if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}:
         lmbda = None
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
index 0efea5e00..24f4a3d7d 100644
--- a/finat/fiat_elements.py
+++ b/finat/fiat_elements.py
@@ -116,9 +116,8 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
                 continue
 
             derivative = sum(alpha)
-            table_roll = fiat_table.reshape(
-                space_dimension, value_size, len(ps.points)
-            ).transpose(1, 2, 0)
+            shp = (space_dimension, value_size, *ps.points.shape[:-1])
+            table_roll = np.moveaxis(fiat_table.reshape(shp), 0, -1)
 
             exprs = []
             for table in table_roll:
diff --git a/finat/finiteelementbase.py b/finat/finiteelementbase.py
index 64a6399d2..b7953510d 100644
--- a/finat/finiteelementbase.py
+++ b/finat/finiteelementbase.py
@@ -77,6 +77,9 @@ def entity_closure_dofs(self):
         element.'''
         return self._entity_closure_dofs
 
+    def is_dg(self):
+        return self.entity_dofs() == self.entity_closure_dofs()
+
     @cached_property
     def _entity_support_dofs(self):
         esd = {}
diff --git a/finat/point_set.py b/finat/point_set.py
index 2a437c9c1..9e91faae4 100644
--- a/finat/point_set.py
+++ b/finat/point_set.py
@@ -33,8 +33,7 @@ def points(self):
     @property
     def dimension(self):
         """Point dimension."""
-        _, dim = self.points.shape
-        return dim
+        return self.points.shape[-1]
 
     @property
     @abc.abstractmethod
@@ -129,8 +128,7 @@ def points(self):
 
     @cached_property
     def indices(self):
-        N, _ = self._points_expr.shape
-        return (gem.Index(extent=N),)
+        return tuple(gem.Index(extent=N) for N in self._points_expr.shape[:-1])
 
     @cached_property
     def expression(self):
@@ -159,7 +157,7 @@ def points(self):
 
     @cached_property
     def indices(self):
-        return (gem.Index(extent=len(self.points)),)
+        return tuple(gem.Index(extent=N) for N in self.points.shape[:-1])
 
     @cached_property
     def expression(self):
@@ -224,3 +222,58 @@ def almost_equal(self, other, tolerance=1e-12):
             len(self.factors) == len(other.factors) and \
             all(s.almost_equal(o, tolerance=tolerance)
                 for s, o in zip(self.factors, other.factors))
+
+
+class FacetPointSet(AbstractPointSet):
+    """A point set on facets.
+
+    A FacetPointSet is constructed by mapping a lower-dimensional PointSet
+    onto all facets of a higher-dimensional cell.
+
+    Parameters
+    ----------
+    cell : FIAT.Cell
+        The cell.
+    ps : PointSet
+        A lower-dimensional point set.
+
+    """
+    def __init__(self, cell, ps):
+        self.cell = cell
+        self.ps = ps
+
+    def __repr__(self):
+        return f"{type(self).__name__}({self.ps!r})"
+
+    @cached_property
+    def entities(self):
+        """The list of all cell entites matching the reference point set dimension."""
+        to_int = lambda x: sum(x) if isinstance(x, tuple) else x
+        top = self.cell.topology
+        return [(dim, entity)
+                for dim in sorted(top)
+                for entity in sorted(top[dim])
+                if to_int(dim) == self.ps.dimension]
+
+    @cached_property
+    def points(self):
+        """The array with the reference points mapped onto each facet."""
+        ref_pts = self.ps.points
+        pts = [self.cell.get_entity_transform(dim, entity)(ref_pts)
+               for dim, entity in self.entities]
+        return numpy.concatenate(pts)
+
+    @cached_property
+    def indices(self):
+        """An Index tuple of the facet index and the reference point indices."""
+        return (gem.Index(extent=len(self.entities)), *self.ps.indices)
+
+    @cached_property
+    def expression(self):
+        raise NotImplementedError("Symbolic point expression is not yet implemented for FacetPointSet.")
+
+    def almost_equal(self, other, tolerance=1e-12):
+        """Approximate numerical equality of point sets"""
+        return type(self) is type(other) and \
+            self.cell == other.cell and \
+            self.ps.almost_equal(other.ps, tolerance=tolerance)
diff --git a/finat/quadrature.py b/finat/quadrature.py
index 10e4048b9..c182a8a9a 100644
--- a/finat/quadrature.py
+++ b/finat/quadrature.py
@@ -1,6 +1,6 @@
 import hashlib
 from abc import ABCMeta, abstractmethod
-from functools import cached_property, reduce
+from functools import cached_property
 
 import gem
 import numpy
@@ -163,4 +163,4 @@ def point_set(self):
 
     @cached_property
     def weight_expression(self):
-        return reduce(gem.Product, (q.weight_expression for q in self.factors))
+        return gem.Product(*(q.weight_expression for q in self.factors))
diff --git a/finat/quadrature_element.py b/finat/quadrature_element.py
index 3f17ec399..c7b784688 100644
--- a/finat/quadrature_element.py
+++ b/finat/quadrature_element.py
@@ -1,5 +1,4 @@
-from finat.point_set import UnknownPointSet
-from functools import reduce
+from finat.point_set import UnknownPointSet, FacetPointSet
 
 import numpy
 
@@ -13,7 +12,7 @@
 from finat.quadrature import make_quadrature, AbstractQuadratureRule
 
 
-def make_quadrature_element(fiat_ref_cell, degree, scheme="default"):
+def make_quadrature_element(fiat_ref_cell, degree, scheme="default", codim=0):
     """Construct a :class:`QuadratureElement` from a given a reference
     element, degree and scheme.
 
@@ -23,9 +22,16 @@ def make_quadrature_element(fiat_ref_cell, degree, scheme="default"):
         integrate exactly.
     :param scheme: The quadrature scheme to use - e.g. "default",
         "canonical" or "KMV".
+    :param codim: The codimension of the quadrature scheme.
     :returns: The appropriate :class:`QuadratureElement`
     """
-    rule = make_quadrature(fiat_ref_cell, degree, scheme=scheme)
+    if codim > 0:
+        sd = fiat_ref_cell.get_spatial_dimension()
+        rule_ref_cell = fiat_ref_cell.construct_subelement(sd - codim)
+    else:
+        rule_ref_cell = fiat_ref_cell
+
+    rule = make_quadrature(rule_ref_cell, degree, scheme=scheme)
     return QuadratureElement(fiat_ref_cell, rule)
 
 
@@ -42,8 +48,6 @@ def __init__(self, fiat_ref_cell, rule):
         self.cell = fiat_ref_cell
         if not isinstance(rule, AbstractQuadratureRule):
             raise TypeError("rule is not an AbstractQuadratureRule")
-        if fiat_ref_cell.get_spatial_dimension() != rule.point_set.dimension:
-            raise ValueError("Cell dimension does not match rule's point set dimension")
         self._rule = rule
 
     @cached_property
@@ -64,10 +68,18 @@ def formdegree(self):
 
     @cached_property
     def _entity_dofs(self):
-        # Inspired by ffc/quadratureelement.py
+        top = self.cell.get_topology()
         entity_dofs = {dim: {entity: [] for entity in entities}
-                       for dim, entities in self.cell.get_topology().items()}
-        entity_dofs[self.cell.get_dimension()] = {0: list(range(self.space_dimension()))}
+                       for dim, entities in top.items()}
+        ps = self._rule.point_set
+        num_pts = len(ps.points)
+        to_int = lambda x: sum(x) if isinstance(x, tuple) else x
+        cur = 0
+        for dim in sorted(top):
+            if to_int(dim) == ps.dimension:
+                for entity in sorted(top[dim]):
+                    entity_dofs[dim][entity].extend(range(cur, cur + num_pts))
+                    cur += num_pts
         return entity_dofs
 
     def entity_dofs(self):
@@ -76,9 +88,15 @@ def entity_dofs(self):
     def space_dimension(self):
         return numpy.prod(self.index_shape, dtype=int)
 
+    @cached_property
+    def _point_set(self):
+        ps = self._rule.point_set
+        sd = self.cell.get_spatial_dimension()
+        return ps if ps.dimension == sd else FacetPointSet(self.cell, ps)
+
     @property
     def index_shape(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         return tuple(index.extent for index in ps.indices)
 
     @property
@@ -87,7 +105,7 @@ def value_shape(self):
 
     @cached_property
     def fiat_equivalent(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         if isinstance(ps, UnknownPointSet):
             raise ValueError("A quadrature element with rule with runtime points has no fiat equivalent!")
         weights = getattr(self._rule, 'weights', None)
@@ -107,8 +125,13 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         :param ps: the point set object.
         :param entity: the cell entity on which to tabulate.
         '''
-        if entity is not None and entity != (self.cell.get_dimension(), 0):
-            raise ValueError('QuadratureElement does not "tabulate" on subentities.')
+        rule_dim = self._rule.point_set.dimension
+        if entity is None:
+            entity = (rule_dim, 0)
+        entity_dim, entity_id = entity
+        if entity_dim != rule_dim:
+            raise ValueError(f"Cannot tabulate QuadratureElement of dimension {rule_dim}"
+                             f" on subentities of dimension {entity_dim}.")
 
         if order:
             raise ValueError("Derivatives are not defined on a QuadratureElement.")
@@ -117,24 +140,25 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
             raise ValueError("Mismatch of quadrature points!")
 
         # Return an outer product of identity matrices
-        multiindex = self.get_indices()
-        product = reduce(gem.Product, [gem.Delta(q, r)
-                                       for q, r in zip(ps.indices, multiindex)])
+        basis_indices = self.get_indices()
+        point_indices = ps.indices
+        if len(basis_indices) > len(point_indices):
+            point_indices = (entity_id, *point_indices)
+        delta = gem.Delta(point_indices, basis_indices)
 
-        dim = self.cell.get_spatial_dimension()
-        return {(0,) * dim: gem.ComponentTensor(product, multiindex)}
+        sd = self.cell.get_spatial_dimension()
+        return {(0,) * sd: gem.ComponentTensor(delta, basis_indices)}
 
     def point_evaluation(self, order, refcoords, entity=None):
         raise NotImplementedError("QuadratureElement cannot do point evaluation!")
 
     @property
     def dual_basis(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         multiindex = self.get_indices()
         # Evaluation matrix is just an outer product of identity
         # matrices, evaluation points are just the quadrature points.
-        Q = reduce(gem.Product, (gem.Delta(q, r)
-                                 for q, r in zip(ps.indices, multiindex)))
+        Q = gem.Delta(ps.indices, multiindex)
         Q = gem.ComponentTensor(Q, multiindex)
         return Q, ps
 
diff --git a/finat/sympy2gem.py b/finat/sympy2gem.py
index 29add8760..9d613a307 100644
--- a/finat/sympy2gem.py
+++ b/finat/sympy2gem.py
@@ -27,13 +27,13 @@ def sympy2gem_expr(node, self):
 @sympy2gem.register(sympy.Add)
 @sympy2gem.register(symengine.Add)
 def sympy2gem_add(node, self):
-    return reduce(gem.Sum, map(self, node.args))
+    return gem.Sum(*map(self, node.args))
 
 
 @sympy2gem.register(sympy.Mul)
 @sympy2gem.register(symengine.Mul)
 def sympy2gem_mul(node, self):
-    return reduce(gem.Product, map(self, node.args))
+    return gem.Product(*map(self, node.args))
 
 
 @sympy2gem.register(sympy.Pow)
diff --git a/finat/tensor_product.py b/finat/tensor_product.py
index f0fd58477..959213d0e 100644
--- a/finat/tensor_product.py
+++ b/finat/tensor_product.py
@@ -1,4 +1,3 @@
-from functools import reduce
 from itertools import chain, product
 from operator import methodcaller
 
@@ -32,11 +31,11 @@ def __init__(self, factors):
 
     @cached_property
     def cell(self):
-        return TensorProductCell(*[fe.cell for fe in self.factors])
+        return TensorProductCell(*(fe.cell for fe in self.factors))
 
     @cached_property
     def complex(self):
-        return TensorProductCell(*[fe.complex for fe in self.factors])
+        return TensorProductCell(*(fe.complex for fe in self.factors))
 
     @property
     def degree(self):
@@ -69,7 +68,7 @@ def space_dimension(self):
 
     @property
     def index_shape(self):
-        return tuple(chain(*[fe.index_shape for fe in self.factors]))
+        return tuple(chain.from_iterable(fe.index_shape for fe in self.factors))
 
     @property
     def value_shape(self):
@@ -117,24 +116,20 @@ def _merge_evaluations(self, factor_results):
         # multiindex describing the value shape of the subelement.
         zetas = [fe.get_value_indices() for fe in self.factors]
 
+        multiindex = tuple(chain(*alphas, *zetas))
         result = {}
         for derivative in range(order + 1):
             for Delta in mis(dimension, derivative):
                 # Split the multiindex for the subelements
                 deltas = [Delta[s] for s in dim_slices]
-                # GEM scalars (can have free indices) for collecting
-                # the contributions from the subelements.
-                scalars = []
-                for fr, delta, alpha, zeta in zip(factor_results, deltas, alphas, zetas):
-                    # Turn basis shape to free indices, select the
-                    # right derivative entry, and collect the result.
-                    scalars.append(gem.Indexed(fr[delta], alpha + zeta))
-                # Multiply the values from the subelements and wrap up
-                # non-point indices into shape.
-                result[Delta] = gem.ComponentTensor(
-                    reduce(gem.Product, scalars),
-                    tuple(chain(*(alphas + zetas)))
-                )
+                # Multiply the values from the subelements
+                # Turn basis shape to free indices, select the
+                # right derivative entry.
+                scalar = gem.Product(*(gem.Indexed(fr[delta], alpha + zeta)
+                                       for fr, delta, alpha, zeta
+                                       in zip(factor_results, deltas, alphas, zetas)))
+                # Wrap up non-point indices into shape.
+                result[Delta] = gem.ComponentTensor(scalar, multiindex)
         return result
 
     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
@@ -179,12 +174,11 @@ def dual_basis(self):
         # Naming as _merge_evaluations above
         alphas = [factor.get_indices() for factor in self.factors]
         zetas = [factor.get_value_indices() for factor in self.factors]
-        # Index the factors by so that we can reshape into index-shape
-        # followed by value-shape
-        qis = [q[alpha + zeta] for q, alpha, zeta in zip(qs, alphas, zetas)]
         Q = gem.ComponentTensor(
-            reduce(gem.Product, qis),
-            tuple(chain(*(alphas + zetas)))
+            # Index the factors by basis function and component so that we can reshape
+            # into index-shape followed by value-shape
+            gem.Product(*(q[alpha + zeta] for q, alpha, zeta in zip(qs, alphas, zetas))),
+            tuple(chain(*alphas, *zetas))
         )
         return Q, ps
 
diff --git a/finat/tensorfiniteelement.py b/finat/tensorfiniteelement.py
index c0a8aa91e..b8c50e158 100644
--- a/finat/tensorfiniteelement.py
+++ b/finat/tensorfiniteelement.py
@@ -1,4 +1,3 @@
-from functools import reduce
 from itertools import chain
 
 import numpy
@@ -134,8 +133,7 @@ def _tensorise(self, scalar_evaluation):
         tensor_vi = tuple(gem.Index(extent=d) for d in self._shape)
 
         # Couple new basis function and value indices
-        deltas = reduce(gem.Product, (gem.Delta(j, k)
-                                      for j, k in zip(tensor_i, tensor_vi)))
+        deltas = gem.Delta(tensor_i, tensor_vi)
 
         if self._transpose:
             index_ordering = tensor_i + scalar_i + tensor_vi + scalar_vi
@@ -163,8 +161,7 @@ def dual_basis(self):
         tensor_i = tuple(gem.Index(extent=d) for d in self._shape)
         tensor_vi = tuple(gem.Index(extent=d) for d in self._shape)
         # Couple new basis function and value indices
-        deltas = reduce(gem.Product, (gem.Delta(j, k)
-                                      for j, k in zip(tensor_i, tensor_vi)))
+        deltas = gem.Delta(tensor_i, tensor_vi)
         if self._transpose:
             index_ordering = tensor_i + scalar_i + tensor_vi + scalar_vi
         else:
diff --git a/gem/gem.py b/gem/gem.py
index 8369b6f75..5e653e7b1 100644
--- a/gem/gem.py
+++ b/gem/gem.py
@@ -16,6 +16,7 @@
 
 from abc import ABCMeta
 from itertools import chain
+from functools import reduce
 from operator import attrgetter
 from numbers import Integral, Number
 
@@ -324,7 +325,12 @@ def __init__(self, name, shape, dtype=None):
 class Sum(Scalar):
     __slots__ = ('children',)
 
-    def __new__(cls, a, b):
+    def __new__(cls, *args):
+        try:
+            a, b = args
+        except ValueError:
+            # Handle more than two arguments
+            return reduce(Sum, args)
         assert not a.shape
         assert not b.shape
 
@@ -345,7 +351,12 @@ def __new__(cls, a, b):
 class Product(Scalar):
     __slots__ = ('children',)
 
-    def __new__(cls, a, b):
+    def __new__(cls, *args):
+        try:
+            a, b = args
+        except ValueError:
+            # Handle more than two arguments
+            return reduce(Product, args)
         assert not a.shape
         assert not b.shape
 
@@ -960,6 +971,9 @@ class Delta(Scalar, Terminal):
     __back__ = ('dtype',)
 
     def __new__(cls, i, j, dtype=None):
+        if isinstance(i, tuple) and isinstance(j, tuple):
+            # Handle multiindices
+            return Product(*map(Delta, i, j))
         assert isinstance(i, IndexBase)
         assert isinstance(j, IndexBase)
 
@@ -968,9 +982,15 @@ def __new__(cls, i, j, dtype=None):
             return one
 
         # Fixed indices
-        if isinstance(i, int) and isinstance(j, int):
+        if isinstance(i, Integral) and isinstance(j, Integral):
             return one if i == j else Zero()
 
+        if isinstance(i, Integral):
+            return Indexed(Literal(numpy.eye(j.extent)[i]), (j,))
+
+        if isinstance(j, Integral):
+            return Indexed(Literal(numpy.eye(i.extent)[j]), (i,))
+
         self = super(Delta, cls).__new__(cls)
         self.i = i
         self.j = j
diff --git a/gem/optimise.py b/gem/optimise.py
index 7d6c8ecd6..c56ba7fbb 100644
--- a/gem/optimise.py
+++ b/gem/optimise.py
@@ -2,7 +2,7 @@
 expressions."""
 
 from collections import OrderedDict, defaultdict
-from functools import singledispatch, partial, reduce
+from functools import singledispatch, partial
 from itertools import combinations, permutations, zip_longest
 from numbers import Integral
 
@@ -374,7 +374,7 @@ def sum_factorise(sum_indices, factors):
 
     # Form groups by free indices
     groups = groupby(factors, key=lambda f: f.free_indices)
-    groups = [reduce(Product, terms) for _, terms in groups]
+    groups = [Product(*terms) for _, terms in groups]
 
     # Sum factorisation
     expression = None
@@ -414,7 +414,7 @@ def sum_factorise(sum_indices, factors):
 def make_sum(summands):
     """Constructs an operation-minimal sum of GEM expressions."""
     groups = groupby(summands, key=lambda f: f.free_indices)
-    summands = [reduce(Sum, terms) for _, terms in groups]
+    summands = [Sum(*terms) for _, terms in groups]
     result, flops = associate(Sum, summands)
     return result
 
@@ -623,7 +623,7 @@ def _replace_delta_delta(node, self):
         return Indexed(Identity(size), (i, j))
     else:
         def expression(index):
-            if isinstance(index, int):
+            if isinstance(index, Integral):
                 return Literal(index)
             elif isinstance(index, VariableIndex):
                 return index.expression
@@ -660,10 +660,8 @@ def _(node, self):
         # Unrolling
         summand = self(node.children[0])
         shape = tuple(index.extent for index in unroll)
-        unrolled = reduce(Sum,
-                          (Indexed(ComponentTensor(summand, unroll), alpha)
-                           for alpha in numpy.ndindex(shape)),
-                          Zero())
+        unrolled = Sum(*(Indexed(ComponentTensor(summand, unroll), alpha)
+                         for alpha in numpy.ndindex(shape)))
         return IndexSum(unrolled, tuple(index for index in node.multiindex
                                         if index not in unroll))
     else: