Test Stokes macroelements (#3795)

* Test Stokes macroelements

* Test GuzmanNeilanSecondKind

* get_embedding_dg_element with superdegree

firedrake/embedding.py

@@ -6,15 +6,16 @@
 
 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"
     if isinstance(cell, ufl.TensorProductCell):
+        degree = element.degree()
         if type(degree) is int:
             scalar_element = finat.ufl.FiniteElement("DQ", cell=cell, degree=degree)
         else:
             scalar_element = finat.ufl.TensorProductElement(*(finat.ufl.FiniteElement(family(c), cell=c, degree=d)
                                                               for (c, d) in zip(cell.sub_cells(), degree)))
     else:
+        degree = element.embedded_superdegree
         scalar_element = finat.ufl.FiniteElement(family(cell), cell=cell, degree=degree)
     if broken_cg:
         scalar_element = finat.ufl.BrokenElement(scalar_element.reconstruct(family="Lagrange"))

tests/macro/test_macro_interp_project.py

@@ -21,6 +21,7 @@ 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())
+
     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)
@@ -123,57 +124,16 @@ def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op):
                          [(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, 'Bernardi-Raugel', 1, 1),
+                          (3, 'Bernardi-Raugel', 1, 1),
                           (2, 'Johnson-Mercier', 1, 1),
                           (3, 'Johnson-Mercier', 1, 1),
+                          (3, 'Guzman-Neilan 1st kind H1', 1, 1),
+                          (3, 'Guzman-Neilan H1(div)', 3, 2),
                           ])
 @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)
-
-    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_stokes_macroelements.py

@@ -0,0 +1,86 @@
+import pytest
+import numpy
+from firedrake import *
+
+
+@pytest.fixture(params=("square", "cube"))
+def mesh(request):
+    n = 2
+    if request.param == "square":
+        return UnitSquareMesh(n, n)
+    elif request.param == "cube":
+        return UnitCubeMesh(n, n, n)
+
+
+@pytest.fixture(params=("SV", "GN", "GN2", "GNH1div"))
+def space(request, mesh):
+    family = request.param
+    dim = mesh.topological_dimension()
+    if family == "GN":
+        V = FunctionSpace(mesh, "GN", 1)
+        Q = FunctionSpace(mesh, "DG", 0)
+    elif family == "GN2":
+        V = FunctionSpace(mesh, "GN2", 1)
+        Q = FunctionSpace(mesh, "DG", 0, variant="alfeld")
+    elif family == "GNH1div":
+        V = FunctionSpace(mesh, "GNH1div", dim)
+        Q = FunctionSpace(mesh, "CG", 1, variant="alfeld")
+    elif family == "SV":
+        V = VectorFunctionSpace(mesh, "CG", dim, variant="alfeld")
+        Q = FunctionSpace(mesh, "DG", dim-1, variant="alfeld")
+    else:
+        raise ValueError("Unexpected family")
+    return V * Q
+
+
+# Test that div(V) is contained in Q
+def test_stokes_complex(mesh, space):
+    Z = space
+    z = Function(Z)
+    u, p = z.subfunctions
+
+    for k in range(len(u.dat.data)):
+        u.dat.data_wo[k] = 1
+        p.interpolate(div(u))
+        assert norm(div(u) - p) < 1E-10
+        u.dat.data_wo[k] = 0
+
+
+# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions
+def test_supersmooth_bcs(mesh):
+    tdim = mesh.topological_dimension()
+    if tdim == 3:
+        V = FunctionSpace(mesh, "GNH1div", 3)
+    else:
+        V = FunctionSpace(mesh, "Alfeld-Sorokina", 2)
+
+    assert V.finat_element.complex.is_macrocell()
+
+    # 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)
+
+    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 tdim == 3:
+            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/multigrid/test_p_multigrid.py

@@ -116,7 +116,7 @@ def test_prolong_low_order_to_restricted(tp_mesh, tp_family, variant):
         P = prolongation_matrix_matfree(uc, v).getPythonContext()
         P._prolong()
 
-    assert norm(ui + uf - uc, "L2") < 2E-14
+    assert norm(ui + uf - uc, "L2") < 1E-13
 
 
 @pytest.fixture(params=["triangles", "quadrilaterals"], scope="module")