Merge pull request #87 from FInAT/ksagiyam/orientation

orientations: define entity_permutations() method for basic finite elements

finat/cube.py

@@ -1,7 +1,7 @@
 from __future__ import absolute_import, print_function, division
 
 from FIAT.reference_element import UFCHexahedron, UFCQuadrilateral
-from FIAT.reference_element import compute_unflattening_map, flatten_entities
+from FIAT.reference_element import compute_unflattening_map, flatten_entities, flatten_permutations
 from FIAT.tensor_product import FlattenedDimensions as FIAT_FlattenedDimensions
 
 from gem.utils import cached_property
@@ -49,6 +49,10 @@ def _entity_support_dofs(self):
     def entity_dofs(self):
         return self._entity_dofs
 
+    @cached_property
+    def entity_permutations(self):
+        return flatten_permutations(self.product.entity_permutations)
+
     def space_dimension(self):
         return self.product.space_dimension()
 

finat/discontinuous.py

@@ -36,6 +36,16 @@ def _entity_dofs(self):
     def entity_dofs(self):
         return self._entity_dofs
 
+    @cached_property
+    def entity_permutations(self):
+        # Return entity_permutations of the base finite element if it only
+        # has cell degrees of freedom; otherwise entity_permutations is not
+        # yet implemented for DiscontinuousElement.
+        if self.element.entity_dofs() == self.element.entity_closure_dofs():
+            return self.element.entity_permutations
+        else:
+            raise NotImplementedError(f"entity_permutations not yet implemented for a general {type(self)}")
+
     def space_dimension(self):
         return self.element.space_dimension()
 

finat/enriched.py

@@ -47,6 +47,12 @@ def entity_dofs(self):
         return concatenate_entity_dofs(self.cell, self.elements,
                                        methodcaller("entity_dofs"))
 
+    @cached_property
+    def entity_permutations(self):
+        '''Return the map of topological entities to the map of
+        orientations to permutation lists for the finite element'''
+        return concatenate_entity_permutations(self.elements)
+
     @cached_property
     def _entity_support_dofs(self):
         return concatenate_entity_dofs(self.cell, self.elements,
@@ -166,3 +172,25 @@ def concatenate_entity_dofs(ref_el, elements, method):
             for ent, off in dofs.items():
                 entity_dofs[dim][ent] += list(map(partial(add, offsets[i]), off))
     return entity_dofs
+
+
+def concatenate_entity_permutations(elements):
+    """For each dimension, for each entity, and for each possible
+    entity orientation, collect the DoF permutation lists from
+    entity_permutations dicts of elements and concatenate them.
+
+    :arg elements: subelements whose DoF permutation lists are concatenated
+    :returns: entity_permutation dict of the :class:`EnrichedElement` object
+        composed of elements.
+    """
+    permutations = {}
+    for element in elements:
+        for dim, e_o_p_map in element.entity_permutations.items():
+            dim_permutations = permutations.setdefault(dim, {})
+            for e, o_p_map in e_o_p_map.items():
+                e_dim_permutations = dim_permutations.setdefault(e, {})
+                for o, p in o_p_map.items():
+                    o_e_dim_permutations = e_dim_permutations.setdefault(o, [])
+                    offset = len(o_e_dim_permutations)
+                    o_e_dim_permutations += list(offset + q for q in p)
+    return permutations

finat/fiat_elements.py

@@ -70,6 +70,10 @@ def entity_dofs(self):
     def entity_closure_dofs(self):
         return self._element.entity_closure_dofs()
 
+    @property
+    def entity_permutations(self):
+        return self._element.entity_permutations()
+
     def space_dimension(self):
         return self._element.space_dimension()
 

finat/finiteelementbase.py

@@ -33,6 +33,31 @@ def entity_dofs(self):
         '''Return the map of topological entities to degrees of
         freedom for the finite element.'''
 
+    @property
+    def entity_permutations(self):
+        '''Returns a nested dictionary that gives, for each dimension,
+        for each entity, and for each possible entity orientation, the
+        DoF permutation array that maps the entity local DoF ordering
+        to the canonical global DoF ordering.
+
+        The entity permutations `dict` for the degree 4 Lagrange finite
+        element on the interval, for instance, is given by:
+
+        .. code-block:: python3
+
+            {0: {0: {0: [0]},
+                 1: {0: [0]}},
+             1: {0: {0: [0, 1, 2],
+                     1: [2, 1, 0]}}}
+
+        Note that there are two entities on dimension ``0`` (vertices),
+        each of which has only one possible orientation, while there is
+        a single entity on dimension ``1`` (interval), which has two
+        possible orientations representing non-reflected and reflected
+        intervals.
+        '''
+        raise NotImplementedError(f"entity_permutations not yet implemented for {type(self)}")
+
     @cached_property
     def _entity_closure_dofs(self):
         # Compute the nodes on the closure of each sub_entity.

finat/hdivcurl.py

@@ -33,6 +33,10 @@ def degree(self):
     def entity_dofs(self):
         return self.wrappee.entity_dofs()
 
+    @property
+    def entity_permutations(self):
+        return self.wrappee.entity_permutations
+
     def entity_closure_dofs(self):
         return self.wrappee.entity_closure_dofs()
 

finat/tensor_product.py

@@ -55,6 +55,10 @@ def _entity_support_dofs(self):
     def entity_dofs(self):
         return self._entity_dofs
 
+    @cached_property
+    def entity_permutations(self):
+        return compose_permutations(self.factors)
+
     def space_dimension(self):
         return numpy.prod([fe.space_dimension() for fe in self.factors])
 
@@ -208,6 +212,40 @@ def productise(factors, method):
     return dofs
 
 
+def compose_permutations(factors):
+    """For the :class:`TensorProductElement` object composed of factors,
+    construct, for each dimension tuple, for each entity, and for each possible
+    entity orientation combination, the DoF permutation list.
+
+    :arg factors: element factors.
+    :returns: entity_permutation dict of the :class:`TensorProductElement` object
+        composed of factors.
+    """
+    permutations = {}
+    for dim in product(*[fe.cell.get_topology().keys()
+                         for fe in factors]):
+        dim_permutations = []
+        e_o_p_maps = [fe.entity_permutations[d]
+                      for fe, d in zip(factors, dim)]
+        for e_tuple in product(*[sorted(e_o_p_map) for e_o_p_map in e_o_p_maps]):
+            o_p_maps = [e_o_p_map[e] for e_o_p_map, e in zip(e_o_p_maps, e_tuple)]
+            o_tuple_perm_map = {}
+            for o_tuple in product(*[o_p_map.keys() for o_p_map in o_p_maps]):
+                ps = [o_p_map[o] for o_p_map, o in zip(o_p_maps, o_tuple)]
+                shape = tuple(len(p) for p in ps)
+                size = numpy.prod(shape)
+                if size == 0:
+                    o_tuple_perm_map[o_tuple] = []
+                else:
+                    a = numpy.arange(size).reshape(shape)
+                    for i, p in enumerate(ps):
+                        a = a.swapaxes(0, i)[p, :].swapaxes(0, i)
+                    o_tuple_perm_map[o_tuple] = a.reshape(-1).tolist()
+            dim_permutations.append((e_tuple, o_tuple_perm_map))
+        permutations[dim] = dict(enumerate(v for k, v in sorted(dim_permutations)))
+    return permutations
+
+
 def factor_point_set(product_cell, product_dim, point_set):
     """Factors a point set for the product element into a point sets for
     each subelement.