Skip to content

Commit

Permalink
Ensure compile_ufl output has correct free indices
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
ReubenHill committed Feb 2, 2021
1 parent c220702 commit c84c379
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 4 deletions.
4 changes: 0 additions & 4 deletions gem/gem.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
13 changes: 13 additions & 0 deletions tsfc/fem.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit c84c379

Please sign in to comment.