diff --git a/finat/enriched.py b/finat/enriched.py index 1530911c..8c94d4a6 100644 --- a/finat/enriched.py +++ b/finat/enriched.py @@ -5,6 +5,7 @@ import FIAT import gem import numpy +from gem.interpreter import evaluate from gem.utils import cached_property from finat.finiteelementbase import FiniteElementBase @@ -243,9 +244,20 @@ def concatenate_entity_permutations(elements): def is_orthogonal(A, B): - """Test whether two elements are orthogonal.""" + """Test whether two elements map into orthogonal components. + + .. todo:: + + Read the components off the products of :class:`gem.Delta` that the + transforms select with, rather than evaluating the contraction + numerically in an element constructor. + """ if isinstance(A, (HCurlElement, HDivElement)) and isinstance(B, (HCurlElement, HDivElement)): - Amap = A.transform(gem.Literal(numpy.ones(A.wrappee.value_shape))) - Bmap = B.transform(gem.Literal(numpy.ones(B.wrappee.value_shape))) - return sum(a * b for a, b in zip(Amap, Bmap)) == gem.Literal(0.0) + dim, = A.value_shape + zeta = gem.Index(extent=dim) + a = A.transform(gem.Literal(numpy.ones(A.wrappee.value_shape)), zeta) + b = B.transform(gem.Literal(numpy.ones(B.wrappee.value_shape)), zeta) + contraction = gem.ComponentTensor(gem.Product(a, b), (zeta,)) + result, = evaluate([contraction]) + return not result.arr.any() return False diff --git a/finat/hdivcurl.py b/finat/hdivcurl.py index 3ca69774..347f25bc 100644 --- a/finat/hdivcurl.py +++ b/finat/hdivcurl.py @@ -1,3 +1,5 @@ +from functools import reduce + from FIAT.hdivcurl import Hdiv, Hcurl from FIAT.reference_element import LINE @@ -20,7 +22,8 @@ def __init__(self, wrappee, transform): self.transform = transform """A transformation applied on the scalar/vector values of the wrapped element to produce an H(div) or H(curl) conforming - element.""" + element. Takes the wrapped value and the value index of the + result, and returns the component that index selects.""" @property def cell(self): @@ -61,11 +64,11 @@ def value_shape(self): def _transform_evaluation(self, core_eval): beta = self.get_indices() zeta = self.get_value_indices() + z, = zeta def promote(table): v = gem.partial_indexed(table, beta) - u = gem.ListTensor(self.transform(v)) - return gem.ComponentTensor(gem.Indexed(u, zeta), beta + zeta) + return gem.ComponentTensor(self.transform(v, z), beta + zeta) return {alpha: promote(table) for alpha, table in core_eval.items()} @@ -83,13 +86,14 @@ def dual_basis(self): Q, x = self.wrappee.dual_basis beta = self.get_indices() zeta = self.get_value_indices() + z, = zeta # Index out the basis indices from wrapee's Q, to get # something of wrappee.value_shape, then promote to new shape # with the same transform as done for basis evaluation - Q = gem.ListTensor(self.transform(gem.partial_indexed(Q, beta))) + Q = self.transform(gem.partial_indexed(Q, beta), z) # Finally wrap up Q in shape again (now with some extra # value_shape indices) - return gem.ComponentTensor(Q[zeta], beta + zeta), x + return gem.ComponentTensor(Q, beta + zeta), x class HDivElement(WrapperElementBase): @@ -148,6 +152,21 @@ def mapping(self): return "covariant piola" +def select_component(zeta, components): + """Select one component of a vector with the vector's value index. + + :arg zeta: the value index of the vector + :arg components: the scalar components of the vector, in order + :returns: the component that ``zeta`` selects + + Selecting with a :class:`gem.Delta` keeps the components apart, so that + contracting two of these over ``zeta`` cancels the components that do not + coincide. + """ + return reduce(gem.Sum, (gem.Product(gem.Delta(zeta, i), v) + for i, v in enumerate(components))) + + def select_hdiv_transformer(element): # Assume: something x interval assert len(element.factors) == 2 @@ -161,27 +180,27 @@ def select_hdiv_transformer(element): if ks == (0, 1): # Make the scalar value the leftward-pointing normal on the # y-aligned edges. - return lambda v: [gem.Product(gem.Literal(-1), v), gem.Zero()] + return lambda v, zeta: select_component(zeta, [gem.Product(gem.Literal(-1), v), gem.Zero()]) elif ks == (1, 0): # Make the scalar value the upward-pointing normal on the # x-aligned edges. - return lambda v: [gem.Zero(), v] + return lambda v, zeta: select_component(zeta, [gem.Zero(), v]) elif ks == (2, 0): # Same for 3D, so z-plane. - return lambda v: [gem.Zero(), gem.Zero(), v] + return lambda v, zeta: select_component(zeta, [gem.Zero(), gem.Zero(), v]) elif ks == (1, 1): if element.mapping == "contravariant piola": # Pad the 2-vector normal on the "base" cell into a # 3-vector, maintaining direction. - return lambda v: [gem.Indexed(v, (0,)), - gem.Indexed(v, (1,)), - gem.Zero()] + return lambda v, zeta: select_component(zeta, [gem.Indexed(v, (0,)), + gem.Indexed(v, (1,)), + gem.Zero()]) elif element.mapping == "covariant piola": # Rotate the 2-vector tangential component on the "base" # cell 90 degrees anticlockwise into a 3-vector and pad. - return lambda v: [gem.Indexed(v, (1,)), - gem.Product(gem.Literal(-1), gem.Indexed(v, (0,))), - gem.Zero()] + return lambda v, zeta: select_component(zeta, [gem.Indexed(v, (1,)), + gem.Product(gem.Literal(-1), gem.Indexed(v, (0,))), + gem.Zero()]) else: assert False, "Unexpected original mapping!" else: @@ -202,24 +221,24 @@ def select_hcurl_transformer(element): if ks == (1, 0): # Can only be 2D. Make the scalar value the # rightward-pointing tangential on the x-aligned edges. - return lambda v: [v, gem.Zero()] + return lambda v, zeta: select_component(zeta, [v, gem.Zero()]) elif ks == (0, 1): # Can be any spatial dimension. Make the scalar value the # upward-pointing tangential. - return lambda v: [gem.Zero()] * (dim - 1) + [v] + return lambda v, zeta: select_component(zeta, [gem.Zero()] * (dim - 1) + [v]) else: assert False elif element.mapping == "covariant piola": # Second factor must be continuous interval. Just padding. - return lambda v: [gem.Indexed(v, (0,)), - gem.Indexed(v, (1,)), - gem.Zero()] + return lambda v, zeta: select_component(zeta, [gem.Indexed(v, (0,)), + gem.Indexed(v, (1,)), + gem.Zero()]) elif element.mapping == "contravariant piola": # Second factor must be continuous interval. Rotate the # 2-vector tangential component on the "base" cell 90 degrees # clockwise into a 3-vector and pad. - return lambda v: [gem.Product(gem.Literal(-1), gem.Indexed(v, (1,))), - gem.Indexed(v, (0,)), - gem.Zero()] + return lambda v, zeta: select_component(zeta, [gem.Product(gem.Literal(-1), gem.Indexed(v, (1,))), + gem.Indexed(v, (0,)), + gem.Zero()]) else: assert False, "Unexpected original mapping!"