Remove repeated points from dual_basis

Author: pbrubeck

FIAT/check_format_variant.py

@@ -103,6 +103,11 @@ def parse_quadrature_scheme(ref_el, degree, quad_scheme=None):
         if opt in supported_splits:
             splitting = supported_splits[opt]
             ref_el = splitting(ref_el)
+        elif opt.startswith("KMV") and opt != "KMV":
+            match = re.match(r"^KMV(?:\((\d+)\))?$", opt)
+            degree, = match.groups()
+            degree = int(degree)
+            scheme = "KMV"
         else:
             scheme = opt
     return create_quadrature(ref_el, degree, scheme or "default")

FIAT/macro.py

@@ -406,6 +406,24 @@ def __init__(self, ref_el, Q_ref, parent_facets=None):
             Q_cur = FacetQuadratureRule(ref_el, parent_dim, entity, Q_ref)
             pts.extend(Q_cur.pts)
             wts.extend(Q_cur.wts)
+
+        # Merge points on facets, if any
+        rtol = 1E-12
+        bary = ref_el.compute_barycentric_coordinates(pts)
+        if numpy.isclose(bary, 0, rtol=rtol).any():
+            iorder = numpy.lexsort(bary.T)
+            iprev = iorder[0]
+            unique_pts = [pts[iprev]]
+            unique_wts = [wts[iprev]]
+            for icur in iorder[1:]:
+                if numpy.allclose(bary[icur], bary[iprev], rtol=rtol):
+                    unique_wts[-1] += wts[icur]
+                else:
+                    unique_pts.append(pts[icur])
+                    unique_wts.append(wts[icur])
+                iprev = icur
+            pts = unique_pts
+            wts = unique_wts
         pts = tuple(pts)
         wts = tuple(wts)
         super().__init__(ref_el, pts, wts)

FIAT/quadrature_schemes.py

@@ -34,7 +34,9 @@
 
 import numpy
 
-from FIAT.quadrature import (QuadratureRule, FacetQuadratureRule, make_quadrature,
+from FIAT.quadrature import (FacetQuadratureRule,
+                             GaussLobattoLegendreQuadratureLineRule,
+                             QuadratureRule, make_quadrature,
                              make_tensor_product_quadrature, map_quadrature)
 from FIAT.reference_element import (HEXAHEDRON, QUADRILATERAL, TENSORPRODUCT,
                                     TETRAHEDRON, TRIANGLE, symmetric_simplex, ufc_simplex)
@@ -118,7 +120,10 @@ def _kmv_lump_scheme(ref_el, degree):
     """Specialized quadrature schemes for P < 6 for KMV simplical elements."""
 
     sd = ref_el.get_spatial_dimension()
-    if sd not in {2, 3}:
+    if sd == 1:
+        num_points = degree + 1
+        return GaussLobattoLegendreQuadratureLineRule(ref_el, num_points)
+    elif sd > 3:
         raise ValueError("Dimension not supported")
 
     T = ufc_simplex(sd)

finat/fiat_elements.py

@@ -213,7 +213,18 @@ def _dual_basis(self):
             Qdense = np.zeros(Qshape, dtype=np.float64)
             for idx, value in Q.items():
                 Qdense[idx] = value
+            # Compress repeated points
+            repeated_pts = tuple(map(tuple, np.round(allpts, decimals=12)))
+            unique_pts = list(dict.fromkeys(repeated_pts))
+            if len(unique_pts) < len(repeated_pts):
+                Qrepeated = Qdense
+                Qshape = (Qshape[0], len(unique_pts), *Qshape[2:])
+                Qdense = np.zeros(Qshape, dtype=np.float64)
+                for j, i in enumerate(map(unique_pts.index, repeated_pts)):
+                    Qdense[:, i, ...] += Qrepeated[:, j, ...]
+                allpts = unique_pts
             Q = gem.Literal(Qdense)
+
         return Q, np.asarray(allpts)
 
     @property