Merge pull request #75 from firedrakeproject/pbrubeck/fix/gll-ordering

GLL: ensure natural ordering of the DOFs

Author: rckirby

FIAT/gauss_lobatto_legendre.py

@@ -14,4 +14,4 @@
 class GaussLobattoLegendre(lagrange.Lagrange):
     """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto-Legendre points."""
     def __init__(self, ref_el, degree):
-        super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll")
+        super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", sort_entities=True)

FIAT/lagrange.py

@@ -14,33 +14,42 @@
 
 class LagrangeDualSet(dual_set.DualSet):
     """The dual basis for Lagrange elements.  This class works for
-    simplices of any dimension.  Nodes are point evaluation at
+    simplicial complexes of any dimension.  Nodes are point evaluation at
     recursively-defined points.
+
+    :arg ref_el: The simplicial complex.
+    :arg degree: The polynomial degree.
+    :arg point_variant: The point distribution variant passed on to recursivenodes.
+    :arg sort_entities: A flag to sort entities by support vertex ids.
+                        If false then entities are sorted first by dimension and then by
+                        entity id. The DOFs are always sorted by the entity ordering
+                        and then lexicographically by lattice multiindex.
     """
-    def __init__(self, ref_el, degree, point_variant="equispaced"):
-        entity_ids = {}
+    def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=False):
         nodes = []
+        entity_ids = {}
         entity_permutations = {}
-
-        # make nodes by getting points
-        # need to do this dimension-by-dimension, facet-by-facet
         top = ref_el.get_topology()
-
-        cur = 0
         for dim in sorted(top):
             entity_ids[dim] = {}
             entity_permutations[dim] = {}
             perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim)
             for entity in sorted(top[dim]):
-                pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant)
-                nodes_cur = [functional.PointEvaluation(ref_el, x)
-                             for x in pts_cur]
-                nnodes_cur = len(nodes_cur)
-                nodes += nodes_cur
-                entity_ids[dim][entity] = list(range(cur, cur + nnodes_cur))
-                cur += nnodes_cur
                 entity_permutations[dim][entity] = perms
 
+        entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])]
+        if sort_entities:
+            # sort the entities by support vertex ids
+            support = [top[dim][entity] for dim, entity in entities]
+            entities = [entity for verts, entity in sorted(zip(support, entities))]
+
+        # make nodes by getting points
+        # need to do this entity-by-entity
+        for dim, entity in entities:
+            cur = len(nodes)
+            pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant)
+            nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur)
+            entity_ids[dim][entity] = list(range(cur, len(nodes)))
         super(LagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
 
 
@@ -58,12 +67,16 @@ class Lagrange(finite_element.CiarletElement):
                            variant='equispaced,Iso(2)' with degree=1 gives the P2:P1 iso element.
                            variant='Alfeld' can be used to obtain a barycentrically refined
                               macroelement for Scott-Vogelius.
+    :arg sort_entities: A flag to sort entities by support vertex ids.
+                        If false then entities are sorted first by dimension and then by
+                        entity id. The DOFs are always sorted by the entity ordering
+                        and then lexicographically by lattice multiindex.
     """
-    def __init__(self, ref_el, degree, variant="equispaced"):
+    def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False):
         splitting, point_variant = parse_lagrange_variant(variant)
         if splitting is not None:
             ref_el = splitting(ref_el)
-        dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant)
+        dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, sort_entities=sort_entities)
         if ref_el.shape == LINE:
             # In 1D we can use the primal basis as the expansion set,
             # avoiding any round-off coming from a basis transformation

test/unit/test_fiat.py

@@ -394,8 +394,8 @@ def test_empty_bubble():
 
 @pytest.mark.parametrize('elements', [
     (Lagrange(I, 2), Lagrange(I, 1), Bubble(I, 2)),
-    (GaussLobattoLegendre(I, 3), Lagrange(I, 1),
-     RestrictedElement(GaussLobattoLegendre(I, 3), restriction_domain="interior")),
+    (Lagrange(I, 3, variant="gll"), Lagrange(I, 1),
+     RestrictedElement(Lagrange(I, 3, variant="gll"), restriction_domain="interior")),
     (RaviartThomas(T, 2),
      RestrictedElement(RaviartThomas(T, 2), restriction_domain='facet'),
      RestrictedElement(RaviartThomas(T, 2), restriction_domain='interior')),

test/unit/test_gauss_lobatto_legendre.py

@@ -59,8 +59,7 @@ def test_edge_dofs(dim, degree):
     # Test that edge DOFs are located at the GLL quadrature points
     line = s if dim == 1 else s.construct_subelement(1)
     lr = quadrature.GaussLobattoLegendreQuadratureLineRule(line, degree + 1)
-    # Edge DOFs are ordered with the two vertex DOFs followed by the interior DOFs
-    quadrature_points = lr.pts[::degree] + lr.pts[1:-1]
+    quadrature_points = lr.get_points()
 
     entity_dofs = fe.entity_closure_dofs()
     edge_dofs = entity_dofs[1]