diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py
index 45b0287e..142dbce9 100644
--- a/FIAT/barycentric_interpolation.py
+++ b/FIAT/barycentric_interpolation.py
@@ -25,22 +25,24 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
     https://doi.org/10.1137/S0036144502417715 Eq. (4.2) & (9.4)
     """
     if pts.dtype == object:
-        from sympy import simplify
-        sp_simplify = numpy.vectorize(simplify)
+        # Do not use barycentric interpolation at unknown points
+        phi = numpy.add.outer(-nodes, pts.flatten())
+        phis = [wi * numpy.prod(phi[:i], axis=0) * numpy.prod(phi[i+1:], axis=0) for i, wi in enumerate(wts)]
+        phi = numpy.asarray(phis)
     else:
-        sp_simplify = lambda x: x
-    phi = numpy.add.outer(-nodes, pts.flatten())
-    with numpy.errstate(divide='ignore', invalid='ignore'):
-        numpy.reciprocal(phi, out=phi)
-        numpy.multiply(phi, wts[:, None], out=phi)
-        numpy.multiply(1.0 / numpy.sum(phi, axis=0), phi, out=phi)
-    phi[phi != phi] = 1.0
-    phi = phi.reshape(-1, *pts.shape[:-1])
+        # Use the second barycentric interpolation formula
+        phi = numpy.add.outer(-nodes, pts.flatten())
+        with numpy.errstate(divide='ignore', invalid='ignore'):
+            numpy.reciprocal(phi, out=phi)
+            numpy.multiply(phi, wts[:, None], out=phi)
+            numpy.multiply(1.0 / numpy.sum(phi, axis=0), phi, out=phi)
+        # Replace nan with one
+        phi[phi != phi] = 1.0
 
-    phi = sp_simplify(phi)
+    phi = phi.reshape(-1, *pts.shape[:-1])
     results = {(0,): phi}
     for r in range(1, order+1):
-        phi = sp_simplify(numpy.dot(dmat, phi))
+        phi = numpy.dot(dmat, phi)
         results[(r,)] = phi
     return results
 
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index b7d846f2..43d6184e 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -86,19 +86,19 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
         scale = -scale
 
     num_members = math.comb(n + dim, dim)
-    results = tuple([None] * num_members for i in range(order+1))
-    phi, dphi, ddphi = results + (None,) * (2-order)
 
     outer = lambda x, y: x[:, None, ...] * y[None, ...]
     sym_outer = lambda x, y: outer(x, y) + outer(y, x)
 
     pad_dim = dim + 2
     dX = pad_jacobian(Jinv, pad_dim)
-    phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), scale)
-    if dphi is not None:
-        dphi[0] = (phi[0] - phi[0]) * dX[0]
-    if ddphi is not None:
-        ddphi[0] = outer(dphi[0], dX[0])
+
+    phi0 = numpy.array([sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 0.0)])
+    results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)
+               for k in range(order+1)]
+
+    phi, dphi, ddphi = results + [None] * (2-order)
+    phi[0] += scale
     if dim == 0 or n == 0:
         return results
     if dim > 3 or dim < 0:
@@ -183,7 +183,7 @@ def C0_basis(dim, n, tabulations):
     # Recover facet bubbles
     for phi in tabulations:
         icur = 0
-        phi[icur] *= -1
+        phi[icur] *= -1.0
         for inext in range(1, dim+1):
             phi[icur] -= phi[inext]
         if dim == 2:
@@ -723,11 +723,15 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
     :kwarg tol: the absolute tolerance.
     :returns: a list of (weighted) characteristic functions for each subcell.
     """
-    from sympy import Piecewise
+    import gem
     sd = ref_el.get_spatial_dimension()
     top = ref_el.get_topology()
     # assert singleton point
     pt = pt.reshape((sd,))
+    if isinstance(pt[0], gem.Node):
+        import gem as backend
+    else:
+        import sympy as backend
 
     # The distance to the nearest cell is equal to the distance to the parent cell
     best = ref_el.get_parent().distance_to_point_l1(pt, rescale=True)
@@ -739,7 +743,7 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
     for cell in sorted(top[sd]):
         # Bin points based on l1 distance
         pt_near_cell = ref_el.distance_to_point_l1(pt, entity=(sd, cell), rescale=True) < tol
-        masks.append(Piecewise(*otherwise, (1.0, pt_near_cell), (0.0, True)))
+        masks.append(backend.Piecewise(*otherwise, (1.0, pt_near_cell), (0.0, True)))
         if unique:
             otherwise.append((0.0, pt_near_cell))
     # If the point is on a facet, divide the characteristic function by the facet multiplicity
diff --git a/finat/argyris.py b/finat/argyris.py
index 019a835d..e907eff8 100644
--- a/finat/argyris.py
+++ b/finat/argyris.py
@@ -124,7 +124,7 @@ def _edge_transform(V, vorder, eorder, fiat_cell, coordinate_mapping, avg=False)
             V[s, v1id] = P1 * Bnt
             V[s, v0id] = P0 * Bnt
             if k > 0:
-                V[s, s + eorder] = -1 * Bnt
+                V[s, s + eorder] = -Bnt
 
 
 class Argyris(PhysicallyMappedElement, ScalarFiatElement):
@@ -167,7 +167,7 @@ def basis_transformation(self, coordinate_mapping):
 
                 # vertex points
                 V[s, v1id] = 15/8 * Bnt
-                V[s, v0id] = -1 * V[s, v1id]
+                V[s, v0id] = -V[s, v1id]
 
                 # vertex derivatives
                 for i in range(sd):
@@ -178,7 +178,7 @@ def basis_transformation(self, coordinate_mapping):
                 tau = [Jt[0]*Jt[0], 2*Jt[0]*Jt[1], Jt[1]*Jt[1]]
                 for i in range(len(tau)):
                     V[s, v1id+3+i] = 1/32 * Bnt * tau[i]
-                    V[s, v0id+3+i] = -1 * V[s, v1id+3+i]
+                    V[s, v0id+3+i] = -V[s, v1id+3+i]
 
         # Patch up conditioning
         h = coordinate_mapping.cell_size()
diff --git a/finat/bell.py b/finat/bell.py
index 8f6b095a..d174fe16 100644
--- a/finat/bell.py
+++ b/finat/bell.py
@@ -44,7 +44,7 @@ def basis_transformation(self, coordinate_mapping):
 
             # vertex points
             V[s, v1id] = 1/21 * Bnt
-            V[s, v0id] = -1 * V[s, v1id]
+            V[s, v0id] = -V[s, v1id]
 
             # vertex derivatives
             for i in range(sd):
@@ -55,7 +55,7 @@ def basis_transformation(self, coordinate_mapping):
             tau = [Jt[0]*Jt[0], 2*Jt[0]*Jt[1], Jt[1]*Jt[1]]
             for i in range(len(tau)):
                 V[s, v1id+3+i] = 1/252 * Bnt * tau[i]
-                V[s, v0id+3+i] = -1 * V[s, v1id+3+i]
+                V[s, v0id+3+i] = -V[s, v1id+3+i]
 
         # Patch up conditioning
         h = coordinate_mapping.cell_size()
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
index 6498298b..22408cf0 100644
--- a/finat/fiat_elements.py
+++ b/finat/fiat_elements.py
@@ -1,12 +1,10 @@
 import FIAT
 import gem
 import numpy as np
-import sympy as sp
 from gem.utils import cached_property
 
 from finat.finiteelementbase import FiniteElementBase
-from finat.point_set import PointSet
-from finat.sympy2gem import sympy2gem
+from finat.point_set import PointSet, PointSingleton
 
 
 class FiatElement(FiniteElementBase):
@@ -67,10 +65,8 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         :param ps: the point set.
         :param entity: the cell entity on which to tabulate.
         '''
-        space_dimension = self._element.space_dimension()
-        value_size = np.prod(self._element.value_shape(), dtype=int)
-        fiat_result = self._element.tabulate(order, ps.points, entity)
-        result = {}
+        fiat_element = self._element
+        fiat_result = fiat_element.tabulate(order, ps.points, entity)
         # In almost all cases, we have
         # self.space_dimension() == self._element.space_dimension()
         # But for Bell, FIAT reports 21 basis functions,
@@ -78,51 +74,52 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         # basis functions, and the additional 3 are for
         # dealing with transformations between physical
         # and reference space).
-        index_shape = (self._element.space_dimension(),)
+        value_shape = self.value_shape
+        space_dimension = fiat_element.space_dimension()
+        if self.space_dimension() == space_dimension:
+            beta = self.get_indices()
+            index_shape = tuple(index.extent for index in beta)
+        else:
+            index_shape = (space_dimension,)
+            beta = tuple(gem.Index(extent=i) for i in index_shape)
+            assert len(beta) == len(self.get_indices())
+
+        zeta = self.get_value_indices()
+        basis_indices = beta + zeta
+
+        result = {}
         for alpha, fiat_table in fiat_result.items():
             if isinstance(fiat_table, Exception):
-                result[alpha] = gem.Failure(self.index_shape + self.value_shape, fiat_table)
+                result[alpha] = gem.Failure(index_shape + value_shape, fiat_table)
                 continue
 
+            point_indices = ()
+            replace_indices = ()
             derivative = sum(alpha)
-            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:
-                if derivative == self.degree and not self.complex.is_macrocell():
-                    # Make sure numerics satisfies theory
-                    exprs.append(gem.Literal(table[0]))
-                elif derivative > self.degree:
-                    # Make sure numerics satisfies theory
-                    assert np.allclose(table, 0.0)
-                    exprs.append(gem.Literal(np.zeros(self.index_shape)))
+            if derivative == self.degree and self.complex.is_simplex():
+                # Ensure a cellwise constant tabulation
+                if fiat_table.dtype == object:
+                    replace_indices = tuple((i, 0) for i in ps.expression.free_indices)
                 else:
-                    point_indices = ps.indices
-                    point_shape = tuple(index.extent for index in point_indices)
-
-                    exprs.append(gem.partial_indexed(
-                        gem.Literal(table.reshape(point_shape + index_shape)),
-                        point_indices
-                    ))
-            if self.value_shape:
-                # As above, this extent may be different from that
-                # advertised by the finat element.
-                beta = tuple(gem.Index(extent=i) for i in index_shape)
-                assert len(beta) == len(self.get_indices())
-
-                zeta = self.get_value_indices()
-                result[alpha] = gem.ComponentTensor(
-                    gem.Indexed(
-                        gem.ListTensor(np.array(
-                            [gem.Indexed(expr, beta) for expr in exprs]
-                        ).reshape(self.value_shape)),
-                        zeta),
-                    beta + zeta
-                )
+                    fiat_table = fiat_table.reshape(*index_shape, *value_shape, -1)
+                    assert np.allclose(fiat_table, fiat_table[..., 0, None])
+                    fiat_table = fiat_table[..., 0]
+            elif derivative > self.degree:
+                # Ensure a zero tabulation
+                if fiat_table.dtype != object:
+                    assert np.allclose(fiat_table, 0.0)
+                fiat_table = np.zeros(index_shape + value_shape)
             else:
-                expr, = exprs
-                result[alpha] = expr
+                point_indices = ps.indices
+
+            point_shape = tuple(i.extent for i in point_indices)
+            fiat_table = fiat_table.reshape(index_shape + value_shape + point_shape)
+            gem_table = gem.as_gem(fiat_table)
+            expr = gem.Indexed(gem_table, basis_indices + point_indices)
+            expr = gem.ComponentTensor(expr, basis_indices)
+            if replace_indices:
+                expr, = gem.optimise.remove_componenttensors((expr,), subst=replace_indices)
+            result[alpha] = expr
         return result
 
     def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
@@ -147,7 +144,20 @@ def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=Non
         esd = self.cell.construct_subelement(entity_dim).get_spatial_dimension()
         assert isinstance(refcoords, gem.Node) and refcoords.shape == (esd,)
 
-        return point_evaluation(self._element, order, refcoords, (entity_dim, entity_i))
+        # Coordinates on the reference entity (GEM)
+        Xi = tuple(gem.Indexed(refcoords, i) for i in np.ndindex(refcoords.shape))
+        ps = PointSingleton(Xi)
+        result = self.basis_evaluation(order, ps, entity=entity,
+                                       coordinate_mapping=coordinate_mapping)
+
+        # Apply symbolic simplification
+        vals = result.values()
+        vals = map(gem.optimise.ffc_rounding, vals, [1E-13]*len(vals))
+        vals = gem.optimise.constant_fold_zero(vals)
+        vals = map(gem.optimise.aggressive_unroll, vals)
+        vals = gem.optimise.remove_componenttensors(vals)
+        result = dict(zip(result.keys(), vals))
+        return result
 
     @cached_property
     def _dual_basis(self):
@@ -260,49 +270,6 @@ def mapping(self):
             return result
 
 
-def point_evaluation(fiat_element, order, refcoords, entity):
-    # Coordinates on the reference entity (SymPy)
-    esd, = refcoords.shape
-    Xi = sp.symbols('X Y Z')[:esd]
-
-    space_dimension = fiat_element.space_dimension()
-    value_size = np.prod(fiat_element.value_shape(), dtype=int)
-    fiat_result = fiat_element.tabulate(order, [Xi], entity)
-    result = {}
-    for alpha, fiat_table in fiat_result.items():
-        if isinstance(fiat_table, Exception):
-            result[alpha] = gem.Failure((space_dimension,) + fiat_element.value_shape(), fiat_table)
-            continue
-
-        # Convert SymPy expression to GEM
-        mapper = gem.node.Memoizer(sympy2gem)
-        mapper.bindings = {s: gem.Indexed(refcoords, (i,))
-                           for i, s in enumerate(Xi)}
-        gem_table = np.vectorize(mapper)(fiat_table)
-
-        table_roll = gem_table.reshape(space_dimension, value_size).transpose()
-
-        exprs = []
-        for table in table_roll:
-            exprs.append(gem.ListTensor(table.reshape(space_dimension)))
-        if fiat_element.value_shape():
-            beta = (gem.Index(extent=space_dimension),)
-            zeta = tuple(gem.Index(extent=d)
-                         for d in fiat_element.value_shape())
-            result[alpha] = gem.ComponentTensor(
-                gem.Indexed(
-                    gem.ListTensor(np.array(
-                        [gem.Indexed(expr, beta) for expr in exprs]
-                    ).reshape(fiat_element.value_shape())),
-                    zeta),
-                beta + zeta
-            )
-        else:
-            expr, = exprs
-            result[alpha] = expr
-    return result
-
-
 class Regge(FiatElement):  # symmetric matrix valued
     def __init__(self, cell, degree, **kwargs):
         super().__init__(FIAT.Regge(cell, degree, **kwargs))
diff --git a/finat/hct.py b/finat/hct.py
index f374be76..72da016c 100644
--- a/finat/hct.py
+++ b/finat/hct.py
@@ -71,7 +71,7 @@ def basis_transformation(self, coordinate_mapping):
 
             # vertex points
             V[s, v0id] = 1/5 * Bnt
-            V[s, v1id] = -1 * V[s, v0id]
+            V[s, v1id] = -V[s, v0id]
 
             # vertex derivatives
             for i in range(sd):
diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py
index df0d598a..511d3e14 100644
--- a/finat/physically_mapped.py
+++ b/finat/physically_mapped.py
@@ -89,10 +89,6 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         result = super().basis_evaluation(order, ps, entity=entity)
         return self.map_tabulation(result, coordinate_mapping)
 
-    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
-        result = super().point_evaluation(order, refcoords, entity=entity)
-        return self.map_tabulation(result, coordinate_mapping)
-
     def dual_transformation(self, Q, coordinate_mapping=None):
         M = self.basis_transformation(coordinate_mapping)
 
diff --git a/finat/point_set.py b/finat/point_set.py
index b0442de4..068f7659 100644
--- a/finat/point_set.py
+++ b/finat/point_set.py
@@ -76,7 +76,7 @@ def points(self):
 
     @cached_property
     def expression(self):
-        return gem.Literal(self.point)
+        return gem.as_gem(self.point)
 
 
 class UnknownPointsArray():
diff --git a/finat/sympy2gem.py b/finat/sympy2gem.py
index 9d613a30..440585cc 100644
--- a/finat/sympy2gem.py
+++ b/finat/sympy2gem.py
@@ -1,6 +1,5 @@
 from functools import singledispatch, reduce
 
-import numpy
 import sympy
 try:
     import symengine
@@ -130,18 +129,7 @@ def sympy2gem_le(node, self):
 @sympy2gem.register(sympy.Piecewise)
 @sympy2gem.register(symengine.Piecewise)
 def sympy2gem_conditional(node, self):
-    expr = None
-    pieces = []
-    for v, c in node.args:
-        if isinstance(c, (bool, numpy.bool, sympy.logic.boolalg.BooleanTrue)) and c:
-            expr = self(v)
-            break
-        pieces.append((v, c))
-    if expr is None:
-        expr = gem.Literal(float("nan"))
-    for v, c in reversed(pieces):
-        expr = gem.Conditional(self(c), self(v), expr)
-    return expr
+    return gem.Piecewise(*[(self(v), self(c)) for v, c in node.args])
 
 
 @sympy2gem.register(sympy.ITE)
diff --git a/gem/gem.py b/gem/gem.py
index 02aeaaf5..5c9231dc 100644
--- a/gem/gem.py
+++ b/gem/gem.py
@@ -16,7 +16,7 @@
 
 from abc import ABCMeta
 from itertools import chain, repeat
-from functools import reduce
+from functools import partial, reduce
 from operator import attrgetter
 from numbers import Integral, Number
 
@@ -36,7 +36,7 @@
            'IndexSum', 'ListTensor', 'Concatenate', 'Delta', 'OrientationVariableIndex',
            'index_sum', 'partial_indexed', 'reshape', 'view',
            'indices', 'as_gem', 'FlexiblyIndexed',
-           'Inverse', 'Solve', 'extract_type', 'uint_type']
+           'Inverse', 'Solve', 'extract_type', 'uint_type', 'Piecewise']
 
 
 uint_type = numpy.dtype(numpy.uintc)
@@ -88,6 +88,9 @@ def __getitem__(self, indices):
             indices = (indices, )
         return Indexed(self, indices)
 
+    def __neg__(self):
+        return componentwise(Product, minus, self)
+
     def __add__(self, other):
         return componentwise(Sum, self, as_gem(other))
 
@@ -95,9 +98,7 @@ def __radd__(self, other):
         return as_gem(other).__add__(self)
 
     def __sub__(self, other):
-        return componentwise(
-            Sum, self,
-            componentwise(Product, Literal(-1), as_gem(other)))
+        return componentwise(Sum, self, -as_gem(other))
 
     def __rsub__(self, other):
         return as_gem(other).__sub__(self)
@@ -124,6 +125,24 @@ def __matmul__(self, other):
     def __rmatmul__(self, other):
         return as_gem(other).__matmul__(self)
 
+    def __abs__(self):
+        return componentwise(partial(MathFunction, "abs"), self)
+
+    def __pow__(self, other):
+        return componentwise(Power, self, as_gem(other))
+
+    def __lt__(self, other):
+        return componentwise(partial(Comparison, "<"), self, as_gem(other))
+
+    def __gt__(self, other):
+        return componentwise(partial(Comparison, ">"), self, as_gem(other))
+
+    def __le__(self, other):
+        return componentwise(partial(Comparison, "<="), self, as_gem(other))
+
+    def __ge__(self, other):
+        return componentwise(partial(Comparison, ">="), self, as_gem(other))
+
     @property
     def T(self):
         i = indices(len(self.shape))
@@ -272,7 +291,6 @@ class Literal(Constant):
     __back__ = ('dtype',)
 
     def __new__(cls, array, dtype=None):
-        array = asarray(array)
         return super(Literal, cls).__new__(cls)
 
     def __init__(self, array, dtype=None):
@@ -307,6 +325,9 @@ def value(self):
     def shape(self):
         return self.array.shape
 
+    def __bool__(self):
+        return bool(self.value)
+
 
 class Variable(Terminal):
     """Symbolic variable tensor"""
@@ -553,8 +574,8 @@ def __init__(self, a, b):
         self.children = a, b
 
 
-class Conditional(Node):
-    __slots__ = ('children', 'shape')
+class Conditional(Scalar):
+    __slots__ = ('children',)
 
     def __new__(cls, condition, then, else_):
         assert not condition.shape
@@ -567,7 +588,6 @@ def __new__(cls, condition, then, else_):
 
         self = super(Conditional, cls).__new__(cls)
         self.children = condition, then, else_
-        self.shape = then.shape
         self.dtype = Node.inherit_dtype_from_children((then, else_))
         return self
 
@@ -1269,6 +1289,7 @@ def view(expression, *slices):
 
 # Static one object for quicker constant folding
 one = Literal(1)
+minus = Literal(-1)
 
 
 # Syntax sugar
@@ -1326,6 +1347,14 @@ def as_gem(expr):
         return expr
     elif isinstance(expr, Number):
         return Literal(expr)
+    elif isinstance(expr, (bool, numpy.bool)):
+        return Literal(bool(expr))
+    elif isinstance(expr, numpy.ndarray):
+        if expr.dtype == object:
+            expr = numpy.vectorize(as_gem)(expr)
+            return ListTensor(expr)
+        else:
+            return Literal(expr)
     else:
         raise ValueError("Do not know how to convert %r to GEM" % expr)
 
@@ -1360,3 +1389,31 @@ def as_gem_uint(expr):
 def extract_type(expressions, klass):
     """Collects objects of type klass in expressions."""
     return tuple(node for node in traversal(expressions) if isinstance(node, klass))
+
+
+def Piecewise(*args):
+    """Represents a piecewise function.
+
+    Parameters
+    ----------
+    *args
+        Each argument is a 2-tuple defining an expression and condition.
+
+    Returns
+    -------
+    Node
+        A nested Conditional.
+
+    """
+    expr = None
+    pieces = []
+    for v, c in args:
+        if isinstance(c, (bool, numpy.bool, Literal)) and c:
+            expr = as_gem(v)
+            break
+        pieces.append((as_gem(v), as_gem(c)))
+    if expr is None:
+        expr = Literal(float("nan"))
+    for v, c in reversed(pieces):
+        expr = Conditional(c, v, expr)
+    return expr
diff --git a/gem/optimise.py b/gem/optimise.py
index ccf0e2c1..caf254e0 100644
--- a/gem/optimise.py
+++ b/gem/optimise.py
@@ -178,10 +178,10 @@ def filtered_replace_indices(node, self, subst):
     return replace_indices(node, self, filtered_subst)
 
 
-def remove_componenttensors(expressions):
+def remove_componenttensors(expressions, subst=()):
     """Removes all ComponentTensors in multi-root expression DAG."""
     mapper = MemoizerArg(filtered_replace_indices)
-    return [mapper(expression, ()) for expression in expressions]
+    return [mapper(expression, subst) for expression in expressions]
 
 
 @singledispatch
diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py
index 17bf870c..9ad45a95 100644
--- a/test/FIAT/unit/test_macro.py
+++ b/test/FIAT/unit/test_macro.py
@@ -444,10 +444,11 @@ def test_macro_sympy(cell, element):
     variant = "spectral,alfeld"
     K = IsoSplit(cell)
     ebig = element(K, 3, variant=variant)
-    pts = get_lagrange_points(ebig.dual_basis())
+    pts = numpy.asarray(get_lagrange_points(ebig.dual_basis()))
+
+    X = tuple(sympy.Symbol(f"X[{i}]") for i in range(pts.shape[1]))
 
     dim = cell.get_spatial_dimension()
-    X = tuple(sympy.Symbol(f"X[{i}]") for i in range(dim))
     degrees = range(1, 3) if element is Lagrange else range(3)
     for degree in degrees:
         fe = element(cell, degree, variant=variant)
@@ -459,6 +460,34 @@ def test_macro_sympy(cell, element):
         assert numpy.allclose(results, tab_numpy)
 
 
+@pytest.mark.parametrize("element", (DiscontinuousLagrange, Lagrange))
+def test_macro_gem(cell, element):
+    import gem
+    from gem.interpreter import evaluate
+
+    variant = "spectral,alfeld"
+    K = IsoSplit(cell)
+    ebig = element(K, 3, variant=variant)
+    pts = numpy.asarray(get_lagrange_points(ebig.dual_basis()))
+
+    coords = gem.Variable('X', pts.shape)
+    bindings = {coords: pts}
+    index = gem.Index()
+    point = gem.partial_indexed(coords, (index,))
+    X = tuple(gem.Indexed(point, i) for i in numpy.ndindex(point.shape))
+
+    dim = cell.get_spatial_dimension()
+    degrees = range(1, 3) if element is Lagrange else range(3)
+    for degree in degrees:
+        fe = element(cell, degree, variant=variant)
+        tab_gem = fe.tabulate(0, X)[(0,) * dim]
+
+        results, = evaluate((gem.as_gem(tab_gem),), bindings=bindings)
+        results = results.arr.T
+        tab_numpy = fe.tabulate(0, pts)[(0,) * dim]
+        assert numpy.allclose(results, tab_numpy)
+
+
 @pytest.mark.parametrize("element,degree", [
     (Lagrange, 1), (Nedelec, 1), (RaviartThomas, 1), (DiscontinuousLagrange, 0),
     (Regge, 0), (HellanHerrmannJohnson, 0),