diff --git a/FIAT/orientation_utils.py b/FIAT/orientation_utils.py index 92c678991..2852f33f4 100644 --- a/FIAT/orientation_utils.py +++ b/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 diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 61702c54c..51e4f45e2 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -8,10 +8,12 @@ # Modified by David A. Ham (david.ham@imperial.ac.uk), 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): diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index f09751578..95f0d8f8a 100644 --- a/FIAT/reference_element.py +++ b/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