diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 3583fac18..a6935548c 100644 --- a/FIAT/__init__.py +++ b/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 diff --git a/FIAT/crouzeix_falk.py b/FIAT/crouzeix_falk.py new file mode 100644 index 000000000..82c94625a --- /dev/null +++ b/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) diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index b8fa3f80e..227d46e21 100644 --- a/FIAT/crouzeix_raviart.py +++ b/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) diff --git a/FIAT/symfem.py b/FIAT/symfem.py new file mode 100644 index 000000000..22d1c2d14 --- /dev/null +++ b/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