Merge pull request #83 from firedrakeproject/ksagiyam/hex_interior_facet

hex: enable interior facet integration

FIAT/orientation_utils.py

@@ -4,6 +4,27 @@
 import numpy as np
 
 
+class Orientation(object):
+    """Base class representing unsigned integer orientations.
+
+    Orientations represented by this class are to be consistent with those used in
+    `make_entity_permutations_simplex` and `make_entity_permutations_tensorproduct`.
+
+    """
+
+    def __floordiv__(self, other):
+        raise NotImplementedError("Subclass must implement __floordiv__")
+
+    def __rfloordiv__(self, other):
+        raise NotImplementedError("Subclass must implement __rfloordiv__")
+
+    def __mod__(self, other):
+        raise NotImplementedError("Subclass must implement __mod__")
+
+    def __rmod__(self, other):
+        raise NotImplementedError("Subclass must implement __rmod__")
+
+
 def make_entity_permutations_simplex(dim, npoints):
     r"""Make orientation-permutation map for the given
     simplex dimension, dim, and the number of points along

FIAT/quadrature.py

@@ -8,10 +8,12 @@
 # Modified by David A. Ham ([email protected]), 2015
 
 import itertools
+from math import factorial
 import numpy
 from recursivenodes.quadrature import gaussjacobi, lobattogaussjacobi, simplexgausslegendre
 
 from FIAT import reference_element
+from FIAT.orientation_utils import make_entity_permutations_simplex
 
 
 def pseudo_determinant(A):
@@ -51,6 +53,7 @@ def __init__(self, ref_el, pts, wts):
         self.ref_el = ref_el
         self.pts = pts
         self.wts = wts
+        self._intrinsic_orientation_permutation_map_tuple = (None, )
 
     def get_points(self):
         return numpy.array(self.pts)
@@ -61,6 +64,32 @@ def get_weights(self):
     def integrate(self, f):
         return sum(w * f(x) for x, w in zip(self.pts, self.wts))
 
+    @property
+    def extrinsic_orientation_permutation_map(self):
+        """A map from extrinsic orientations to corresponding axis permutation matrices.
+
+        Notes
+        -----
+        result[eo] gives the physical axis-reference axis permutation matrix corresponding to
+        eo (extrinsic orientation).
+
+        """
+        return self.ref_el.extrinsic_orientation_permutation_map
+
+    @property
+    def intrinsic_orientation_permutation_map_tuple(self):
+        """A tuple of maps from intrinsic orientations to corresponding point permutations for each reference cell axis.
+
+        Notes
+        -----
+        result[axis][io] gives the physical point-reference point permutation array corresponding to
+        io (intrinsic orientation) on ``axis``.
+
+        """
+        if any(m is None for m in self._intrinsic_orientation_permutation_map_tuple):
+            raise ValueError("Must set _intrinsic_orientation_permutation_map_tuple")
+        return self._intrinsic_orientation_permutation_map_tuple
+
 
 class GaussJacobiQuadratureLineRule(QuadratureRule):
     """Gauss-Jacobi quadature rule determined by Jacobi weights a and b
@@ -71,6 +100,12 @@ def __init__(self, ref_el, m, a=0, b=0):
         pts_ref, wts_ref = gaussjacobi(m, a, b)
         pts, wts = map_quadrature(pts_ref, wts_ref, Ref1, ref_el)
         QuadratureRule.__init__(self, ref_el, pts, wts)
+        # Set _intrinsic_orientation_permutation_map_tuple.
+        dim = 1
+        a = numpy.zeros((factorial(dim + 1), m), dtype=int)
+        for io, perm in make_entity_permutations_simplex(dim, m).items():
+            a[io, perm] = range(m)
+        self._intrinsic_orientation_permutation_map_tuple = (a, )
 
 
 class GaussLobattoLegendreQuadratureLineRule(QuadratureRule):

FIAT/reference_element.py

@@ -29,8 +29,11 @@
 from recursivenodes.nodes import _decode_family, _recursive
 
 from FIAT.orientation_utils import (
+    Orientation,
     make_cell_orientation_reflection_map_simplex,
-    make_cell_orientation_reflection_map_tensorproduct)
+    make_cell_orientation_reflection_map_tensorproduct,
+    make_entity_permutations_simplex,
+)
 
 POINT = 0
 LINE = 1
@@ -256,6 +259,52 @@ def cell_orientation_reflection_map(self):
         """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``)."""
         raise NotImplementedError("Should be implemented in a subclass.")
 
+    def extract_extrinsic_orientation(self, o):
+        """Extract extrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation.
+
+        Returns
+        -------
+        Orientation
+            Extrinsic orientation.
+
+        """
+        raise NotImplementedError("Should be implemented in a subclass.")
+
+    def extract_intrinsic_orientation(self, o, axis):
+        """Extract intrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation.
+        axis : int
+            Reference cell axis for which intrinsic orientation is computed.
+
+        Returns
+        -------
+        Orientation
+            Intrinsic orientation.
+
+        """
+        raise NotImplementedError("Should be implemented in a subclass.")
+
+    @property
+    def extrinsic_orientation_permutation_map(self):
+        """A map from extrinsic orientations to corresponding axis permutation matrices.
+
+        Notes
+        -----
+        result[eo] gives the physical axis-reference axis permutation matrix corresponding to
+        eo (extrinsic orientation).
+
+        """
+        raise NotImplementedError("Should be implemented in a subclass.")
+
     def is_simplex(self):
         return False
 
@@ -725,6 +774,58 @@ def contains_point(self, point, epsilon=0.0, entity=None):
         """
         return self.distance_to_point_l1(point, entity=entity) <= epsilon
 
+    def extract_extrinsic_orientation(self, o):
+        """Extract extrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation.
+
+        Returns
+        -------
+        Orientation
+            Extrinsic orientation.
+
+        """
+        if not isinstance(o, Orientation):
+            raise TypeError(f"Expecting an instance of Orientation : got {o}")
+        return 0
+
+    def extract_intrinsic_orientation(self, o, axis):
+        """Extract intrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation.
+        axis : int
+            Reference cell axis for which intrinsic orientation is computed.
+
+        Returns
+        -------
+        Orientation
+            Intrinsic orientation.
+
+        """
+        if not isinstance(o, Orientation):
+            raise TypeError(f"Expecting an instance of Orientation : got {o}")
+        if axis != 0:
+            raise ValueError(f"axis ({axis}) != 0")
+        return o
+
+    @property
+    def extrinsic_orientation_permutation_map(self):
+        """A map from extrinsic orientations to corresponding axis permutation matrices.
+
+        Notes
+        -----
+        result[eo] gives the physical axis-reference axis permutation matrix corresponding to
+        eo (extrinsic orientation).
+
+        """
+        return numpy.diag((1, )).astype(int).reshape((1, 1, 1))
+
 
 class Simplex(SimplicialComplex):
     r"""Abstract class for a reference simplex.
@@ -1162,6 +1263,75 @@ def __ge__(self, other):
     def __le__(self, other):
         return self.compare(operator.le, other)
 
+    def extract_extrinsic_orientation(self, o):
+        """Extract extrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation.
+
+        Returns
+        -------
+        Orientation
+            Extrinsic orientation.
+
+        Notes
+        -----
+        The difinition of orientations used here must be consistent with
+        that used in make_entity_permutations_tensorproduct.
+
+        """
+        if not isinstance(o, Orientation):
+            raise TypeError(f"Expecting an instance of Orientation : got {o}")
+        dim = len(self.cells)
+        size_io = 2  # Number of possible intrinsic orientations along each axis.
+        return o // size_io**dim
+
+    def extract_intrinsic_orientation(self, o, axis):
+        """Extract intrinsic orientation.
+
+        Parameters
+        ----------
+        o : Orientation
+            Total orientation. ``//`` and ``%`` must be overloaded in type(o).
+        axis : int
+            Reference cell axis for which intrinsic orientation is computed.
+
+        Returns
+        -------
+        Orientation
+            Intrinsic orientation.
+
+        Notes
+        -----
+        Must be consistent with make_entity_permutations_tensorproduct.
+
+        """
+        if not isinstance(o, Orientation):
+            raise TypeError(f"Expecting an instance of Orientation : got {o}")
+        dim = len(self.cells)
+        if axis >= dim:
+            raise ValueError(f"Must give 0 <= axis < {dim} : got {axis}")
+        size_io = 2  # Number of possible intrinsic orientations along each axis.
+        return o % size_io**dim // size_io**(dim - 1 - axis) % size_io
+
+    @property
+    def extrinsic_orientation_permutation_map(self):
+        """A map from extrinsic orientations to corresponding axis permutation matrices.
+
+        Notes
+        -----
+        result[eo] gives the physical axis-reference axis permutation matrix corresponding to
+        eo (extrinsic orientation).
+
+        """
+        dim = len(self.cells)
+        a = numpy.zeros((factorial(dim), dim, dim), dtype=int)
+        ai = numpy.array(list(make_entity_permutations_simplex(dim - 1, 2).values()), dtype=int).reshape((factorial(dim), dim, 1))
+        numpy.put_along_axis(a, ai, 1, axis=2)
+        return a
+
 
 class UFCQuadrilateral(Cell):
     r"""This is the reference quadrilateral with vertices