WIP

Signed-off-by: Umberto Zerbinati <[email protected]>

FIAT/crouzeix_raviart.py

@@ -48,6 +48,43 @@ def __init__(self, cell, degree):
         # Initialize super-class
         super(CrouzeixRaviartDualSet, self).__init__(nodes, cell, entity_ids)
 
+class CrouzeixRaviartThreeDualSet(dual_set.DualSet):
+    """Dual basis for Crouzeix-Raviart element (linears continuous at
+    boundary midpoints)."""
+
+    def __init__(self, cell, degree):
+
+        # Get topology dictionary
+        d = cell.get_spatial_dimension()
+        topology = cell.get_topology()
+
+        # Initialize empty nodes and entity_ids
+        entity_ids = _initialize_entity_ids(topology)
+        nodes = [None for i in range(10)]
+        # Construct nodes and entity_ids
+        dof_counter = 0
+        for i in topology[d - 1]:
+
+            # Construct midpoint
+            x = cell.make_points(d-1,i,4)
+            print(x)
+            # Degree of freedom number i is evaluation at midpoint
+            nodes[dof_counter] = functional.PointEvaluation(cell, x[0])
+            nodes[dof_counter+1] = functional.PointEvaluation(cell, x[1])
+            nodes[dof_counter+2] = functional.PointEvaluation(cell, x[2])
+            entity_ids[d - 1][i] += [dof_counter]
+            entity_ids[d - 1][i] += [dof_counter+1]
+            entity_ids[d - 1][i] += [dof_counter+2]
+            dof_counter += 3
+        x = cell.make_points(d,0,3)
+        nodes[-1] = functional.PointEvaluation(cell,x[0])
+        entity_ids[2][0] += [9]
+        print("N",nodes)
+        print("E",entity_ids)
+        # Initialize super-class
+        super(CrouzeixRaviartThreeDualSet, self).__init__(nodes, cell, entity_ids)
+
+
 
 class CrouzeixRaviart(finite_element.CiarletElement):
     """The Crouzeix-Raviart finite element:
@@ -60,11 +97,16 @@ class CrouzeixRaviart(finite_element.CiarletElement):
     def __init__(self, cell, degree):
 
         # Crouzeix Raviart is only defined for polynomial degree == 1
-        if not (degree == 1):
+        if not (degree in [1,3]):
             raise Exception("Crouzeix-Raviart only defined for degree 1")
-
         # Construct polynomial spaces, dual basis and initialize
         # FiniteElement
-        space = polynomial_set.ONPolynomialSet(cell, 1)
-        dual = CrouzeixRaviartDualSet(cell, 1)
-        super(CrouzeixRaviart, self).__init__(space, dual, 1)
+        if degree == 1:
+            space = polynomial_set.ONPolynomialSet(cell, 1)
+            dual = CrouzeixRaviartDualSet(cell, 1)
+            super(CrouzeixRaviart, self).__init__(space, dual, 1)
+        elif degree == 3:
+            space = polynomial_set.ONPolynomialSet(cell, 3)
+            dual = CrouzeixRaviartThreeDualSet(cell, 3)
+            super(CrouzeixRaviart, self).__init__(space, dual, 3)
+