From c84c379dda53ef505877f734498248ff422a3110 Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Tue, 2 Feb 2021 18:38:37 +0000 Subject: [PATCH] Ensure compile_ufl output has correct free indices When an expression given to compile_ufl does not have points in the input PointSet in its expression tree, the resulting gem expression has missing free indices. This attempts to fix that. --- gem/gem.py | 4 ---- tsfc/fem.py | 13 +++++++++++++ 2 files changed, 13 insertions(+), 4 deletions(-) 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