CrouzeixRaviart: integral variant (#118)

* CrouzeixRaviart: integral variant

* High-order CR variants

FIAT/crouzeix_raviart.py

@@ -9,62 +9,84 @@
 #
 # Last changed: 2010-01-28
 
+import numpy
 from FIAT import finite_element, polynomial_set, dual_set, functional
-
-
-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.check_format_variant import check_format_variant
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
 
 
 class CrouzeixRaviartDualSet(dual_set.DualSet):
-    """Dual basis for Crouzeix-Raviart element (linears continuous at
-    boundary midpoints)."""
-
-    def __init__(self, cell, degree):
 
+    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()
+
+        if degree > 1 and sd != 2:
+            raise NotImplementedError("High-order Crouzeix-Raviart is only implemented on triangles.")
 
         # Initialize empty nodes and entity_ids
-        entity_ids = _initialize_entity_ids(topology)
-        nodes = [None for i in list(topology[d - 1].keys())]
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
         # Construct nodes and entity_ids
-        for i in topology[d - 1]:
-
-            # Construct midpoint
-            x = cell.make_points(d - 1, i, d)[0]
-
-            # Degree of freedom number i is evaluation at midpoint
-            nodes[i] = functional.PointEvaluation(cell, x)
-            entity_ids[d - 1][i] += [i]
-
-        # Initialize super-class
-        super().__init__(nodes, cell, entity_ids)
+        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 + interpolant_deg-1)
+                    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)
+                    poly_set = polynomial_set.ONPolynomialSet(facet, k)
+                    Phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
+
+                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:
+            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)
+                    if dim == sd-1 and dim != 0:
+                        pts = ref_el.make_points(dim, i, degree-1, variant="gl", interior=0)
+                    else:
+                        pts = ref_el.make_points(dim, i, degree, variant="gll")
+                    nodes.extend(functional.PointEvaluation(ref_el, x) for x 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):
+    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)
-        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)

FIAT/reference_element.py

@@ -510,7 +510,7 @@ def compute_face_edge_tangents(self, dim, entity_id):
                 v1.append(dest)
         return vert_coords[v1] - vert_coords[v0]
 
-    def make_points(self, dim, entity_id, order, variant=None):
+    def make_points(self, dim, entity_id, order, variant=None, interior=1):
         """Constructs a lattice of points on the entity_id:th
         facet of dimension dim.  Order indicates how many points to
         include in each direction."""
@@ -520,7 +520,7 @@ def make_points(self, dim, entity_id, order, variant=None):
             entity_verts = \
                 self.get_vertices_of_subcomplex(
                     self.get_topology()[dim][entity_id])
-            return make_lattice(entity_verts, order, 1, variant=variant)
+            return make_lattice(entity_verts, order, interior=interior, variant=variant)
         else:
             raise ValueError("illegal dimension")
 

finat/fiat_elements.py

@@ -356,8 +356,8 @@ def __init__(self, cell, degree):
 
 
 class CrouzeixRaviart(ScalarFiatElement):
-    def __init__(self, cell, degree):
-        super().__init__(FIAT.CrouzeixRaviart(cell, degree))
+    def __init__(self, cell, degree, variant=None):
+        super().__init__(FIAT.CrouzeixRaviart(cell, degree, variant=variant))
 
 
 class Lagrange(ScalarFiatElement):

finat/ufl/elementlist.py

@@ -94,7 +94,7 @@ def show_elements():
 # Elements not in the periodic table
 # TODO: Implement generic Tear operator for elements instead of this:
 register_element("Brezzi-Douglas-Fortin-Marini", "BDFM", 1, HDiv, "contravariant Piola", (1, None), simplices[1:])
-register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, 1), simplices[1:])
+register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, None), simplices[1:])
 register_element("Discontinuous Raviart-Thomas", "DRT", 1, L2, "contravariant Piola", (1, None), simplices[1:])
 register_element("Kong-Mulder-Veldhuizen", "KMV", 0, H1, "identity", (1, None), simplices[1:])
 

gem/utils.py

@@ -1,23 +1,5 @@
 import collections
-
-
-# This is copied from PyOP2, and it is here to be available for both
-# FInAT and TSFC without depending on PyOP2.
-class cached_property(object):
-    """A read-only @property that is only evaluated once. The value is cached
-    on the object itself rather than the function or class; this should prevent
-    memory leakage."""
-    def __init__(self, fget, doc=None):
-        self.fget = fget
-        self.__doc__ = doc or fget.__doc__
-        self.__name__ = fget.__name__
-        self.__module__ = fget.__module__
-
-    def __get__(self, obj, cls):
-        if obj is None:
-            return self
-        obj.__dict__[self.__name__] = result = self.fget(obj)
-        return result
+from functools import cached_property  # noqa: F401
 
 
 def groupby(iterable, key=None):

test/FIAT/regression/test_regression.py

@@ -281,7 +281,7 @@ def create_data(family, dim, degree):
         '''Create the reference data.
         '''
         kwargs = {}
-        if family in {"Regge", "Hellan-Herrmann-Johnson"}:
+        if family in {"Crouzeix-Raviart", "Regge", "Hellan-Herrmann-Johnson"}:
             kwargs["variant"] = "point"
         # Get domain and element class
         domain = ufc_simplex(dim)

test/FIAT/unit/test_fiat.py

@@ -146,6 +146,13 @@ def __init__(self, a, b):
     "CrouzeixRaviart(I, 1)",
     "CrouzeixRaviart(T, 1)",
     "CrouzeixRaviart(S, 1)",
+    "CrouzeixRaviart(T, 3)",
+    "CrouzeixRaviart(T, 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)",