diff --git a/gem/gem.py b/gem/gem.py
index 2859ded3..23719f1e 100644
--- a/gem/gem.py
+++ b/gem/gem.py
@@ -792,10 +792,6 @@ def __new__(cls, i, j):
         assert isinstance(i, IndexBase)
         assert isinstance(j, IndexBase)
 
-        # \delta_{i,i} = 1
-        if i == j:
-            return one
-
         # Fixed indices
         if isinstance(i, int) and isinstance(j, int):
             return Literal(int(i == j))
diff --git a/tsfc/fem.py b/tsfc/fem.py
index 2411a1dd..9a9bb6aa 100644
--- a/tsfc/fem.py
+++ b/tsfc/fem.py
@@ -699,4 +699,17 @@ def compile_ufl(expression, interior_facet=False, point_sum=False, **kwargs):
     result = map_expr_dags(context.translator, expressions)
     if point_sum:
         result = [gem.index_sum(expr, context.point_indices) for expr in result]
+    # Check that any indices from the PointSet have made it into the resulting expression
+    for i in range(len(result)):
+        expansion_indices = ()
+        for index in context.point_indices:
+            if index not in result[i].free_indices:
+                expansion_indices = expansion_indices + (index,)
+        if expansion_indices != ():
+            # blow up shape to include expansion_indices
+            deltas = gem.Delta(expansion_indices[0], expansion_indices[0])
+            # Use deltas to increase the possible rank of our expression
+            for index in expansion_indices[1:]:
+                deltas *= gem.Delta(index, index)
+            result[i] = gem.Indexed(gem.ComponentTensor(result[i] * deltas, expansion_indices), expansion_indices)
     return result