Skip to content
This repository has been archived by the owner on Feb 21, 2022. It is now read-only.

Commit

Permalink
make the same changes for BDM and N2curl
Browse files Browse the repository at this point in the history
  • Loading branch information
@FabianL1908 @wence-
FabianL1908 authored and wence- committed Jun 10, 2020
1 parent 57ca6db commit cb9b8db
Show file tree
Hide file tree
Showing 2 changed files with 140 additions and 47 deletions.
102 changes: 77 additions & 25 deletions FIAT/brezzi_douglas_marini.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@


class BDMDualSet(dual_set.DualSet):
def __init__(self, ref_el, degree):
def __init__(self, ref_el, degree, variant, quad_deg):

# Initialize containers for map: mesh_entity -> dof number and
# dual basis
Expand All @@ -20,28 +20,55 @@ def __init__(self, ref_el, degree):
sd = ref_el.get_spatial_dimension()
t = ref_el.get_topology()

# Define each functional for the dual set
# codimension 1 facets
for i in range(len(t[sd - 1])):
pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
for j in range(len(pts_cur)):
pt_cur = pts_cur[j]
f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
nodes.append(f)

# internal nodes
if degree > 1:
Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
qpts = Q.get_points()
Nedel = nedelec.Nedelec(ref_el, degree - 1)
Nedfs = Nedel.get_nodal_basis()
zero_index = tuple([0 for i in range(sd)])
Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]

for i in range(len(Ned_at_qpts)):
phi_cur = Ned_at_qpts[i, :]
l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
nodes.append(l_cur)
if variant == "integral":
facet = ref_el.get_facet_element()
# Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
# degree is q - 1
Q = quadrature.make_quadrature(facet, quad_deg)
Pq = polynomial_set.ONPolynomialSet(facet, degree)
Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))]
for f in range(len(t[sd - 1])):
for i in range(Pq_at_qpts.shape[0]):
phi = Pq_at_qpts[i, :]
nodes.append(functional.IntegralMomentOfScaledNormalEvaluation(ref_el, Q, phi, f))

# internal nodes
if degree > 1:
Q = quadrature.make_quadrature(ref_el, quad_deg)
qpts = Q.get_points()
Nedel = nedelec.Nedelec(ref_el, degree - 1, variant)
Nedfs = Nedel.get_nodal_basis()
zero_index = tuple([0 for i in range(sd)])
Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]

for i in range(len(Ned_at_qpts)):
phi_cur = Ned_at_qpts[i, :]
l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
nodes.append(l_cur)

elif variant == "point":
# Define each functional for the dual set
# codimension 1 facets
for i in range(len(t[sd - 1])):
pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
for j in range(len(pts_cur)):
pt_cur = pts_cur[j]
f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
nodes.append(f)

# internal nodes
if degree > 1:
Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
qpts = Q.get_points()
Nedel = nedelec.Nedelec(ref_el, degree - 1, variant)
Nedfs = Nedel.get_nodal_basis()
zero_index = tuple([0 for i in range(sd)])
Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]

for i in range(len(Ned_at_qpts)):
phi_cur = Ned_at_qpts[i, :]
l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
nodes.append(l_cur)

# sets vertices (and in 3d, edges) to have no nodes
for i in range(sd - 1):
Expand Down Expand Up @@ -73,14 +100,39 @@ def __init__(self, ref_el, degree):
class BrezziDouglasMarini(finite_element.CiarletElement):
"""The BDM element"""

def __init__(self, ref_el, degree):
def __init__(self, ref_el, degree, variant=None):

if variant is None:
variant = "point"
print('Warning: Variant of BDM element will change from point evaluation to integral evaluation.'
'You should project into variant="integral"')
#Replace by the following in a month time
#variant = "integral"

if not (variant == "point" or "integral" in variant):
raise ValueError('Choose either variant="point" or variant="integral"'
'or variant="integral(Quadrature degree)"')

if variant == "integral":
quad_deg = 5 * (degree + 1)
variant = "integral"
elif "integral" in variant:
try:
quad_deg = int(''.join(filter(str.isdigit, variant)))
except:
raise ValueError("Wrong format for variant")
if quad_deg < degree + 1:
raise ValueError("Warning, quadrature degree should be at least %s" % (degree + 1))
variant = "integral"
elif variant == "point":
quad_deg = None

if degree < 1:
raise Exception("BDM_k elements only valid for k >= 1")

sd = ref_el.get_spatial_dimension()
poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, ))
dual = BDMDualSet(ref_el, degree)
dual = BDMDualSet(ref_el, degree, variant, quad_deg)
formdegree = sd - 1 # (n-1)-form
super(BrezziDouglasMarini, self).__init__(poly_set, dual, degree, formdegree,
mapping="contravariant piola")
85 changes: 63 additions & 22 deletions FIAT/nedelec_second_kind.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@
from FIAT.quadrature import make_quadrature, UFCTetrahedronFaceQuadratureRule
from FIAT.reference_element import UFCTetrahedron

from FIAT import (polynomial_set, expansions, quadrature, dual_set,
finite_element, functional)


class NedelecSecondKindDual(DualSet):
r"""
Expand Down Expand Up @@ -46,15 +49,14 @@ class NedelecSecondKindDual(DualSet):
these elements coincide with the CG_k elements.)
"""

def __init__(self, cell, degree):
def __init__(self, cell, degree, variant, quad_deg):

# Define degrees of freedom
(dofs, ids) = self.generate_degrees_of_freedom(cell, degree)

(dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, quad_deg)
# Call init of super-class
super(NedelecSecondKindDual, self).__init__(dofs, cell, ids)

def generate_degrees_of_freedom(self, cell, degree):
def generate_degrees_of_freedom(self, cell, degree, variant, quad_deg):
"Generate dofs and geometry-to-dof maps (ids)."

dofs = []
Expand All @@ -68,47 +70,61 @@ def generate_degrees_of_freedom(self, cell, degree):
ids[0] = dict(list(zip(list(range(d + 1)), ([] for i in range(d + 1)))))

# (d+1) degrees of freedom per entity of codimension 1 (edges)
(edge_dofs, edge_ids) = self._generate_edge_dofs(cell, degree, 0)
(edge_dofs, edge_ids) = self._generate_edge_dofs(cell, degree, 0, variant, quad_deg)
dofs.extend(edge_dofs)
ids[1] = edge_ids

# Include face degrees of freedom if 3D
if d == 3:
(face_dofs, face_ids) = self._generate_face_dofs(cell, degree,
len(dofs))
len(dofs), variant, quad_deg)
dofs.extend(face_dofs)
ids[2] = face_ids

# Varying degrees of freedom (possibly zero) per cell
(cell_dofs, cell_ids) = self._generate_cell_dofs(cell, degree, len(dofs))
(cell_dofs, cell_ids) = self._generate_cell_dofs(cell, degree, len(dofs), variant, quad_deg)
dofs.extend(cell_dofs)
ids[d] = cell_ids

return (dofs, ids)

def _generate_edge_dofs(self, cell, degree, offset):
def _generate_edge_dofs(self, cell, degree, offset, variant, quad_deg):
"""Generate degrees of freedoms (dofs) for entities of
codimension 1 (edges)."""

# (degree+1) tangential component point evaluation degrees of
# freedom per entity of codimension 1 (edges)
dofs = []
ids = {}
for edge in range(len(cell.get_topology()[1])):

# Create points for evaluation of tangential components
points = cell.make_points(1, edge, degree + 2)
if variant == "integral":
edge = cell.get_facet_element().get_facet_element()
Q = quadrature.make_quadrature(edge, quad_deg)
Pq = polynomial_set.ONPolynomialSet(edge, degree)
Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(1))]
for e in range(len(cell.get_topology()[1])):
for i in range(Pq_at_qpts.shape[0]):
phi = Pq_at_qpts[i, :]
dofs.append(functional.IntegralMomentOfEdgeTangentEvaluation(cell, Q, phi, e))
jj = Pq_at_qpts.shape[0] * e
ids[e] = list(range(offset + jj, offset + jj + Pq_at_qpts.shape[0]))

elif variant == "point":
for edge in range(len(cell.get_topology()[1])):

# Create points for evaluation of tangential components
points = cell.make_points(1, edge, degree + 2)

# A tangential component evaluation for each point
dofs += [Tangent(cell, edge, point) for point in points]
# A tangential component evaluation for each point
dofs += [Tangent(cell, edge, point) for point in points]

# Associate these dofs with this edge
i = len(points) * edge
ids[edge] = list(range(offset + i, offset + i + len(points)))
# Associate these dofs with this edge
i = len(points) * edge
ids[edge] = list(range(offset + i, offset + i + len(points)))

return (dofs, ids)

def _generate_face_dofs(self, cell, degree, offset):
def _generate_face_dofs(self, cell, degree, offset, variant, quad_deg):
"""Generate degrees of freedoms (dofs) for faces."""

# Initialize empty dofs and identifiers (ids)
Expand All @@ -134,7 +150,7 @@ def _generate_face_dofs(self, cell, degree, offset):
# Construct Raviart-Thomas of (degree - 1) on the
# reference face
reference_face = Q_face.reference_rule().ref_el
RT = RaviartThomas(reference_face, degree - 1)
RT = RaviartThomas(reference_face, degree - 1, variant)
num_rts = RT.space_dimension()

# Evaluate RT basis functions at reference quadrature
Expand Down Expand Up @@ -168,7 +184,7 @@ def _generate_face_dofs(self, cell, degree, offset):

return (dofs, ids)

def _generate_cell_dofs(self, cell, degree, offset):
def _generate_cell_dofs(self, cell, degree, offset, variant, quad_deg):
"""Generate degrees of freedoms (dofs) for entities of
codimension d (cells)."""

Expand All @@ -182,7 +198,7 @@ def _generate_cell_dofs(self, cell, degree, offset):
qs = Q.get_points()

# Create Raviart-Thomas nodal basis
RT = RaviartThomas(cell, degree + 1 - d)
RT = RaviartThomas(cell, degree + 1 - d, variant)
phi = RT.get_nodal_basis()

# Evaluate Raviart-Thomas basis at quadrature points
Expand All @@ -205,7 +221,32 @@ class NedelecSecondKind(CiarletElement):
H(curl) conformity.
"""

def __init__(self, cell, degree):
def __init__(self, cell, degree, variant=None):

if variant is None:
variant = "point"
print('Warning: Variant of Nedelec 2nd kind element will change from point evaluation to integral evaluation.'
'You should project into variant="integral"')
#Replace by the following in a month time
#variant = "integral"

if not (variant == "point" or "integral" in variant):
raise ValueError('Choose either variant="point" or variant="integral"'
'or variant="integral(Quadrature degree)"')

if variant == "integral":
quad_deg = 5 * (degree + 1)
variant = "integral"
elif "integral" in variant:
try:
quad_deg = int(''.join(filter(str.isdigit, variant)))
except:
raise ValueError("Wrong format for variant")
if quad_deg < degree + 1:
raise ValueError("Warning, quadrature degree should be at least %s" % (degree + 1))
variant = "integral"
elif variant == "point":
quad_deg = None

# Check degree
assert degree >= 1, "Second kind Nedelecs start at 1!"
Expand All @@ -217,7 +258,7 @@ def __init__(self, cell, degree):
Ps = ONPolynomialSet(cell, degree, (d, ))

# Construct dual space
Ls = NedelecSecondKindDual(cell, degree)
Ls = NedelecSecondKindDual(cell, degree, variant, quad_deg)

# Set form degree
formdegree = 1 # 1-form
Expand Down

0 comments on commit cb9b8db

Please sign in to comment.