Select H(div)/H(curl) components with a Delta

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

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!"