Testing stuff to understand Crouzeix Raviart and Crouzeix Falk element.

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

FIAT/__init__.py

@@ -35,6 +35,7 @@
 from FIAT.P0 import P0
 from FIAT.raviart_thomas import RaviartThomas
 from FIAT.crouzeix_raviart import CrouzeixRaviart
+from FIAT.crouzeix_falk import CrouzeixFalk
 from FIAT.regge import Regge
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
 from FIAT.arnold_winther import ArnoldWinther

FIAT/crouzeix_falk.py

@@ -0,0 +1,21 @@
+from FIAT.finite_element import FiniteElement
+from FIAT.dual_set import DualSet
+from FIAT.polynomial_set import mis
+from FIAT.symfem import *
+from FIAT import functional,quadrature
+import symfem
+import numpy
+
+class CrouzeixFalk(SymFEM):
+    """An implementation fo the Crouzeix Falk finite element."""
+
+    def __init__(self, ref_el, degree):
+        k = 0  # 0-formi
+        topology = ref_el.get_topology()
+        print("Degree {}".format(degree))
+        self.sym_el = symfem.create_element("triangle", "CF", degree)
+
+        entity_ids = SymFEM_initialize_entity_ids(topology,self.sym_el.entity_dofs)
+        dual = SymFEMDualSet(ref_el, degree,self.sym_el)
+        self.generateBasis();
+        super(CrouzeixFalk, self).__init__(ref_el, entity_ids, dual, degree, k)

FIAT/crouzeix_raviart.py

@@ -10,7 +10,8 @@
 # Last changed: 2010-01-28
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
-
+import math
+import numpy as np
 
 def _initialize_entity_ids(topology):
     entity_ids = {}
@@ -20,6 +21,24 @@ def _initialize_entity_ids(topology):
             entity_ids[i][j] = []
     return entity_ids
 
+gauss_points =  {
+    1: {
+            0: [(1/2, 1/2)],
+            1: [(0  , 1/2)],
+            2: [(1/2  , 0)]
+        },
+    3: {
+        0: [(1-(1/2)+(1/2)*math.sqrt((3/5)),(1/2)-(1/2)*math.sqrt((3/5))),
+            ((1/2), (1/2)),
+            (1-(1/2)-(1/2)*math.sqrt((3/5)),(1/2)+(1/2)*math.sqrt((3/5)))],
+        1: [(0,(1/2)-(1/2)*math.sqrt((3/5))),
+            (0, (1/2)),
+            (0,(1/2)+(1/2)*math.sqrt((3/5)))],
+        2: [((1/2)-(1/2)*math.sqrt((3/5)),0),
+            ((1/2),0),
+            ((1/2)+(1/2)*math.sqrt((3/5)),0)]
+    }
+}
 
 class CrouzeixRaviartDualSet(dual_set.DualSet):
     """Dual basis for Crouzeix-Raviart element (linears continuous at
@@ -33,19 +52,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 = [None for _ in range(int(0.5*(degree+1)*(degree+2)))]
 
         # Construct nodes and entity_ids
-        for i in topology[d - 1]:
-
+        dof_counter = 0
+        for i in sorted(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]
-
+            #import pdb; pdb.set_trace()
+            #x = cell.make_points(d - 1, i, d)[0]
+            #print(x)
+            for pt_idx in range(len(gauss_points[degree][i])):
+                print(gauss_points[degree][i][pt_idx])
+                x = gauss_points[degree][i][pt_idx]
+                # Degree of freedom number i is evaluation at midpoint
+                nodes[dof_counter] = functional.PointEvaluation(cell, x)
+                entity_ids[d - 1][i] += [dof_counter]
+                dof_counter += 1
+        if degree == 3:
+            nodes[-1] = functional.PointEvaluation(cell, (1/3, 1/3))
+            entity_ids[d][0] += [dof_counter]
         # Initialize super-class
+        print(nodes, entity_ids)
         super(CrouzeixRaviartDualSet, self).__init__(nodes, cell, entity_ids)
 
 
@@ -60,11 +87,11 @@ 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)
+        space = polynomial_set.ONPolynomialSet(cell, degree)
+        dual = CrouzeixRaviartDualSet(cell, degree)
         super(CrouzeixRaviart, self).__init__(space, dual, 1)

FIAT/symfem.py

@@ -0,0 +1,120 @@
+from FIAT.finite_element import FiniteElement
+from FIAT.dual_set import DualSet
+from FIAT.polynomial_set import mis
+from FIAT import functional
+import symfem
+from sympy import lambdify,diff,Expr
+from sympy.abc import x,y,z
+import numpy
+
+def SymFEM_initialize_entity_ids(topology,ref_map):
+    entity_ids = {}
+    for (i, entity) in list(topology.items()):
+        entity_ids[i] = {}
+        for j in entity:
+            entity_ids[i][j] = ref_map(i,j)
+    return entity_ids
+
+def SymFEMRefToUFLRef(P):
+    return tuple(map(lambda i: 2*i-1,P))
+
+class SymFEMDualSet(DualSet):
+    """The dual basis for Bernstein elements."""
+
+    def __init__(self, ref_el, degree,sym_el):
+        topology = ref_el.get_topology()
+        entity_ids = SymFEM_initialize_entity_ids(topology,sym_el.entity_dofs)
+        nodes =  [None for i in range(len(sym_el.dofs))]
+        i=0
+        for dof in sym_el.dofs:
+            if dof.__class__ == symfem.functionals.PointEvaluation:
+                #Wrap DoF of Point Evaluation type
+                P = tuple(map(numpy.float64,dof.point));
+                nodes[i] = functional.PointEvaluation(ref_el,P)#SymFEMRefToUFLRef(P)) 
+
+            else:
+                raise RuntimeError('This type of DoF has not yet been wrapped from SymFEM.') 
+            i=i+1
+        super(SymFEMDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+class SymFEM(FiniteElement):
+    """A finite element generated from symfem."""
+
+    def __init__(self, ref_el,entity_ids,dual,degree,k,mapping="affine"):
+        self.order = degree
+        self.formdegree = k
+        self.ref_el = ref_el
+        self.dual = dual
+        self.entity_ids = entity_ids
+        self._mapping = mapping
+        self.ref_complex = ref_el
+    def generateBasis(self):
+        self.sym_basis = self.sym_el.get_basis_functions()
+        self.numberBases = len(self.sym_basis)
+        self.B = {}
+        for i in range(self.numberBases):
+            self.B[i] = {};
+            self.B[i][(0,0)] = lambdify([(x,y)],self.sym_basis[i]._sympy_())
+            self.B[i][(1,0)] = lambdify([(x,y)],diff(self.sym_basis[i]._sympy_(),x))
+            self.B[i][(0,1)] = lambdify([(x,y)],diff(self.sym_basis[i]._sympy_(),y))
+
+    def degree(self):
+        """The degree of the polynomial space."""
+        return self.get_order()
+
+    def value_shape(self):
+        """The value shape of the finite element functions."""
+        return ()
+
+    def entity_dofs(self):
+        """Return the map of topological entities to degrees of
+        freedom for the finite element."""
+        return self.entity_ids
+
+    def tabulate(self, order, points, entity=None):
+        """Return tabulated values of derivatives up to given order of
+        basis functions at given points.
+
+        :arg order: The maximum order of derivative.
+        :arg points: An iterable of points.
+        :arg entity: Optional (dimension, entity number) pair
+                     indicating which topological entity of the
+                     reference element to tabulate on.  If ``None``,
+                     default cell-wise tabulation is performed.
+        """
+
+        if order > 1:
+            raise RuntimeError("SymFEM are only implemented for first derivative operator.");
+        # Transform points to reference cell coordinates
+        ref_el = self.get_reference_element()
+        if entity is None:
+            entity = (ref_el.get_spatial_dimension(), 0)
+
+        entity_dim, entity_id = entity
+        entity_transform = ref_el.get_entity_transform(entity_dim, entity_id)
+        cell_points = list(map(entity_transform, points))
+
+        # Evaluate everything
+        deg = self.degree()
+        dim = ref_el.get_spatial_dimension()
+
+        # Rearrange result
+        space_dim = self.space_dimension()
+        dtype = numpy.float64#numpy.array(list(raw_result.values())).dtype
+        result = {alpha: numpy.zeros((space_dim, len(cell_points)), dtype=dtype)
+                  for o in range(order + 1)
+                  for alpha in mis(dim, o)}
+        for i in range(self.numberBases):
+            for o in range(order+1):
+                for alpha, vec in self.basisFunction(cell_points, i, o).items():
+                    result[alpha][i, :] = vec
+        return result
+
+    def basisFunction(self,points,i, order):
+        points = numpy.asarray(points);
+        N, d_1 = points.shape
+        result = {}
+        for alpha in mis(d_1, order):
+            values = self.B[i][alpha](points.T);
+            result[alpha] = values
+        return result