Skip to content

Commit

Permalink
Integral variant for Regge/HHJ (#96)
Browse files Browse the repository at this point in the history 
* docstring

* cleanup test_fiat.py

* Regge integral dofs

* HHJ integral dofs in 2d/3d

* support variant="integral(q)"

* BidirectionalMoment dofs

* address review comments

* Cleanup H(div, S) elements (#98)

* implementation of Hu-Zhang element

---------

Co-authored-by: Francis Aznaran <[email protected]>
  • Loading branch information
@pbrubeck @FAznaran
pbrubeck and FAznaran authored Nov 19, 2024
1 parent c6369c7 commit b063854
Show file tree
Hide file tree
Showing 20 changed files with 909 additions and 671 deletions.
33 changes: 17 additions & 16 deletions FIAT/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,17 @@
evaluating arbitrary order Lagrange and many other elements.
Simplices in one, two, and three dimensions are supported."""

import pkg_resources
# Important functionality
from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401
from FIAT.quadrature import make_quadrature # noqa: F401
from FIAT.quadrature_schemes import create_quadrature # noqa: F401
from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401
from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401
from FIAT.mixed import MixedElement # noqa: F401
from FIAT.restricted import RestrictedElement # noqa: F401
from FIAT.quadrature_element import QuadratureElement # noqa: F401

# Import finite element classes
from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401
from FIAT.argyris import Argyris
from FIAT.bernardi_raugel import BernardiRaugel
from FIAT.bernstein import Bernstein
Expand All @@ -17,8 +24,7 @@
from FIAT.christiansen_hu import ChristiansenHu
from FIAT.johnson_mercier import JohnsonMercier
from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401
from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401
from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace # noqa: F401
from FIAT.SminusDiv import TrimmedSerendipityDiv # noqa: F401
from FIAT.SminusCurl import TrimmedSerendipityCurl # noqa: F401
from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini
Expand All @@ -42,30 +48,22 @@
from FIAT.raviart_thomas import RaviartThomas
from FIAT.crouzeix_raviart import CrouzeixRaviart
from FIAT.regge import Regge
from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind
from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind
from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
from FIAT.arnold_winther import ArnoldWinther
from FIAT.arnold_winther import ArnoldWintherNC
from FIAT.hu_zhang import HuZhang
from FIAT.mardal_tai_winther import MardalTaiWinther
from FIAT.bubble import Bubble, FacetBubble
from FIAT.tensor_product import TensorProductElement
from FIAT.enriched import EnrichedElement
from FIAT.nodal_enriched import NodalEnrichedElement
from FIAT.discontinuous import DiscontinuousElement
from FIAT.hdiv_trace import HDivTrace
from FIAT.mixed import MixedElement # noqa: F401
from FIAT.restricted import RestrictedElement # noqa: F401
from FIAT.quadrature_element import QuadratureElement # noqa: F401
from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401
from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen
from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401

# Important functionality
from FIAT.quadrature import make_quadrature # noqa: F401
from FIAT.quadrature_schemes import create_quadrature # noqa: F401
from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401
from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401

__version__ = pkg_resources.get_distribution("fenics-fiat").version

# List of supported elements and mapping to element classes
supported_elements = {"Argyris": Argyris,
"Bell": Bell,
Expand Down Expand Up @@ -112,8 +110,11 @@
"BrokenElement": DiscontinuousElement,
"HDiv Trace": HDivTrace,
"Hellan-Herrmann-Johnson": HellanHerrmannJohnson,
"Gopalakrishnan-Lederer-Schoberl 1st kind": GopalakrishnanLedererSchoberlFirstKind,
"Gopalakrishnan-Lederer-Schoberl 2nd kind": GopalakrishnanLedererSchoberlSecondKind,
"Conforming Arnold-Winther": ArnoldWinther,
"Nonconforming Arnold-Winther": ArnoldWintherNC,
"Hu-Zhang": HuZhang,
"Mardal-Tai-Winther": MardalTaiWinther}

# List of extra elements
Expand Down
228 changes: 98 additions & 130 deletions FIAT/arnold_winther.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,157 +9,125 @@
# SPDX-License-Identifier: LGPL-3.0-or-later


from FIAT.finite_element import CiarletElement
from FIAT.dual_set import DualSet
from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet
from FIAT.functional import (
PointwiseInnerProductEvaluation as InnerProduct,
FrobeniusIntegralMoment as FIM,
IntegralMomentOfTensorDivergence,
IntegralLegendreNormalNormalMoment,
IntegralLegendreNormalTangentialMoment)

from FIAT.quadrature import make_quadrature
from FIAT import finite_element, dual_set, polynomial_set
from FIAT.reference_element import TRIANGLE
from FIAT.quadrature_schemes import create_quadrature
from FIAT.functional import (ComponentPointEvaluation,
TensorBidirectionalIntegralMoment,
IntegralMomentOfTensorDivergence,
IntegralLegendreNormalNormalMoment,
IntegralLegendreNormalTangentialMoment)

import numpy


class ArnoldWintherNCDual(DualSet):
def __init__(self, cell, degree=2):
if not degree == 2:
raise ValueError("Nonconforming Arnold-Winther elements are"
"only defined for degree 2.")
dofs = []
dof_ids = {}
dof_ids[0] = {0: [], 1: [], 2: []}
dof_ids[1] = {0: [], 1: [], 2: []}
dof_ids[2] = {0: []}
class ArnoldWintherNCDual(dual_set.DualSet):
def __init__(self, ref_el, degree=2):
if degree != 2:
raise ValueError("Nonconforming Arnold-Winther elements are only defined for degree 2.")
top = ref_el.get_topology()
sd = ref_el.get_spatial_dimension()
entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
nodes = []

dof_cur = 0
# no vertex dofs
# proper edge dofs now (not the contraints)
# moments of normal . sigma against constants and linears.
for entity_id in range(3): # a triangle has 3 edges
for order in (0, 1):
dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
dof_cur += 4
# edge dofs: bidirectional nn and nt moments against P1.
qdegree = degree + 2
for entity in sorted(top[1]):
cur = len(nodes)
for order in range(2):
nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree))
nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree))
entity_ids[1][entity].extend(range(cur, len(nodes)))

# internal dofs: constant moments of three unique components
Q = make_quadrature(cell, 2)

e1 = numpy.array([1.0, 0.0]) # euclidean basis 1
e2 = numpy.array([0.0, 1.0]) # euclidean basis 2
basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices
for (v1, v2) in basis:
v1v2t = numpy.outer(v1, v2)
fatqp = numpy.zeros((2, 2, len(Q.pts)))
for i, y in enumerate(v1v2t):
for j, x in enumerate(y):
for k in range(len(Q.pts)):
fatqp[i, j, k] = x
dofs.append(FIM(cell, Q, fatqp))
dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
dof_cur += 3
cur = len(nodes)
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_el, degree)
phi = numpy.full(Q.get_weights().shape, 1/ref_el.volume())
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for i in range(sd) for j in range(i, sd))
entity_ids[2][0].extend(range(cur, len(nodes)))

# put the constraint dofs last.
for entity_id in range(3):
dof = IntegralLegendreNormalNormalMoment(cell, entity_id, 2, 6)
dofs.append(dof)
dof_ids[1][entity_id].append(dof_cur)
dof_cur += 1
for entity in sorted(top[1]):
cur = len(nodes)
nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, 2, qdegree))
entity_ids[1][entity].append(cur)

super().__init__(dofs, cell, dof_ids)
super().__init__(nodes, ref_el, entity_ids)


class ArnoldWintherNC(CiarletElement):
class ArnoldWintherNC(finite_element.CiarletElement):
"""The definition of the nonconforming Arnold-Winther element.
"""
def __init__(self, cell, degree=2):
assert degree == 2, "Only defined for degree 2"
Ps = ONSymTensorPolynomialSet(cell, degree)
Ls = ArnoldWintherNCDual(cell, degree)
def __init__(self, ref_el, degree=2):
if ref_el.shape != TRIANGLE:
raise ValueError(f"{type(self).__name__} only defined on triangles")
Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
Ls = ArnoldWintherNCDual(ref_el, degree)
formdegree = ref_el.get_spatial_dimension() - 1
mapping = "double contravariant piola"
super().__init__(Ps, Ls, degree, formdegree, mapping=mapping)

super().__init__(Ps, Ls, degree, mapping=mapping)


class ArnoldWintherDual(DualSet):
def __init__(self, cell, degree=3):
if not degree == 3:
raise ValueError("Arnold-Winther elements are"
"only defined for degree 3.")
dofs = []
dof_ids = {}
dof_ids[0] = {0: [], 1: [], 2: []}
dof_ids[1] = {0: [], 1: [], 2: []}
dof_ids[2] = {0: []}

dof_cur = 0
class ArnoldWintherDual(dual_set.DualSet):
def __init__(self, ref_el, degree=3):
if degree != 3:
raise ValueError("Arnold-Winther elements are only defined for degree 3.")
top = ref_el.get_topology()
sd = ref_el.get_spatial_dimension()
shp = (sd, sd)
entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
nodes = []

# vertex dofs
vs = cell.get_vertices()
e1 = numpy.array([1.0, 0.0])
e2 = numpy.array([0.0, 1.0])
basis = [(e1, e1), (e1, e2), (e2, e2)]

dof_cur = 0

for entity_id in range(3):
node = tuple(vs[entity_id])
for (v1, v2) in basis:
dofs.append(InnerProduct(cell, v1, v2, node))
dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3))
dof_cur += 3

# edge dofs now
# moments of normal . sigma against constants and linears.
for entity_id in range(3):
for order in (0, 1):
dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
dof_cur += 4

# internal dofs: constant moments of three unique components
Q = make_quadrature(cell, 3)

e1 = numpy.array([1.0, 0.0]) # euclidean basis 1
e2 = numpy.array([0.0, 1.0]) # euclidean basis 2
basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices
for (v1, v2) in basis:
v1v2t = numpy.outer(v1, v2)
fatqp = numpy.zeros((2, 2, len(Q.pts)))
for k in range(len(Q.pts)):
fatqp[:, :, k] = v1v2t
dofs.append(FIM(cell, Q, fatqp))
dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
dof_cur += 3

# Constraint dofs

Q = make_quadrature(cell, 5)

onp = ONPolynomialSet(cell, 2, (2,))
pts = Q.get_points()
onpvals = onp.tabulate(pts)[0, 0]

for i in list(range(3, 6)) + list(range(9, 12)):
dofs.append(IntegralMomentOfTensorDivergence(cell, Q,
onpvals[i, :, :]))

dof_ids[2][0] += list(range(dof_cur, dof_cur+6))

super().__init__(dofs, cell, dof_ids)


class ArnoldWinther(CiarletElement):
for v in sorted(top[0]):
cur = len(nodes)
pt, = ref_el.make_points(0, v, degree)
nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt)
for i in range(sd) for j in range(i, sd))
entity_ids[0][v].extend(range(cur, len(nodes)))

# edge dofs: bidirectional nn and nt moments against P_{k-2}
max_order = degree - 2
qdegree = degree + max_order
for entity in sorted(top[1]):
cur = len(nodes)
for order in range(max_order+1):
nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree))
nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree))
entity_ids[1][entity].extend(range(cur, len(nodes)))

# internal dofs: moments of unique components against P_{k-3}
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_el, 2*(degree-1))
P = polynomial_set.ONPolynomialSet(ref_el, degree-3, scale="L2 piola")
phis = P.tabulate(Q.get_points())[(0,)*sd]
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for phi in phis for i in range(sd) for j in range(i, sd))

# constraint dofs: moments of divergence against P_{k-1} \ P_{k-2}
P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,))
dimPkm1 = P.expansion_set.get_num_members(degree-1)
dimPkm2 = P.expansion_set.get_num_members(degree-2)
PH = P.take([i + j * dimPkm1 for j in range(sd) for i in range(dimPkm2, dimPkm1)])
phis = PH.tabulate(Q.get_points())[(0,)*sd]
nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis)

entity_ids[2][0].extend(range(cur, len(nodes)))
super().__init__(nodes, ref_el, entity_ids)


class ArnoldWinther(finite_element.CiarletElement):
"""The definition of the conforming Arnold-Winther element.
"""
def __init__(self, cell, degree=3):
assert degree == 3, "Only defined for degree 3"
Ps = ONSymTensorPolynomialSet(cell, degree)
Ls = ArnoldWintherDual(cell, degree)
def __init__(self, ref_el, degree=3):
if ref_el.shape != TRIANGLE:
raise ValueError(f"{type(self).__name__} only defined on triangles")
Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
Ls = ArnoldWintherDual(ref_el, degree)
formdegree = ref_el.get_spatial_dimension() - 1
mapping = "double contravariant piola"
super().__init__(Ps, Ls, degree, mapping=mapping)
super().__init__(Ps, Ls, degree, formdegree, mapping=mapping)
13 changes: 7 additions & 6 deletions FIAT/expansions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,16 +279,19 @@ def __init__(self, ref_el, scale=None, variant=None):
self._dmats_cache = {}
self._cell_node_map_cache = {}

def get_scale(self, cell=0):
def get_scale(self, n, cell=0):
scale = self.scale
sd = self.ref_el.get_spatial_dimension()
if isinstance(scale, str):
sd = self.ref_el.get_spatial_dimension()
vol = self.ref_el.volume_of_subcomplex(sd, cell)
scale = scale.lower()
if scale == "orthonormal":
scale = math.sqrt(1.0 / vol)
elif scale == "l2 piola":
scale = 1.0 / vol
elif n == 0 and sd > 1 and len(self.affine_mappings) == 1:
# return 1 for n=0 to make regression tests pass
scale = 1
return scale

def get_num_members(self, n):
Expand All @@ -310,9 +313,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
ref_pts = numpy.add(numpy.dot(pts, A.T), b).T
Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
sd = self.ref_el.get_spatial_dimension()

# Always return 1 for n=0 to make regression tests pass
scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell)
scale = self.get_scale(n, cell=cell)
phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
scale, variant=self.variant)
if self.continuity == "C0":
Expand Down Expand Up @@ -549,7 +550,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
Jinv = A[0, 0] if direction is None else numpy.dot(A, direction)
xs = numpy.add(numpy.dot(pts, A.T), b)
results = {}
scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
scale = self.get_scale(n, cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
for k in range(order+1):
v = numpy.zeros((n + 1, len(xs)), xs.dtype)
if n >= k:
Expand Down
Loading

0 comments on commit b063854

Please sign in to comment.