Test Alfeld-Sorokina (#3771)

* Test Alfeld-Sorokina

* Skip derivative nodes

* Tests for AS, AQ, and AQ-red macroelements


* Tests for AS, AQ, and AQ-red macroelements

firedrake/bcs.py

@@ -134,11 +134,19 @@ def nodes(self):
             raise NotImplementedError("Strong BCs not implemented for element %r, use Nitsche-type methods until we figure this out" % V.finat_element)
 
         def hermite_stride(bcnodes):
-            if isinstance(self._function_space.finat_element, finat.Hermite) and \
-               self._function_space.mesh().topological_dimension() == 1:
-                return bcnodes[::2]  # every second dof is the vertex value
-            else:
-                return bcnodes
+            fe = self._function_space.finat_element
+            tdim = self._function_space.mesh().topological_dimension()
+            if isinstance(fe, finat.Hermite) and tdim == 1:
+                bcnodes = bcnodes[::2]  # every second dof is the vertex value
+            elif fe.complex.is_macrocell() and self._function_space.ufl_element().sobolev_space == ufl.H1:
+                # Skip derivative nodes for supersmooth H1 functions
+                nodes = fe.fiat_equivalent.dual_basis()
+                deriv_nodes = [i for i, node in enumerate(nodes)
+                               if len(node.deriv_dict) != 0]
+                if len(deriv_nodes) > 0:
+                    deriv_ids = self._function_space.cell_node_list[:, deriv_nodes]
+                    bcnodes = np.setdiff1d(bcnodes, deriv_ids)
+            return bcnodes
 
         sub_d = (self.sub_domain, ) if isinstance(self.sub_domain, str) else as_tuple(self.sub_domain)
         sub_d = [s if isinstance(s, str) else as_tuple(s) for s in sub_d]

firedrake/embedding.py

@@ -4,7 +4,7 @@
 import ufl
 
 
-def get_embedding_dg_element(element):
+def get_embedding_dg_element(element, broken_cg=False):
     cell = element.cell
     degree = element.degree()
     family = lambda c: "DG" if c.is_simplex() else "DQ"
@@ -16,6 +16,8 @@ def get_embedding_dg_element(element):
                                                               for (c, d) in zip(cell.sub_cells(), degree)))
     else:
         scalar_element = finat.ufl.FiniteElement(family(cell), cell=cell, degree=degree)
+    if broken_cg:
+        scalar_element = finat.ufl.BrokenElement(scalar_element.reconstruct(family="Lagrange"))
     shape = element.value_shape
     if len(shape) == 0:
         DG = scalar_element

firedrake/mg/embedded.py

@@ -12,15 +12,22 @@
 __all__ = ("TransferManager", )
 
 
-native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ"])
-alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced Hsieh-Clough-Tocher", "Johnson-Mercier"])
-non_native_variants = frozenset(["integral", "fdm"])
+native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement"])
+alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced-Hsieh-Clough-Tocher", "Johnson-Mercier",
+                             "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu",
+                             "Guzman-Neilan", "Guzman-Neilan Bubble"])
+non_native_variants = frozenset(["integral", "fdm", "alfeld"])
 
 
 def get_embedding_element(element):
-    dg_element = get_embedding_dg_element(element)
-    if element.family() in alfeld_families:
-        dg_element = dg_element.reconstruct(variant="alfeld")
+    broken_cg = element.sobolev_space in {ufl.H1, ufl.H2}
+    dg_element = get_embedding_dg_element(element, broken_cg=broken_cg)
+    variant = element.variant() or "default"
+    family = element.family()
+    # Elements on Alfeld splits are embedded onto DG Powell-Sabin.
+    # This yields supermesh projection
+    if (family in alfeld_families) or ("alfeld" in variant.lower() and family != "Discontinuous Lagrange"):
+        dg_element = dg_element.reconstruct(variant="powell-sabin")
     return dg_element
 
 
@@ -67,7 +74,7 @@ def is_native(self, element):
             return True
         if isinstance(element.cell, ufl.TensorProductCell) and len(element.sub_elements) > 0:
             return reduce(and_, map(self.is_native, element.sub_elements))
-        return element.family() in native_families and not element.variant() in non_native_variants
+        return (element.family() in native_families) and not (element.variant() in non_native_variants)
 
     def _native_transfer(self, element, op):
         try:

tests/macro/test_macro_interp_project.py

@@ -1,34 +1,41 @@
 import pytest
+import numpy
 from firedrake import *
 
 
 def interp(u, f):
     u.interpolate(f)
+    return assemble(inner(u-f, u-f)*dx)**0.5
 
 
-def proj(u, f):
-    u.project(f)
+def proj(u, f, bcs=None):
+    u.project(f, bcs=bcs)
+    return assemble(inner(u-f, u-f)*dx)**0.5
 
 
 def proj_bc(u, f):
-    u.project(f, bcs=DirichletBC(u.function_space(), f, "on_boundary"))
+    return proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary"))
 
 
 def h1_proj(u, f, bcs=None):
     # compute h1 projection of f into u's function
     # space, store the result in u.
     v = TestFunction(u.function_space())
-    F = (inner(grad(u-f), grad(v)) * dx
+    d = {H2: grad, H1: grad, HCurl: curl, HDiv: div, HDivDiv: div}[u.ufl_element().sobolev_space]
+    F = (inner(d(u-f), d(v)) * dx
          + inner(u-f, v) * dx)
+    fcp = {"mode": "vanilla"}
     solve(F == 0, u,
           bcs=bcs,
           solver_parameters={"snes_type": "ksponly",
                              "ksp_type": "preonly",
-                             "pc_type": "cholesky"})
+                             "pc_type": "cholesky"},
+          form_compiler_parameters=fcp)
+    return assemble(F(u-f), form_compiler_parameters=fcp)**0.5
 
 
 def h1_proj_bc(u, f):
-    h1_proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary"))
+    return h1_proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary"))
 
 
 @pytest.fixture(params=("square", "cube"))
@@ -53,21 +60,120 @@ def test_projection_scalar_monomial(op, mesh, degree, variant):
     u = Function(V)
     x = SpatialCoordinate(mesh)
     f = sum(x) ** degree
-    op(u, f)
-    error = sqrt(assemble(inner(u - f, u - f) * dx))
+    error = op(u, f)
     assert error < 1E-7
 
 
-@pytest.mark.parametrize(('el', 'accorder'),
-                         [('HCT', 3),
-                          ('HCT-red', 2)])
+@pytest.fixture
+def hierarchy(request):
+    refine = 1
+    mh2 = MeshHierarchy(UnitSquareMesh(5, 5), refine)
+    mh3 = MeshHierarchy(UnitCubeMesh(3, 3, 3), refine)
+    return {2: mh2, 3: mh3}
+
+
+def mesh_sizes(mh):
+    mesh_size = []
+    for msh in mh:
+        DG0 = FunctionSpace(msh, "DG", 0)
+        h = Function(DG0).interpolate(CellDiameter(msh))
+        with h.dat.vec as hvec:
+            _, maxh = hvec.max()
+        mesh_size.append(maxh)
+    return mesh_size
+
+
+def conv_rates(x, h):
+    x = numpy.asarray(x)
+    h = numpy.asarray(h)
+    return numpy.log2(x[:-1] / x[1:]) / numpy.log2(h[:-1] / h[1:])
+
+
+def run_convergence(mh, el, deg, convrate, op):
+    errors = []
+    for msh in mh:
+        V = FunctionSpace(msh, el, deg)
+        u = Function(V)
+        x = SpatialCoordinate(msh)
+        f = sum(x) ** (deg+1)
+        if u.ufl_shape != ():
+            f = f * Constant(numpy.ones(u.ufl_shape))
+        errors.append(op(u, f))
+
+    conv = conv_rates(errors, mesh_sizes(mh))
+    assert numpy.all(conv > convrate - 0.25)
+
+
+# Test L2/H1 convergence on C1 elements
+@pytest.mark.parametrize(('dim', 'el', 'deg', 'convrate'),
+                         [(2, 'PS6', 2, 2),
+                          (2, 'PS12', 2, 2),
+                          (2, 'HCT', 3, 3),
+                          (2, 'HCT-red', 3, 2),
+                          ])
 @pytest.mark.parametrize('op', (proj, h1_proj))
-def test_projection_hct(el, accorder, op):
-    msh = UnitSquareMesh(1, 1)
-    V = FunctionSpace(msh, el, 3)
+def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op):
+    if op == proj:
+        convrate += 1
+    run_convergence(hierarchy[dim], el, deg, convrate, op)
+
+
+# Test L2/H1 convergence on Stokes elements
+@pytest.mark.parametrize(('dim', 'el', 'deg', 'convrate'),
+                         [(2, 'Alfeld-Sorokina', 2, 2),
+                          (3, 'Alfeld-Sorokina', 2, 2),
+                          (2, 'Reduced-Arnold-Qin', 2, 1),
+                          (3, 'Guzman-Neilan', 3, 1),
+                          (3, 'Christiansen-Hu', 1, 1),
+                          (2, 'Bernardi-Raugel', 2, 1),
+                          (3, 'Bernardi-Raugel', 3, 1),
+                          (2, 'Johnson-Mercier', 1, 1),
+                          (3, 'Johnson-Mercier', 1, 1),
+                          ])
+@pytest.mark.parametrize('op', (proj, h1_proj))
+def test_piola_convergence(hierarchy, dim, el, deg, convrate, op):
+    if op == proj:
+        convrate += 1
+    run_convergence(hierarchy[dim], el, deg, convrate, op)
+
+
+# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions
+@pytest.mark.parametrize('enriched', (False, True))
+def test_supersmooth_bcs(mesh, enriched):
+    tdim = mesh.topological_dimension()
+    if enriched:
+        cell = mesh.ufl_cell()
+        element = EnrichedElement(FiniteElement("AS", cell=cell, degree=2),
+                                  FiniteElement("GNB", cell=cell, degree=tdim))
+        V = FunctionSpace(mesh, element)
+    else:
+        V = FunctionSpace(mesh, "Alfeld-Sorokina", 2)
+
+    # check that V in H1
+    assert V.ufl_element().sobolev_space == H1
+
+    # check that V is supersmooth
+    nodes = V.finat_element.fiat_equivalent.dual.nodes
+    deriv_nodes = [i for i, node in enumerate(nodes) if len(node.deriv_dict)]
+    assert len(deriv_nodes) == tdim + 1
+
+    deriv_ids = V.cell_node_list[:, deriv_nodes]
     u = Function(V)
-    x = SpatialCoordinate(msh)
-    f = sum(x) ** accorder
-    op(u, f)
-    error = sqrt(assemble(inner(u - f, u - f) * dx))
-    assert error < 1E-9
+
+    CG = FunctionSpace(mesh, "Lagrange", 2)
+    RT = FunctionSpace(mesh, "RT", 1)
+    for sub in [1, (1, 2), "on_boundary"]:
+        bc = DirichletBC(V, 0, sub)
+
+        # check that we have the expected number of bc nodes
+        nnodes = len(bc.nodes)
+        expected = tdim * len(DirichletBC(CG, 0, sub).nodes)
+        if enriched:
+            expected += len(DirichletBC(RT, 0, sub).nodes)
+        assert nnodes == expected
+
+        # check that the bc does not set the derivative nodes
+        u.assign(111)
+        u.dat.data_wo[deriv_ids] = 42
+        bc.zero(u)
+        assert numpy.allclose(u.dat.data_ro[deriv_ids], 42)

tests/macro/test_macro_multigrid.py

@@ -152,7 +152,12 @@ def test_macro_multigrid_poisson(hierarchy, degree, variant):
     uh = Function(V)
     problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
     solver = LinearVariationalSolver(problem, solver_parameters=mg_params)
-    solver.solve()
+    if complex_mode and variant == "alfeld":
+        with pytest.raises(NotImplementedError):
+            solver.solve()
+    else:
+        solver.solve()
+
     expected = 10
     if mesh.geometric_dimension() == 3 and variant == "alfeld":
         expected = 14

tests/macro/test_macro_quadrature.py

@@ -7,7 +7,7 @@
 def alfeld_split(msh):
     dim = msh.geometric_dimension()
     coords = msh.coordinates.dat.data.reshape((-1, dim))
-    coords = numpy.row_stack((coords, numpy.average(coords, 0)))
+    coords = numpy.vstack((coords, numpy.average(coords, 0)))
     cells = [list(map(lambda i: dim+1 if i == j else i, range(dim+1))) for j in range(dim+1)]
     plex = plex_from_cell_list(dim, cells, coords, msh.comm)
     return Mesh(plex)
@@ -68,8 +68,8 @@ def test_macro_quadrature_piecewise(degree, variant, meshes):
             c = Constant(numpy.arange(1, gdim+1))
             expr = abs(dot(c, x - x0)) ** degree
         elif variant == "iso":
-            vecs = list(map(Constant, numpy.row_stack([numpy.eye(gdim),
-                                                       numpy.ones((max(degree-gdim, 0), gdim))])))
+            vecs = list(map(Constant, numpy.vstack([numpy.eye(gdim),
+                                                    numpy.ones((max(degree-gdim, 0), gdim))])))
             expr = numpy.prod([abs(dot(c, x)) for c in vecs[:degree]])
         else:
             raise ValueError("Unexpected variant")

tests/macro/test_macro_solve.py

@@ -33,23 +33,35 @@ def mixed_element(mh, variant):
     return Vel, Pel
 
 
-def conv_rates(x):
+def mesh_sizes(mh):
+    mesh_size = []
+    for msh in mh:
+        DG0 = FunctionSpace(msh, "DG", 0)
+        h = Function(DG0).interpolate(CellDiameter(msh))
+        with h.dat.vec as hvec:
+            _, maxh = hvec.max()
+        mesh_size.append(maxh)
+    return mesh_size
+
+
+def conv_rates(x, h):
     x = np.asarray(x)
-    return np.log2(x[:-1] / x[1:])
+    h = np.asarray(h)
+    return np.log2(x[:-1] / x[1:]) / np.log2(h[:-1] / h[1:])
 
 
 @pytest.fixture
 def convergence_test(variant):
     if variant == "iso":
-        def check(uerr, perr):
-            return (conv_rates(uerr)[-1] >= 1.9
+        def check(uerr, perr, h):
+            return (conv_rates(uerr, h)[-1] >= 1.9
                     and np.allclose(perr, 0, atol=1.e-8))
     elif variant == "alfeld":
-        def check(uerr, perr):
+        def check(uerr, perr, h):
             return (np.allclose(uerr, 0, atol=5.e-9)
                     and np.allclose(perr, 0, atol=5.e-7))
     elif variant == "th":
-        def check(uerr, perr):
+        def check(uerr, perr, h):
             return (np.allclose(uerr, 0, atol=1.e-10)
                     and np.allclose(perr, 0, atol=1.e-8))
     return check
@@ -100,7 +112,7 @@ def test_riesz(mh, variant, mixed_element, convergence_test):
         u_err.append(errornorm(as_vector(zexact[:dim]), uh))
         p_err.append(errornorm(zexact[-1], ph))
 
-    assert convergence_test(u_err, p_err)
+    assert convergence_test(u_err, p_err, mesh_sizes(mh))
 
 
 def stokes_mms(Z, zexact):
@@ -128,7 +140,6 @@ def test_stokes(mh, variant, mixed_element, convergence_test):
     dim = mh[0].geometric_dimension()
     if variant == "iso" and dim == 3:
         pytest.xfail("P2:P1 iso x P1 is not inf-sup stable in 3D")
-
     u_err = []
     p_err = []
     el1, el2 = mixed_element
@@ -155,7 +166,7 @@ def test_stokes(mh, variant, mixed_element, convergence_test):
         u_err.append(errornorm(as_vector(zexact[:dim]), uh))
         p_err.append(errornormL2_0(zexact[-1], ph))
 
-    assert convergence_test(u_err, p_err)
+    assert convergence_test(u_err, p_err, mesh_sizes(mh))
 
 
 def test_div_free(mh, variant, mixed_element, div_test):

tests/regression/test_helmholtz_zany.py

@@ -54,11 +54,7 @@ def helmholtz(n, el_type, degree, perturb):
 # It is, somewhat oddly, a suitable C^1 nonconforming element but
 # not a suitable C^0 nonconforming one.
 @pytest.mark.parametrize(('el', 'deg', 'convrate'),
-                         [('PS6', 2, 2),
-                          ('PS12', 2, 2),
-                          ('Hermite', 3, 3),
-                          ('HCT', 3, 3),
-                          ('HCT', 4, 4),
+                         [('Hermite', 3, 3),
                           ('Bell', 5, 4),
                           ('Argyris', 5, 5),
                           ('Argyris', 6, 6)])