From 10a4f98195574cc43d0e57e8695dbce0f3f894af Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Wed, 11 Dec 2024 21:35:30 +0000
Subject: [PATCH] CrouzeixRaviart: integral variant

---
 FIAT/crouzeix_raviart.py                | 40 +++++++++++++++++--------
 gem/utils.py                            | 20 +------------
 test/FIAT/regression/test_regression.py |  2 +-
 3 files changed, 30 insertions(+), 32 deletions(-)

diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py
index b809f7a87..149cfa03a 100644
--- a/FIAT/crouzeix_raviart.py
+++ b/FIAT/crouzeix_raviart.py
@@ -9,7 +9,11 @@
 #
 # Last changed: 2010-01-28
 
+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.quadrature import FacetQuadratureRule
 
 
 def _initialize_entity_ids(topology):
@@ -25,7 +29,7 @@ class CrouzeixRaviartDualSet(dual_set.DualSet):
     """Dual basis for Crouzeix-Raviart element (linears continuous at
     boundary midpoints)."""
 
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree, variant, interpolant_deg):
 
         # Get topology dictionary
         d = cell.get_spatial_dimension()
@@ -33,17 +37,27 @@ def __init__(self, cell, degree):
 
         # Initialize empty nodes and entity_ids
         entity_ids = _initialize_entity_ids(topology)
-        nodes = [None for i in list(topology[d - 1].keys())]
+        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]
+        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)
+        else:
+            facet = cell.construct_subelement(d-1)
+            Q_facet = create_quadrature(facet, degree + 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)
@@ -57,7 +71,9 @@ class CrouzeixRaviart(finite_element.CiarletElement):
     Dual basis:        Evaluation at facet midpoints
     """
 
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree, variant=None):
+
+        variant, interpolant_deg = check_format_variant(variant, degree)
 
         # Crouzeix Raviart is only defined for polynomial degree == 1
         if not (degree == 1):
@@ -66,5 +82,5 @@ def __init__(self, cell, degree):
         # Construct polynomial spaces, dual basis and initialize
         # FiniteElement
         space = polynomial_set.ONPolynomialSet(cell, 1)
-        dual = CrouzeixRaviartDualSet(cell, 1)
+        dual = CrouzeixRaviartDualSet(cell, 1, variant, interpolant_deg)
         super().__init__(space, dual, 1)
diff --git a/gem/utils.py b/gem/utils.py
index 12e0e0f6c..855859f55 100644
--- a/gem/utils.py
+++ b/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
 
 
 def groupby(iterable, key=None):
diff --git a/test/FIAT/regression/test_regression.py b/test/FIAT/regression/test_regression.py
index a4101d0cb..fc5ee01e5 100644
--- a/test/FIAT/regression/test_regression.py
+++ b/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)