CrouzeixRaviart: support quad_scheme

Author: pbrubeck

FIAT/check_format_variant.py

@@ -39,6 +39,8 @@ def check_format_variant(variant, degree):
         extra_degree, = match.groups()
         extra_degree = int(extra_degree) if extra_degree is not None else 0
         interpolant_degree = degree + extra_degree
+        if interpolant_degree < degree:
+            raise ValueError(f"Quadrature degree should be at least {degree}")
 
     if variant not in {"point", "integral"}:
         raise ValueError('Choose either variant="point" or variant="integral"'

FIAT/crouzeix_raviart.py

@@ -11,14 +11,13 @@
 
 import numpy
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.check_format_variant import check_format_variant
-from FIAT.quadrature_schemes import create_quadrature
+from FIAT.check_format_variant import check_format_variant, parse_quadrature_scheme
 from FIAT.quadrature import FacetQuadratureRule
 
 
 class CrouzeixRaviartDualSet(dual_set.DualSet):
 
-    def __init__(self, ref_el, degree, variant, interpolant_deg):
+    def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
         # Get topology dictionary
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
@@ -38,13 +37,14 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                     continue
                 facet = ref_el.construct_subelement(dim)
                 if dim == 0:
-                    Q_facet = create_quadrature(facet, degree + interpolant_deg-1)
+                    Q_facet = parse_quadrature_scheme(facet, degree + interpolant_deg-1, quad_scheme)
                     Phis = numpy.ones((1, len(Q_facet.pts)))
                 else:
                     k = degree - 1 if dim == sd-1 else degree - (1+dim)
                     if k < 0:
                         continue
-                    Q_facet = create_quadrature(facet, k + interpolant_deg)
+                    Q_facet = parse_quadrature_scheme(facet, k + interpolant_deg, quad_scheme)
+
                     poly_set = polynomial_set.ONPolynomialSet(facet, k)
                     Phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
 
@@ -80,7 +80,7 @@ class CrouzeixRaviart(finite_element.CiarletElement):
     Dual basis:        Evaluation at points or integral moments
     """
 
-    def __init__(self, ref_el, degree, variant=None):
+    def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
         if degree % 2 != 1:
             raise ValueError("Crouzeix-Raviart only defined for odd degree")
 
@@ -89,5 +89,5 @@ def __init__(self, ref_el, degree, variant=None):
             ref_el = splitting(ref_el)
 
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
-        dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg)
+        dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg, quad_scheme)
         super().__init__(poly_set, dual, degree)

FIAT/hierarchical.py

@@ -25,18 +25,12 @@ def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme
     if interpolant_deg is None:
         interpolant_deg = degree
 
-    if interpolant_deg >= degree:
-        Q = parse_quadrature_scheme(ref_el, degree + interpolant_deg, quad_scheme)
-        B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal")
-        P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points())
-        M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T)
-        phis = numpy.linalg.solve(M, P_at_qpts)
-        phis = numpy.dot(B.get_coeffs(), phis)
-    else:
-        k = max(degree - dim - 1, 0)
-        Q = parse_quadrature_scheme(ref_el, k + degree, quad_scheme)
-        P = ONPolynomialSet(ref_el, k, scale="L2 piola")
-        phis = P.tabulate(Q.get_points())[(0,) * dim]
+    Q = parse_quadrature_scheme(ref_el, degree + interpolant_deg, quad_scheme)
+    B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal")
+    P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points())
+    M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T)
+    phis = numpy.linalg.solve(M, P_at_qpts)
+    phis = numpy.dot(B.get_coeffs(), phis)
     return Q, phis
 
 

finat/fiat_elements.py

@@ -296,23 +296,23 @@ def point_evaluation(fiat_element, order, refcoords, entity):
 
 
 class Regge(FiatElement):  # symmetric matrix valued
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.Regge(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.Regge(cell, degree, **kwargs))
 
 
 class HellanHerrmannJohnson(FiatElement):  # symmetric matrix valued
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.HellanHerrmannJohnson(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.HellanHerrmannJohnson(cell, degree, **kwargs))
 
 
 class GopalakrishnanLedererSchoberlFirstKind(FiatElement):  # traceless matrix valued
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.GopalakrishnanLedererSchoberlFirstKind(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.GopalakrishnanLedererSchoberlFirstKind(cell, degree, **kwargs))
 
 
 class GopalakrishnanLedererSchoberlSecondKind(FiatElement):  # traceless matrix valued
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.GopalakrishnanLedererSchoberlSecondKind(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.GopalakrishnanLedererSchoberlSecondKind(cell, degree, **kwargs))
 
 
 class ScalarFiatElement(FiatElement):
@@ -328,28 +328,28 @@ def __init__(self, cell, degree):
 
 
 class Bubble(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.Bubble(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.Bubble(cell, degree, **kwargs))
 
 
 class FacetBubble(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.FacetBubble(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.FacetBubble(cell, degree, **kwargs))
 
 
 class CrouzeixRaviart(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.CrouzeixRaviart(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.CrouzeixRaviart(cell, degree, **kwargs))
 
 
 class Lagrange(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.Lagrange(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.Lagrange(cell, degree, *kwargs))
 
 
 class DiscontinuousLagrange(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.DiscontinuousLagrange(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.DiscontinuousLagrange(cell, degree, **kwargs))
 
 
 class Histopolation(ScalarFiatElement):
@@ -377,8 +377,8 @@ def __init__(self, cell, degree):
 
 
 class HDivTrace(ScalarFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.HDivTrace(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.HDivTrace(cell, degree, **kwargs))
 
 
 class VectorFiatElement(FiatElement):
@@ -388,8 +388,8 @@ def value_shape(self):
 
 
 class RaviartThomas(VectorFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.RaviartThomas(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.RaviartThomas(cell, degree, **kwargs))
 
 
 class TrimmedSerendipityFace(VectorFiatElement):
@@ -429,8 +429,8 @@ def entity_permutations(self):
 
 
 class BrezziDouglasMarini(VectorFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.BrezziDouglasMarini(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.BrezziDouglasMarini(cell, degree, *kwargs))
 
 
 class BrezziDouglasMariniCubeEdge(VectorFiatElement):
@@ -457,10 +457,10 @@ def __init__(self, cell, degree):
 
 
 class Nedelec(VectorFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.Nedelec(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.Nedelec(cell, degree, **kwargs))
 
 
 class NedelecSecondKind(VectorFiatElement):
-    def __init__(self, cell, degree, variant=None):
-        super().__init__(FIAT.NedelecSecondKind(cell, degree, variant=variant))
+    def __init__(self, cell, degree, **kwargs):
+        super().__init__(FIAT.NedelecSecondKind(cell, degree, **kwargs))