Test the Johnson-Mercier macroelement (#3548)

---------

Co-authored-by: Robert Kirby <robert.c.kirby@gmail.com>

tests/regression/test_projection_zany.py

@@ -13,33 +13,42 @@
 from firedrake import *
 
 
-def do_projection(n, el_type, degree):
-    # Create mesh and define function space
-    mesh = UnitSquareMesh(2**n, 2**n)
+@pytest.fixture
+def hierarchy(request):
+    msh = UnitSquareMesh(2, 2)
+    return MeshHierarchy(msh, 2)
 
+
+def do_projection(mesh, el_type, degree):
     V = FunctionSpace(mesh, el_type, degree)
 
     # Define variational problem
     u = TrialFunction(V)
     v = TestFunction(V)
-    x, y = SpatialCoordinate(mesh)
-    f = sin(x*pi)*sin(2*pi*y)
+    x = SpatialCoordinate(mesh)
+
+    f = np.prod([sin((1+i) * x[i]*pi) for i in range(len(x))])
+    f = f * Constant(np.ones(u.ufl_shape))
+
     a = inner(u, v) * dx
     L = inner(f, v) * dx
 
     # Compute solution
-    x = Function(V)
-    solve(a == L, x, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'})
+    fh = Function(V)
+    solve(a == L, fh, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'})
 
-    return sqrt(assemble(inner(x - f, x - f) * dx))
+    return sqrt(assemble(inner(fh - f, fh - f) * dx))
 
 
 @pytest.mark.parametrize(('el', 'deg', 'convrate'),
-                         [('Morley', 2, 2.4),
+                         [('Johnson-Mercier', 1, 1.8),
+                          ('Morley', 2, 2.4),
+                          ('HCT-red', 3, 2.7),
+                          ('HCT', 3, 3),
                           ('Hermite', 3, 3),
                           ('Bell', 5, 4),
                           ('Argyris', 5, 4.9)])
-def test_firedrake_projection_scalar_convergence(el, deg, convrate):
-    diff = np.array([do_projection(i, el, deg) for i in range(1, 4)])
+def test_projection_zany_convergence_2d(hierarchy, el, deg, convrate):
+    diff = np.array([do_projection(m, el, deg) for m in hierarchy])
     conv = np.log2(diff[:-1] / diff[1:])
     assert (np.array(conv) > convrate).all()

tests/regression/test_aw.py → tests/regression/test_stress_els.py

@@ -7,12 +7,14 @@
 convergence_orders = lambda x: -np.log2(relative_magnitudes(x))
 
 
-@pytest.fixture(scope='module', params=["conforming", "nonconforming"])
+@pytest.fixture(scope='module', params=["conforming", "nonconforming", "macro"])
 def stress_element(request):
     if request.param == "conforming":
         return FiniteElement("AWc", triangle, 3)
     elif request.param == "nonconforming":
         return FiniteElement("AWnc", triangle, 2)
+    elif request.param == "macro":
+        return FiniteElement("JM", triangle, 1)
     else:
         raise ValueError("Unknown family")
 
@@ -25,7 +27,7 @@ def mesh_hierarchy(request):
     return mh
 
 
-def test_aw_convergence(stress_element, mesh_hierarchy):
+def test_stress_displacement_convergence(stress_element, mesh_hierarchy):
 
     mesh = mesh_hierarchy[0]
     V = FunctionSpace(mesh, mesh.coordinates.ufl_element())
@@ -95,7 +97,7 @@ def epsilon(u):
             )  # noqa: E123
 
         # Augmented Lagrangian preconditioner
-        gamma = Constant(1E2)
+        gamma = Constant(1E4)
         Jp = inner(A(sig), tau)*dx + inner(div(sig) * gamma, div(tau))*dx + inner(u * (1/gamma), v) * dx
 
         params = {"mat_type": "matfree",
@@ -132,11 +134,15 @@ def epsilon(u):
         assert min(convergence_orders(l2_u)) > 1.9
         assert min(convergence_orders(l2_sigma)) > 1
         assert min(convergence_orders(l2_div_sigma)) > 1.9
+    elif stress_element.family().startswith("Johnson-Mercier"):
+        assert min(convergence_orders(l2_u)) > 1.9
+        assert min(convergence_orders(l2_sigma)) > 1
+        assert min(convergence_orders(l2_div_sigma)) > 0.9
     else:
         raise ValueError("Don't know what the convergence should be")
 
 
-def test_aw_conditioning(stress_element, mesh_hierarchy):
+def test_mass_conditioning(stress_element, mesh_hierarchy):
     mass_cond = []
     for msh in mesh_hierarchy[:3]:
         Sig = FunctionSpace(msh, stress_element)