BR Bubbles, and non-Piola mapped Arnold-Qin

FIAT/arnold_qin.py

@@ -4,7 +4,8 @@
 #
 # Written by Pablo D. Brubeck ([email protected]), 2024
 
-from FIAT import finite_element, polynomial_set
+from FIAT import finite_element, polynomial_set, dual_set
+from FIAT.functional import ComponentPointEvaluation
 from FIAT.bernardi_raugel import BernardiRaugelDualSet
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.reference_element import TRIANGLE
@@ -55,12 +56,35 @@ def ArnoldQinSpace(ref_el, degree, reduced=False):
  C0.get_expansion_set(), coeffs)
 
 
+class VectorLagrangeDualSet(dual_set.DualSet):
+
+ def __init__(self, ref_el, degree, variant=None):
+ sd = ref_el.get_spatial_dimension()
+ top = ref_el.get_topology()
+ entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+ # Point evaluation at lattice points
+ nodes = []
+ for dim in sorted(top):
+ for entity in sorted(top[dim]):
+ cur = len(nodes)
+ pts = ref_el.make_points(dim, entity, degree, variant=variant)
+ nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt)
+ for pt in pts for comp in range(sd))
+ entity_ids[dim][entity].extend(range(cur, len(nodes)))
+ super().__init__(nodes, ref_el, entity_ids)
+
+
 class ArnoldQin(finite_element.CiarletElement):
  """The Arnold-Qin C^0(Alfeld) quadratic macroelement with divergence in P0.
  This element belongs to a Stokes complex, and is paired with unsplit DG0."""
- def __init__(self, ref_el, degree=2):
+ def __init__(self, ref_el, degree=2, reduced=False):
  poly_set = ArnoldQinSpace(ref_el, degree)
- ref_complex = poly_set.get_reference_element()
- dual = BernardiRaugelDualSet(ref_complex, degree)
+ if reduced:
+ dual = BernardiRaugelDualSet(ref_el, degree)
+ mapping = "contravariant piola"
+ else:
+ dual = VectorLagrangeDualSet(ref_el, degree)
+ mapping = "affine"
  formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form
- super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")
+ super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)

FIAT/bernardi_raugel.py

@@ -16,37 +16,42 @@
 import numpy
 
 
-def ExtendedBernardiRaugelSpace(ref_el, degree):
+def ExtendedBernardiRaugelSpace(ref_el, subdegree):
  r"""Return a basis for the extended Bernardi-Raugel space.
- P_1^d + (P_{d} \ P_{d-1})^d"""
+ P_k^d + (P_{dim} \ P_{dim-1})^d"""
  sd = ref_el.get_spatial_dimension()
- Pk = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), scale=1, variant="bubble")
- dimPk = expansions.polynomial_dimension(ref_el, degree, continuity="C0")
- entity_ids = expansions.polynomial_entity_ids(ref_el, degree, continuity="C0")
- ids = [i + j * dimPk
- for dim in (0, sd-1)
+ if subdegree >= sd:
+ raise ValueError("The Bernardi-Raugel space is only defined for subdegree < dim")
+ Pd = polynomial_set.ONPolynomialSet(ref_el, sd, shape=(sd,), scale=1, variant="bubble")
+ dimPd = expansions.polynomial_dimension(ref_el, sd, continuity="C0")
+ entity_ids = expansions.polynomial_entity_ids(ref_el, sd, continuity="C0")
+ ids = [i + j * dimPd
+ for dim in (*tuple(range(subdegree)), sd-1)
  for f in sorted(entity_ids[dim])
  for i in entity_ids[dim][f]
  for j in range(sd)]
- return Pk.take(ids)
+ return Pd.take(ids)
 
 
 class BernardiRaugelDualSet(dual_set.DualSet):
- def __init__(self, ref_el, degree, reduced=False):
-  ref_complex = ref_el
+ """The Bernardi-Raugel dual set."""
+ def __init__(self, ref_complex, degree, subdegree=1, reduced=False):
  ref_el = ref_complex.get_parent() or ref_complex
  sd = ref_el.get_spatial_dimension()
  top = ref_el.get_topology()
  entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
 
+ # Point evaluation at lattice points
  nodes = []
- for v in sorted(top[0]):
- cur = len(nodes)
- pt, = ref_el.make_points(0, v, degree)
- nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt)
- for k in range(sd))
- entity_ids[0][v].extend(range(cur, len(nodes)))
+ for dim in range(subdegree):
+ for entity in sorted(top[dim]):
+ cur = len(nodes)
+ pts = ref_el.make_points(dim, entity, subdegree)
+ nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt)
+ for pt in pts for comp in range(sd))
+ entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
+ # Face moments of normal/tangential components against mean-free bubbles
  facet = ref_complex.construct_subcomplex(sd-1)
  Q = create_quadrature(facet, 2*degree)
  if degree == 1 and facet.is_macrocell():
@@ -55,7 +60,6 @@ def __init__(self, ref_el, degree, reduced=False):
  else:
  ref_facet = facet.get_parent() or facet
  f_at_qpts = ref_facet.compute_bubble(Q.get_points())
-
  f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume()
 
  Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q)
@@ -83,13 +87,13 @@ def __init__(self, ref_el, degree, reduced=False):
 
 class BernardiRaugel(finite_element.CiarletElement):
  """The Bernardi-Raugel extended element."""
- def __init__(self, ref_el, degree=None):
+ def __init__(self, ref_el, degree=None, subdegree=1):
  sd = ref_el.get_spatial_dimension()
  if degree is None:
  degree = sd
  if degree != sd:
  raise ValueError("Bernardi-Raugel only defined for degree = dim")
- poly_set = ExtendedBernardiRaugelSpace(ref_el, degree)
- dual = BernardiRaugelDualSet(ref_el, degree)
- formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form
+ poly_set = ExtendedBernardiRaugelSpace(ref_el, subdegree)
+ dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree)
+ formdegree = sd - 1 # (n-1)-form
  super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

FIAT/guzman_neilan.py

@@ -20,33 +20,35 @@
 import math
 
 
-def ExtendedGuzmanNeilanSpace(ref_el, degree, reduced=False):
+def ExtendedGuzmanNeilanSpace(ref_el, subdegree, reduced=False):
  """Return a basis for the extended Guzman-Neilan space."""
  sd = ref_el.get_spatial_dimension()
  ref_complex = AlfeldSplit(ref_el)
- C0 = polynomial_set.ONPolynomialSet(ref_complex, degree, shape=(sd,), scale=1, variant="bubble")
+ C0 = polynomial_set.ONPolynomialSet(ref_complex, sd, shape=(sd,), scale=1, variant="bubble")
  B = take_interior_bubbles(C0)
  if sd > 2:
  B = modified_bubble_subspace(B)
+
  if reduced:
- BR = BernardiRaugel(ref_el, degree).get_nodal_basis().take(list(range((sd+1)**2)))
+ BR = BernardiRaugel(ref_el, sd, subdegree=subdegree).get_nodal_basis()
+ reduced_dim = BR.get_num_members() - (sd-1) * (sd+1)
+ BR = BR.take(list(range(reduced_dim)))
  else:
- BR = ExtendedBernardiRaugelSpace(ref_el, degree)
-
+ BR = ExtendedBernardiRaugelSpace(ref_el, subdegree)
  GN = constant_div_projection(BR, C0, B)
  return GN
 
 
 class GuzmanNeilan(finite_element.CiarletElement):
  """The Guzman-Neilan extended element."""
- def __init__(self, ref_el, degree=None):
+ def __init__(self, ref_el, degree=None, subdegree=1):
  sd = ref_el.get_spatial_dimension()
  if degree is None:
  degree = sd
  if degree != sd:
  raise ValueError("Guzman-Neilan only defined for degree = dim")
- poly_set = ExtendedGuzmanNeilanSpace(ref_el, degree)
- dual = BernardiRaugelDualSet(ref_el, degree)
+ poly_set = ExtendedGuzmanNeilanSpace(ref_el, subdegree)
+ dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree)
  formdegree = sd - 1 # (n-1)-form
  super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")
 

test/unit/test_hct.py

@@ -2,11 +2,6 @@
 import numpy
 
 from FIAT import HsiehCloughTocher as HCT
-from FIAT import AlfeldSorokina as AS
-from FIAT import ArnoldQin as AQ
-from FIAT import Lagrange as CG
-from FIAT import DiscontinuousLagrange as DG
-from FIAT.restricted import RestrictedElement
 from FIAT.reference_element import ufc_simplex, make_lattice
 from FIAT.functional import PointEvaluation
 
@@ -38,44 +33,3 @@ def test_hct_constant(cell, reduced):
  expected = 1 if sum(alpha) == 0 else 0
  vals = numpy.dot(coefs, tab[alpha])
  assert numpy.allclose(vals, expected)
-
-
-@pytest.mark.parametrize("reduced", (False, True))
-def test_stokes_complex(cell, reduced):
- # Test that we have the lowest-order discrete Stokes complex, verifying
- # that the range of the exterior derivative of each space is contained by
- # the next space in the sequence
- H2 = HCT(cell, reduced=reduced)
- ref_complex = H2.get_reference_complex()
- if reduced:
- # HCT-red(3) -curl-> AQ-red(2) -div-> DG(0)
- H2 = RestrictedElement(H2, restriction_domain="vertex")
- H1 = RestrictedElement(AQ(cell), indices=list(range(9)))
- L2 = DG(cell, 0)
- else:
- # HCT(3) -curl-> AS(2) -div-> CG(1, Alfeld)
- H1 = AS(cell)
- L2 = CG(ref_complex, 1)
-
- pts = []
- top = ref_complex.get_topology()
- for dim in top:
- for entity in top[dim]:
- pts.extend(ref_complex.make_points(dim, entity, 4))
-
- H2tab = H2.tabulate(1, pts)
- H1tab = H1.tabulate(1, pts)
- L2tab = L2.tabulate(0, pts)
-
- L2val = L2tab[(0, 0)]
- H1val = H1tab[(0, 0)]
- H1div = sum(H1tab[alpha][:, alpha.index(1), :] for alpha in H1tab if sum(alpha) == 1)
- H2curl = numpy.stack([H2tab[(0, 1)], -H2tab[(1, 0)]], axis=1)
-
- H2dim = H2.space_dimension()
- H1dim = H1.space_dimension()
- _, residual, *_ = numpy.linalg.lstsq(H1val.reshape(H1dim, -1).T, H2curl.reshape(H2dim, -1).T)
- assert numpy.allclose(residual, 0)
-
- _, residual, *_ = numpy.linalg.lstsq(L2val.T, H1div.T)
- assert numpy.allclose(residual, 0)

test/unit/test_macro.py

@@ -10,9 +10,9 @@
 from FIAT.barycentric_interpolation import get_lagrange_points
 
 
-@pytest.fixture(params=("I", "T", "S"))
+@pytest.fixture(params=(1, 2, 3), ids=("I", "T", "S"))
 def cell(request):
- dim = {"I": 1, "T": 2, "S": 3}[request.param]
+ dim = request.param
  return ufc_simplex(dim)
 
 
@@ -358,93 +358,6 @@ def test_Ck_basis(cell, order, degree, variant):
  assert numpy.allclose(local_phis, phis[:, ipts])
 
 
-@pytest.mark.parametrize("degree", (2, 3, 4))
-def test_AlfeldSorokinaSpace(cell, degree):
- # Test that the divergence of the Alfeld-Sorokina space is spanned by a C0 basis
- from FIAT.alfeld_sorokina import AlfeldSorokinaSpace
-
- P1 = AlfeldSorokinaSpace(cell, degree)
- A = P1.get_reference_element()
- top = A.get_topology()
- sd = A.get_spatial_dimension()
-
- pts = []
- for dim in top:
- for entity in top[dim]:
- pts.extend(A.make_points(dim, entity, degree))
-
- P1_tab = P1.tabulate(pts, 1)
- divP1_tab = sum(P1_tab[alpha][:, alpha.index(1), :]
- for alpha in P1_tab if sum(alpha) == 1)
-
- P2 = CkPolynomialSet(A, degree-1, order=0, variant="bubble")
- P2_tab = P2.tabulate(pts)[(0,)*sd]
- _, residual, *_ = numpy.linalg.lstsq(P2_tab.T, divP1_tab.T)
- assert numpy.allclose(residual, 0)
-
-
-@pytest.mark.parametrize("family", ("AQ", "CH", "GN"))
-def test_minimal_stokes_space(cell, family):
- # Test that the C0 Stokes space is spanned by a C0 basis
- # Also test that its divergence is constant
- from FIAT.arnold_qin import ArnoldQinSpace
- from FIAT.christiansen_hu import ChristiansenHuSpace
- from FIAT.guzman_neilan import ExtendedGuzmanNeilanSpace
- sd = cell.get_spatial_dimension()
- if sd == 1:
- return
- if family == "GN":
- degree = sd
- space = ExtendedGuzmanNeilanSpace
- elif family == "CH":
- degree = 1
- space = ChristiansenHuSpace
- elif family == "AQ":
- if sd != 2:
- return
- degree = 2
- space = ArnoldQinSpace
-
- W = space(cell, degree)
- V = space(cell, degree, reduced=True)
- Wdim = W.get_num_members()
- Vdim = V.get_num_members()
- K = W.get_reference_element()
- sd = K.get_spatial_dimension()
- top = K.get_topology()
-
- pts = []
- for dim in top:
- for entity in top[dim]:
- pts.extend(K.make_points(dim, entity, degree))
-
- C0 = CkPolynomialSet(K, degree, order=0, variant="bubble")
- C0_tab = C0.tabulate(pts)
- Wtab = W.tabulate(pts, 1)
- Vtab = V.tabulate(pts, 1)
- z = (0,)*sd
- for tab in (Vtab, Wtab):
- # Test that the space is full rank
- _, sig, _ = numpy.linalg.svd(tab[z].reshape(-1, sd*len(pts)).T, full_matrices=True)
- assert all(sig > 1E-10)
-
- # Test that the space is C0
- for k in range(sd):
- _, residual, *_ = numpy.linalg.lstsq(C0_tab[z].T, tab[z][:, k, :].T)
- assert numpy.allclose(residual, 0)
-
- # Test that divergence is in P0
- div = sum(tab[alpha][:, alpha.index(1), :]
- for alpha in tab if sum(alpha) == 1)[:Vdim]
- assert numpy.allclose(div, div[:, 0][:, None])
-
- # Test that the full space includes the reduced space
- assert Wdim > Vdim
- _, residual, *_ = numpy.linalg.lstsq(Wtab[z].reshape(Wdim, -1).T,
- Vtab[z].reshape(Vdim, -1).T)
- assert numpy.allclose(residual, 0)
-
-
 def test_distance_to_point_l1(cell):
  A = AlfeldSplit(cell)
  dim = A.get_spatial_dimension()