Remove entity permutations for non-Lagrange elements

FIAT/fdm_element.py

@@ -11,10 +11,9 @@
 
 from FIAT import dual_set, finite_element, functional, quadrature
 from FIAT.barycentric_interpolation import LagrangePolynomialSet
-from FIAT.hierarchical import IntegratedLegendre
-from FIAT.orientation_utils import make_entity_permutations_simplex
-from FIAT.P0 import P0Dual
+from FIAT.polynomial_set import ONPolynomialSet
 from FIAT.reference_element import LINE
+from FIAT.P0 import P0
 
 
 def sym_eig(A, B):
@@ -40,23 +39,24 @@ def tridiag_eig(A, B):
     numpy.reciprocal(Z, out=Z)
     numpy.multiply(numpy.sqrt(Z), V, out=V)
     numpy.multiply(V, a[:, None], out=V)
-    return Z, V
+    # Reorder by increasing eigenvalue
+    return Z[::-1], V[:, ::-1]
 
 
 class FDMDual(dual_set.DualSet):
     """The dual basis for 1D elements with FDM shape functions."""
     def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False):
         # Define the generalized eigenproblem on a reference element
-        embedded_degree = degree + formdegree
-        embedded = IntegratedLegendre(ref_el, embedded_degree)
-        edim = embedded.space_dimension()
-        self.embedded = embedded
+        P = ONPolynomialSet(ref_el, degree + formdegree, variant="bubble")
+        Pdim = len(P)
+        # Apply even / odd reordering on edge bubbles
+        P = P.take([*range(2), *range(2, Pdim, 2), *range(3, Pdim, 2)])
+        self.poly_set = P
 
-        vertices = ref_el.get_vertices()
         if bc_order == 1 and formdegree == 0:
-            rule = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, edim+1)
+            rule = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, Pdim+1)
         else:
-            rule = quadrature.GaussLegendreQuadratureLineRule(ref_el, edim)
+            rule = quadrature.GaussLegendreQuadratureLineRule(ref_el, Pdim)
         self.rule = rule
 
         solve_eig = sym_eig
@@ -65,21 +65,21 @@ def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False
 
         # Tabulate the BC nodes
         if bc_order == 0:
-            C = numpy.empty((0, edim), "d")
+            C = numpy.empty((0, Pdim), "d")
         else:
-            constraints = embedded.tabulate(bc_order-1, vertices)
+            constraints = P.tabulate(ref_el.get_vertices(), bc_order-1)
             C = numpy.transpose(numpy.column_stack(list(constraints.values())))
         bdof = slice(None, C.shape[0])
         idof = slice(C.shape[0], None)
 
-        # Tabulate the basis that splits the DOFs into interior and bcs
-        E = numpy.eye(edim)
+        # Coefficients of the vertex and interior modes
+        E = numpy.eye(Pdim)
         E[bdof, idof] = -C[:, idof]
         E[bdof, :] = numpy.linalg.solve(C[:, bdof], E[bdof, :])
 
         # Assemble the constrained Galerkin matrices on the reference cell
         k = max(1, bc_order)
-        phi = embedded.tabulate(k, rule.get_points())
+        phi = P.tabulate(rule.get_points(), k)
         wts = rule.get_weights()
         E0 = numpy.dot(E.T, phi[(0, )])
         Ek = numpy.dot(E.T, phi[(k, )])
@@ -113,34 +113,33 @@ def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False
             numpy.multiply(S, lam, out=S)
             basis = numpy.dot(S.T, Ek)
 
-        bc_nodes = []
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
         if formdegree == 0:
             if orthogonalize:
                 idof = slice(None)
             elif bc_order > 0:
-                bc_nodes += [functional.PointEvaluation(ref_el, x) for x in vertices]
-                for alpha in range(1, bc_order):
-                    bc_nodes += [functional.PointDerivative(ref_el, x, [alpha]) for x in vertices]
+                # Vertex dofs -- jet evaluation
+                for v in sorted(top[0]):
+                    cur = len(nodes)
+                    x, = ref_el.make_points(0, v, 0)
+                    nodes.append(functional.PointEvaluation(ref_el, x))
+                    nodes.extend(functional.PointDerivative(ref_el, x, (alpha, ))
+                                 for alpha in range(1, bc_order))
+                    entity_ids[0][v].extend(range(cur, len(nodes)))
 
         elif bc_order > 0:
-            basis[bdof, :] = numpy.sqrt(1.0E0 / ref_el.volume())
+            basis[bdof] = numpy.sqrt(1.0E0 / ref_el.volume())
             idof = slice(formdegree, None)
 
-        nodes = bc_nodes + [functional.IntegralMoment(ref_el, rule, f) for f in basis[idof]]
+        # Interior dofs -- moments against eigenfunctions
+        cur = len(nodes)
+        nodes.extend(functional.IntegralMoment(ref_el, rule, f) for f in basis[idof])
+        entity_ids[sd][0].extend(range(cur, len(nodes)))
 
-        if len(bc_nodes) > 0:
-            entity_ids = {0: {0: [0], 1: [1]},
-                          1: {0: list(range(2, degree+1))}}
-            entity_permutations = {}
-            entity_permutations[0] = {0: {0: [0]}, 1: {0: [0]}}
-            entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree - 1)}
-        else:
-            entity_ids = {0: {0: [], 1: []},
-                          1: {0: list(range(0, degree+1))}}
-            entity_permutations = {}
-            entity_permutations[0] = {0: {0: []}, 1: {0: []}}
-            entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)}
-        super().__init__(nodes, ref_el, entity_ids, entity_permutations)
+        super().__init__(nodes, ref_el, entity_ids)
 
 
 class FDMFiniteElement(finite_element.CiarletElement):
@@ -158,16 +157,18 @@ def _bc_order(self):
     def _formdegree(self):
         pass
 
+    def __new__(cls, ref_el, degree):
+        if cls._formdegree == 1 and degree == 0:
+            return P0(ref_el)
+        return super().__new__(cls)
+
     def __init__(self, ref_el, degree):
         if ref_el.shape != LINE:
             raise ValueError("%s is only defined in one dimension." % type(self))
-        if degree == 0:
-            dual = P0Dual(ref_el)
-        else:
-            dual = FDMDual(ref_el, degree, bc_order=self._bc_order,
-                           formdegree=self._formdegree, orthogonalize=self._orthogonalize)
+        dual = FDMDual(ref_el, degree, bc_order=self._bc_order,
+                       formdegree=self._formdegree, orthogonalize=self._orthogonalize)
         if self._formdegree == 0:
-            poly_set = dual.embedded.poly_set
+            poly_set = dual.poly_set
         else:
             lr = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1)
             poly_set = LagrangePolynomialSet(ref_el, lr.get_points())

FIAT/hierarchical.py

@@ -8,17 +8,19 @@
 
 import numpy
 
-from FIAT import finite_element, dual_set, functional, P0
+from FIAT import finite_element, dual_set, functional
 from FIAT.reference_element import symmetric_simplex
-from FIAT.orientation_utils import make_entity_permutations_simplex
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.polynomial_set import ONPolynomialSet, make_bubbles
-from FIAT.check_format_variant import parse_lagrange_variant
+from FIAT.check_format_variant import check_format_variant
+from FIAT.P0 import P0
 
 
 def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None):
     """Tabulate the L2-duals of the hierarchical C0 basis."""
+    if ref_el.get_spatial_dimension() == 0:
+        degree = 0
     if interpolant_deg is None:
         interpolant_deg = degree
     Q = create_quadrature(ref_el, degree + interpolant_deg)
@@ -32,105 +34,86 @@ def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None):
 
 class LegendreDual(dual_set.DualSet):
     """The dual basis for Legendre elements."""
-    def __init__(self, ref_el, degree, codim=0):
-        nodes = []
-        entity_ids = {}
-        entity_permutations = {}
-
+    def __init__(self, ref_el, degree, codim=0, interpolant_deg=None):
+        if interpolant_deg is None:
+            interpolant_deg = degree
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
-        for dim in sorted(top):
-            npoints = degree + 1 if dim == sd - codim else 0
-            perms = make_entity_permutations_simplex(dim, npoints)
-            entity_permutations[dim] = {}
-            entity_ids[dim] = {}
-            if npoints == 0:
-                for entity in sorted(top[dim]):
-                    entity_ids[dim][entity] = []
-                    entity_permutations[dim][entity] = perms
-                continue
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
-            ref_facet = ref_el.construct_subelement(dim)
-            poly_set = ONPolynomialSet(ref_facet, degree)
-            Q_ref = create_quadrature(ref_facet, 2 * degree)
-            phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
-            for entity in sorted(top[dim]):
-                cur = len(nodes)
-                Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
-                entity_ids[dim][entity] = list(range(cur, len(nodes)))
-                entity_permutations[dim][entity] = perms
+        dim = sd - codim
+        ref_facet = ref_el.construct_subelement(dim)
+        poly_set = ONPolynomialSet(ref_facet, degree, scale="L2 piola")
+        Q_ref = create_quadrature(ref_facet, degree + interpolant_deg)
+        Phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
+        for entity in sorted(top[dim]):
+            cur = len(nodes)
+            Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
+            # phis must transform like a d-form to undo the measure transformation
+            scale = 1 / Q_facet.jacobian_determinant()
+            phis = scale * Phis
+            nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
+            entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
-        super().__init__(nodes, ref_el, entity_ids, entity_permutations)
+        super().__init__(nodes, ref_el, entity_ids)
 
 
 class Legendre(finite_element.CiarletElement):
     """Simplicial discontinuous element with Legendre polynomials."""
     def __new__(cls, ref_el, degree, variant=None):
         if degree == 0:
-            splitting, _ = parse_lagrange_variant(variant, integral=True)
-            if splitting is None:
+            splitting, variant, interpolant_deg = check_format_variant(variant, degree)
+            if splitting is None and interpolant_deg == 0:
                 # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
-                return P0.P0(ref_el)
+                return P0(ref_el)
         return super().__new__(cls)
 
     def __init__(self, ref_el, degree, variant=None):
-        splitting, _ = parse_lagrange_variant(variant, integral=True)
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
         poly_set = ONPolynomialSet(ref_el, degree)
-        dual = LegendreDual(ref_el, degree)
+        dual = LegendreDual(ref_el, degree, interpolant_deg=interpolant_deg)
         formdegree = ref_el.get_spatial_dimension()  # n-form
         super().__init__(poly_set, dual, degree, formdegree)
 
 
 class IntegratedLegendreDual(dual_set.DualSet):
     """The dual basis for integrated Legendre elements."""
-    def __init__(self, ref_el, degree):
+    def __init__(self, ref_el, degree, interpolant_deg=None):
+        if interpolant_deg is None:
+            interpolant_deg = degree
+        top = ref_el.get_topology()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
         nodes = []
-        entity_ids = {}
-        entity_permutations = {}
 
-        top = ref_el.get_topology()
         for dim in sorted(top):
-            perms = make_entity_permutations_simplex(dim, degree - dim)
-            entity_ids[dim] = {}
-            entity_permutations[dim] = {}
-            if dim == 0 or degree <= dim:
-                for entity in sorted(top[dim]):
-                    cur = len(nodes)
-                    pts = ref_el.make_points(dim, entity, degree)
-                    nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts)
-                    entity_ids[dim][entity] = list(range(cur, len(nodes)))
-                    entity_permutations[dim][entity] = perms
+            if degree <= dim:
                 continue
-
             ref_facet = symmetric_simplex(dim)
-            Q_ref, phis = make_dual_bubbles(ref_facet, degree)
+            Q_ref, Phis = make_dual_bubbles(ref_facet, degree, interpolant_deg=interpolant_deg)
             for entity in sorted(top[dim]):
                 cur = len(nodes)
                 Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-
                 # phis must transform like a d-form to undo the measure transformation
                 scale = 1 / Q_facet.jacobian_determinant()
-                Jphis = scale * phis
-
-                nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in Jphis)
-                entity_ids[dim][entity] = list(range(cur, len(nodes)))
-                entity_permutations[dim][entity] = perms
+                phis = scale * Phis
+                nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
+                entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
-        super().__init__(nodes, ref_el, entity_ids, entity_permutations)
+        super().__init__(nodes, ref_el, entity_ids)
 
 
 class IntegratedLegendre(finite_element.CiarletElement):
     """Simplicial continuous element with integrated Legendre polynomials."""
     def __init__(self, ref_el, degree, variant=None):
-        splitting, _ = parse_lagrange_variant(variant, integral=True)
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
         if degree < 1:
             raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
         poly_set = ONPolynomialSet(ref_el, degree, variant="bubble")
-        dual = IntegratedLegendreDual(ref_el, degree)
+        dual = IntegratedLegendreDual(ref_el, degree, interpolant_deg=interpolant_deg)
         formdegree = 0  # 0-form
         super().__init__(poly_set, dual, degree, formdegree)

FIAT/polynomial_set.py

@@ -291,21 +291,15 @@ def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"):
     """Construct a polynomial set with codim bubbles up to the given degree.
     """
     poly_set = ONPolynomialSet(ref_el, degree, shape=shape, scale=scale, variant="bubble")
+    if ref_el.get_spatial_dimension() == 0:
+        return poly_set
+
     entity_ids = expansions.polynomial_entity_ids(ref_el, degree, continuity="C0")
     sd = ref_el.get_spatial_dimension()
     dim = sd - codim
-    if dim == 1:
-        # Apply even / odd reordering on edge bubbles
-        indices = []
-        for entity in entity_ids[dim]:
-            ids = entity_ids[dim][entity]
-            indices.extend(ids[::2])
-            indices.extend(ids[1::2])
-    else:
-        indices = list(chain(*entity_ids[dim].values()))
-
+    indices = list(chain(*entity_ids[dim].values()))
     if shape != ():
-        ncomp = numpy.prod(shape)
+        ncomp = numpy.prod(shape, dtype=int)
         dimPk = poly_set.get_num_members() // ncomp
         indices = list((numpy.array(indices)[:, None] + dimPk * numpy.arange(ncomp)[None, :]).flat)
     poly_set = poly_set.take(indices)

FIAT/restricted.py

@@ -24,6 +24,7 @@ def __init__(self, dual, indices):
                                          for dof in dofs if dof in indices]
         nodes = [nodes_old[i] for i in indices]
         self._dual = dual
+
         super().__init__(nodes, ref_el, entity_ids)
 
     def get_indices(self, restriction_domain, take_closure=True):
@@ -69,4 +70,4 @@ def __init__(self, element, indices=None, restriction_domain=None, take_closure=
         assert all(e_mapping == mapping_new[0] for e_mapping in mapping_new)
 
         # Call constructor of CiarletElement
-        super().__init__(poly_set, dual, 0, element.get_formdegree(), mapping_new[0])
+        super().__init__(poly_set, dual, element.degree(), element.get_formdegree(), mapping_new[0])

finat/restricted.py

@@ -32,11 +32,19 @@ def restrict(element, domain, take_closure):
 @restrict.register(FiatElement)
 def restrict_fiat(element, domain, take_closure):
     try:
-        return FiatElement(FIAT.RestrictedElement(element._element,
-                           restriction_domain=domain, take_closure=take_closure))
+        re = FIAT.RestrictedElement(element._element,
+                                    restriction_domain=domain,
+                                    take_closure=take_closure)
     except ValueError:
         return null_element
 
+    if element.space_dimension() == re.space_dimension():
+        # FIAT.RestrictedElement wipes out entity_permuations.
+        # In case the restriction is trivial we return the original element
+        # to avoid reconstructing the space with an undesired permutation.
+        return element
+    return FiatElement(re)
+
 
 @restrict.register(PhysicallyMappedElement)
 def restrict_physically_mapped(element, domain, take_closure):