Test the GN space

FIAT/guzman_neilan.py

@@ -10,7 +10,7 @@
 # transformation theory.
 
 from FIAT import finite_element, polynomial_set, expansions
-from FIAT.bernardi_raugel import ExtendedBernardiRaugelSpace, BernardiRaugelDualSet
+from FIAT.bernardi_raugel import ExtendedBernardiRaugelSpace, BernardiRaugelDualSet, BernardiRaugel
 from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
 from FIAT.macro import AlfeldSplit
 from FIAT.quadrature_schemes import create_quadrature
@@ -20,15 +20,19 @@
 import math
 
 
-def ExtendedGuzmanNeilanSpace(ref_el, degree):
+def ExtendedGuzmanNeilanSpace(ref_el, degree, 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")
  B = take_interior_bubbles(C0)
  if sd > 2:
  B = modified_bubble_subspace(B)
- BR = ExtendedBernardiRaugelSpace(ref_el, degree)
+ if reduced:
+ BR = BernardiRaugel(ref_el, degree).get_nodal_basis().take(list(range((sd+1)**2)))
+ else:
+ BR = ExtendedBernardiRaugelSpace(ref_el, degree)
+
  GN = constant_div_projection(BR, C0, B)
  return GN
 
@@ -129,13 +133,13 @@ def constant_div_projection(BR, C0, M):
  # Invert the divergence
  B = inner(P, div(U), qwts)
  g = inner(P, div(X), qwts)
- u = numpy.linalg.solve(B, g)
+ w = numpy.linalg.solve(B, g)
 
  # Add correction to BR bubbles
  v = C0.tabulate(qpts)[(0,)*sd]
  coeffs = numpy.linalg.solve(inner(v, v, qwts), inner(v, X[(0,)*sd], qwts))
  coeffs = coeffs.T.reshape(BR.get_num_members(), sd, -1)
- coeffs -= numpy.tensordot(u, M.get_coeffs(), axes=(0, 0))
+ coeffs -= numpy.tensordot(w, M.get_coeffs(), axes=(0, 0))
  GN = polynomial_set.PolynomialSet(ref_complex, degree, degree,
  C0.get_expansion_set(), coeffs)
  return GN

test/unit/test_macro.py

@@ -383,20 +383,25 @@ def test_AlfeldSorokinaSpace(cell, degree):
  assert numpy.allclose(residual, 0)
 
 
-def test_minimal_stokes_space(cell):
+@pytest.mark.parametrize("family", ("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 == 3:
- degree = 1
- space = ChristiansenHuSpace
+ if sd == 1:
+ return
+ if family == "GN":
+ degree = sd
+ space = ExtendedGuzmanNeilanSpace
  elif sd == 2:
  degree = 2
  space = ArnoldQinSpace
- else:
- return
+ elif sd == 3:
+ degree = 1
+ space = ChristiansenHuSpace
 
  W = space(cell, degree)
  V = space(cell, degree, reduced=True)