Implement DG variants, remove topological DOF ordering in DG

FIAT/discontinuous_lagrange.py

@@ -7,8 +7,11 @@
 
 import itertools
 import numpy as np
-from FIAT import finite_element, polynomial_set, dual_set, functional, P0
 from FIAT.polynomial_set import mis
+from FIAT import finite_element, polynomial_set, dual_set, functional, P0
+from FIAT.orientation_utils import make_entity_permutations_simplex
+from FIAT.barycentric_interpolation import LagrangePolynomialSet
+from FIAT.reference_element import LINE, make_lattice
 
 
 def make_entity_permutations(dim, npoints):
@@ -140,45 +143,46 @@ def make_entity_permutations(dim, npoints):
 
 
 class DiscontinuousLagrangeDualSet(dual_set.DualSet):
-    """The dual basis for Lagrange elements.  This class works for
+    """The dual basis for discontinuous Lagrange elements.  This class works for
     simplices of any dimension.  Nodes are point evaluation at
-    equispaced points.  This is the discontinuous version where
+    lattice points. This is the discontinuous version where
     all nodes are topologically associated with the cell itself"""
 
-    def __init__(self, ref_el, degree):
+    def __init__(self, ref_el, degree, variant="equispaced"):
         entity_ids = {}
-        nodes = []
         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 = make_entity_permutations(dim, degree + 1 if dim == len(top) - 1 else -1)
+            perms = make_entity_permutations_simplex(dim, degree + 1 if dim == len(top)-1 else -1)
             for entity in sorted(top[dim]):
-                pts_cur = ref_el.make_points(dim, entity, degree)
-                nodes_cur = [functional.PointEvaluation(ref_el, x)
-                             for x in pts_cur]
-                nnodes_cur = len(nodes_cur)
-                nodes += nodes_cur
                 entity_ids[dim][entity] = []
-                cur += nnodes_cur
                 entity_permutations[dim][entity] = perms
-        entity_ids[dim][0] = list(range(len(nodes)))
 
+        # make nodes by getting points
+        pts = make_lattice(ref_el.get_vertices(), degree, variant=variant)
+        nodes = [functional.PointEvaluation(ref_el, x) for x in pts]
+        entity_ids[dim][0] = list(range(len(nodes)))
         super(DiscontinuousLagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
 
 
 class HigherOrderDiscontinuousLagrange(finite_element.CiarletElement):
     """The discontinuous Lagrange finite element.  It is what it is."""
 
-    def __init__(self, ref_el, degree):
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
-        dual = DiscontinuousLagrangeDualSet(ref_el, degree)
+    def __init__(self, ref_el, degree, variant="equispaced"):
+        dual = DiscontinuousLagrangeDualSet(ref_el, degree, variant=variant)
+        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
+            points = []
+            for node in dual.nodes:
+                # Assert singleton point for each node.
+                pt, = node.get_point_dict().keys()
+                points.append(pt)
+            poly_set = LagrangePolynomialSet(ref_el, points)
+        else:
+            poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
         formdegree = ref_el.get_spatial_dimension()  # n-form
         super(HigherOrderDiscontinuousLagrange, self).__init__(poly_set, dual, degree, formdegree)
 

FIAT/gauss_legendre.py

@@ -6,54 +6,12 @@
 #
 # Written by David A. Ham ([email protected]), 2015
 #
-# Modified by Pablo D. Brubeck ([email protected]), 2021
+# Modified by Pablo D. Brubeck ([email protected]), 2023
 
-from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.reference_element import POINT, LINE, TRIANGLE, TETRAHEDRON
-from FIAT.orientation_utils import make_entity_permutations_simplex
-from FIAT.barycentric_interpolation import LagrangePolynomialSet
-from FIAT.reference_element import make_lattice
+from FIAT import discontinuous_lagrange
 
 
-class GaussLegendreDualSet(dual_set.DualSet):
-    """The dual basis for discontinuous elements with nodes at the
-    (recursive) Gauss-Legendre points."""
-
-    def __init__(self, ref_el, degree):
-        entity_ids = {}
-        entity_permutations = {}
-        top = ref_el.get_topology()
-        for dim in sorted(top):
-            entity_ids[dim] = {}
-            entity_permutations[dim] = {}
-            perms = make_entity_permutations_simplex(dim, degree + 1 if dim == len(top)-1 else -1)
-            for entity in sorted(top[dim]):
-                entity_ids[dim][entity] = []
-                entity_permutations[dim][entity] = perms
-
-        # make nodes by getting points
-        pts = make_lattice(ref_el.get_vertices(), degree, variant="gl")
-        nodes = [functional.PointEvaluation(ref_el, x) for x in pts]
-        entity_ids[dim][0] = list(range(len(nodes)))
-        super(GaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
-
-
-class GaussLegendre(finite_element.CiarletElement):
+class GaussLegendre(discontinuous_lagrange.HigherOrderDiscontinuousLagrange):
     """Simplicial discontinuous element with nodes at the (recursive) Gauss-Legendre points."""
     def __init__(self, ref_el, degree):
-        if ref_el.shape not in {POINT, LINE, TRIANGLE, TETRAHEDRON}:
-            raise ValueError("Gauss-Legendre elements are only defined on simplices.")
-        dual = GaussLegendreDualSet(ref_el, degree)
-        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
-            points = []
-            for node in dual.nodes:
-                # Assert singleton point for each node.
-                pt, = node.get_point_dict().keys()
-                points.append(pt)
-            poly_set = LagrangePolynomialSet(ref_el, points)
-        else:
-            poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
-        formdegree = ref_el.get_spatial_dimension()  # n-form
-        super(GaussLegendre, self).__init__(poly_set, dual, degree, formdegree)
+        super(GaussLegendre, self).__init__(ref_el, degree, variant="gl")

FIAT/lagrange.py

@@ -14,7 +14,7 @@
 class LagrangeDualSet(dual_set.DualSet):
     """The dual basis for Lagrange elements.  This class works for
     simplices of any dimension.  Nodes are point evaluation at
-    equispaced points."""
+    lattice points."""
 
     def __init__(self, ref_el, degree, variant="equispaced"):
         entity_ids = {}
@@ -47,6 +47,8 @@ class Lagrange(finite_element.CiarletElement):
     """The Lagrange finite element.  It is what it is."""
 
     def __init__(self, ref_el, degree, variant="equispaced"):
+        if variant not in ["equispaced", "gll"]:
+            raise ValueError
         dual = LagrangeDualSet(ref_el, degree, variant=variant)
         if ref_el.shape == LINE:
             # In 1D we can use the primal basis as the expansion set,