-
Notifications
You must be signed in to change notification settings - Fork 25
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Switch to FInAT #86
Switch to FInAT #86
Conversation
Nascent Delta support for FInAT
Implement affine index groups and an iterator for multi-indices over sets of Index and AffineIndex.
AffineIndex and IndexIterator
Conflicts: - tsfc/driver.py - tsfc/kernel_interface.py
Conflicts: tsfc/fem.py
Conflicts: tsfc/driver.py tsfc/fem.py
Conflicts: tsfc/fem.py tsfc/geometric.py tsfc/kernel_interface.py
Conflicts: gem/gem.py gem/optimise.py tsfc/kernel_interface.py
Update UFC kernel interface
self.j = j | ||
# Set up free indices | ||
free_indices = tuple(index for index in (i, j) if isinstance(index, Index)) | ||
self.free_indices = tuple(unique(free_indices)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The constant-folding of i == j
makes the uniquifying unnecessary? Unless indices hash differently from equality?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh, you care about ordering, making Delta(j, i)
and Delta(i, j)
have the same ordering of free_indices
by calling unique here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This has evolved from numpy.unique
which sorts and removes duplicates. This idea is to make x.free_indices == y.free_indices
work (without having to cast to sets first).
if len(types) == 1: | ||
cls, = types | ||
if cls.__front__ or cls.__back__: | ||
raise NotImplementedError |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Descriptive error message?
for nth_children in zip(*[e.children | ||
for e in expressions])]) | ||
|
||
raise NotImplementedError |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
And again?
|
||
|
||
def sum_factorise(sum_indices, factors): | ||
"""Optimise a tensor product throw sum factorisation. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
through?
:returns: optimised GEM expression | ||
""" | ||
if len(sum_indices) > 5: | ||
raise NotImplementedError("Too many indices for sum factorisation!") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess if we care, we can probably figure out a way to cast this as a shortest path problem.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not sure what you mean. And the number of indices to sum over is never more than three after delta elimination. 3! = 6 which is quite small.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, so we don't care. What I meant was, you can probably recast this minimisation over permutations as some kind of graph problem...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
doi:10.1109/18.910572 does something like that, but that approach only yields an easy optimum if the factor graph is cycle-free which isn't the case for us. Cycle-breaking techniques can produce correct factorisations in this case, with no guarantee about optimality, so it's not really helpful for us.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK.
@convert.register(ufl.VectorElement) | ||
def convert_vectorelement(element): | ||
scalar_element = create_element(element.sub_elements()[0]) | ||
return finat.TensorFiniteElement(scalar_element, (element.num_sub_elements(),)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Presumably element.reference_value_shape()
would work in this case too?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, because if IIRC ufl.TensorElement
asserts that the inner element is scalar, while ufl.VectorElement
gladly accepts an RT element for example.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh yeah
# This destroys the structure of the quadrature points, but since | ||
# this code path is only used to implement CellCoordinate in facet | ||
# integrals, hopefully it does not matter much. | ||
point_shape = tuple(index.extent for index in ps.indices) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we think of a test case where doing so would mean that evaluating CellCoordinate
would become the leading-order cost?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't think so. Leading order terms contain both basis function indices and quadrature indices; this one has no basis functions, so it involves quadrature indices only.
"Q": None, | ||
} | ||
"""A :class:`.dict` mapping UFL element family names to their | ||
FIAT-equivalent constructors. If the value is ``None``, the UFL |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
finat?
|
||
|
||
def create_element(element): | ||
"""Create a FIAT element (suitable for tabulating with) given a UFL element. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
finat?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Overall looks good.
Closes #75. Requires firedrakeproject/firedrake#964.