Use ONPolynomialSet for equispaced Lagrange

FIAT/lagrange.py

@@ -7,14 +7,16 @@
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
 from FIAT.orientation_utils import make_entity_permutations_simplex
+from FIAT.barycentric_interpolation import LagrangePolynomialSet
+from FIAT.reference_element import LINE
 
 
 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."""
 
-    def __init__(self, ref_el, degree, variant=None):
+    def __init__(self, ref_el, degree, variant="equispaced"):
         entity_ids = {}
         nodes = []
         entity_permutations = {}
@@ -44,8 +46,18 @@ def __init__(self, ref_el, degree, variant=None):
 class Lagrange(finite_element.CiarletElement):
     """The Lagrange finite element.  It is what it is."""
 
-    def __init__(self, ref_el, degree):
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
-        dual = LagrangeDualSet(ref_el, degree)
+    def __init__(self, ref_el, degree, variant="equispaced"):
+        dual = LagrangeDualSet(ref_el, degree, variant=variant)
+        if ref_el.shape == LINE and variant != "equispaced":
+            # 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 = 0  # 0-form
         super(Lagrange, self).__init__(poly_set, dual, degree, formdegree)