High-order CR variants

FIAT/crouzeix_raviart.py

@@ -11,76 +11,77 @@
 
 import numpy
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.check_format_variant import check_format_variant
+from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
-
-
-def _initialize_entity_ids(topology):
-    entity_ids = {}
-    for (i, entity) in list(topology.items()):
-        entity_ids[i] = {}
-        for j in entity:
-            entity_ids[i][j] = []
-    return entity_ids
+from FIAT.hierarchical import make_dual_bubbles
 
 
 class CrouzeixRaviartDualSet(dual_set.DualSet):
-    """Dual basis for Crouzeix-Raviart element (linears continuous at
-    boundary midpoints)."""
-
-    def __init__(self, cell, degree, variant, interpolant_deg):
 
+    def __init__(self, ref_el, degree, variant, interpolant_deg):
         # Get topology dictionary
-        d = cell.get_spatial_dimension()
-        topology = cell.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
 
         # Initialize empty nodes and entity_ids
-        entity_ids = _initialize_entity_ids(topology)
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
         nodes = []
 
         # Construct nodes and entity_ids
-        if variant == "point":
-            for i in sorted(topology[d - 1]):
-                # Construct midpoint
-                pt, = cell.make_points(d - 1, i, d)
-                # Degree of freedom number i is evaluation at midpoint
-                nodes.append(functional.PointEvaluation(cell, pt))
-                entity_ids[d - 1][i].append(i)
+        if variant == "integral":
+            for dim in sorted(top):
+                if dim == 0 and dim != sd-1:
+                    # Skip vertex dofs
+                    continue
+                facet = ref_el.construct_subelement(dim)
+                if dim == 0:
+                    Q_facet = create_quadrature(facet, degree-1 + interpolant_deg)
+                    Phis = numpy.empty((0, len(Q_facet.pts)))
+                else:
+                    Q_facet, Phis = make_dual_bubbles(facet, degree, interpolant_deg=interpolant_deg)
+                if dim == sd-1:
+                    # Include constant moment on codim 1 facets
+                    Phis = numpy.vstack((numpy.ones((1, Phis.shape[1])), Phis))
+                for i in sorted(top[dim]):
+                    cur = len(nodes)
+                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet)
+                    scale = 1 / Q.jacobian_determinant()
+                    phis = scale * Phis
+                    nodes.extend(functional.IntegralMoment(ref_el, Q, phi) for phi in phis)
+                    entity_ids[dim][i].extend(range(cur, len(nodes)))
         else:
-            facet = cell.construct_subelement(d-1)
-            Q_facet = create_quadrature(facet, degree-1 + interpolant_deg)
-            for i in sorted(topology[d - 1]):
-                # Map quadrature
-                Q = FacetQuadratureRule(cell, d-1, i, Q_facet)
-                f = 1 / Q.jacobian_determinant()
-                f_at_qpts = numpy.full(Q.get_weights().shape, f)
-                # Degree of freedom number i is integral moment on facet
-                nodes.append(functional.IntegralMoment(cell, Q, f_at_qpts))
-                entity_ids[d - 1][i].append(i)
-
-        # Initialize super-class
-        super().__init__(nodes, cell, entity_ids)
+            for dim in sorted(top):
+                if dim == 0 and dim != sd-1:
+                    # Skip vertex dofs
+                    continue
+                for i in sorted(top[dim]):
+                    cur = len(nodes)
+                    pts = ref_el.make_points(dim, i, degree)
+                    if dim == sd-1 and dim != 0:
+                        # Include midpoint on codim 1 facets
+                        pts = [*ref_el.make_points(dim, i, dim+1), *pts]
+                    nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts)
+                    entity_ids[dim][i].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
 
 
 class CrouzeixRaviart(finite_element.CiarletElement):
     """The Crouzeix-Raviart finite element:
 
     K:                 Triangle/Tetrahedron
-    Polynomial space:  P_1
-    Dual basis:        Evaluation at facet midpoints
+    Polynomial space:  P_k
+    Dual basis:        Evaluation at points or integral moments
     """
 
-    def __init__(self, cell, degree, variant=None):
+    def __init__(self, ref_el, degree, variant=None):
 
         variant, interpolant_deg = check_format_variant(variant, degree)
 
-        # Crouzeix Raviart is only defined for polynomial degree == 1
-        if not (degree == 1):
-            raise Exception("Crouzeix-Raviart only defined for degree 1")
+        if degree % 2 != 1:
+            raise ValueError("Crouzeix-Raviart only defined for odd degree")
 
-        # Construct polynomial spaces, dual basis and initialize
-        # FiniteElement
-        space = polynomial_set.ONPolynomialSet(cell, 1)
-        dual = CrouzeixRaviartDualSet(cell, 1, variant, interpolant_deg)
-        super().__init__(space, dual, 1)
+        space = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble")
+        dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg)
+        super().__init__(space, dual, degree)

FIAT/hierarchical.py

@@ -17,12 +17,13 @@
 from FIAT.check_format_variant import parse_lagrange_variant
 
 
-def make_dual_bubbles(ref_el, degree, codim=0):
+def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None):
     """Tabulate the L2-duals of the hierarchical C0 basis."""
-    Q = create_quadrature(ref_el, 2 * degree)
+    if interpolant_deg is None:
+        interpolant_deg = degree
+    Q = create_quadrature(ref_el, degree + interpolant_deg)
     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)

test/FIAT/unit/test_fiat.py

@@ -145,7 +145,16 @@ def __init__(self, a, b):
     "DiscontinuousTaylor(S, 2)",
     "CrouzeixRaviart(I, 1)",
     "CrouzeixRaviart(T, 1)",
+    "CrouzeixRaviart(T, 3)",
+    "CrouzeixRaviart(T, 5)",
     "CrouzeixRaviart(S, 1)",
+    "CrouzeixRaviart(S, 3)",
+    "CrouzeixRaviart(S, 5)",
+    "CrouzeixRaviart(I, 1, variant='point')",
+    "CrouzeixRaviart(T, 1, variant='point')",
+    "CrouzeixRaviart(S, 1, variant='point')",
+    "CrouzeixRaviart(T, 3, variant='point')",
+    "CrouzeixRaviart(T, 5, variant='point')",
     "RaviartThomas(I, 1)",
     "RaviartThomas(I, 2)",
     "RaviartThomas(I, 3)",