Test high-order Argyris and HCT (#3658)

* Test high-order Argyris and HCT

tests/regression/test_helmholtz_zany.py

@@ -20,7 +20,7 @@
 def helmholtz(n, el_type, degree, perturb):
     mesh = UnitSquareMesh(2**n, 2**n)
     if perturb:
-        V = FunctionSpace(mesh, mesh.coordinates.ufl_element())
+        V = mesh.coordinates.function_space()
         eps = Constant(1 / 2**(n+1))
 
         x, y = SpatialCoordinate(mesh)
@@ -30,38 +30,40 @@ def helmholtz(n, el_type, degree, perturb):
 
     V = FunctionSpace(mesh, el_type, degree)
     x = SpatialCoordinate(V.mesh())
+    degree = V.ufl_element().degree()
 
     # Define variational problem
-    lmbda = 1
+    k = Constant(3)
+    u_exact = cos(x[0]*pi*k) * cos(x[1]*pi*k)
+
+    lmbda = Constant(1)
     u = TrialFunction(V)
     v = TestFunction(V)
-    f = Function(V)
-    f.project((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2))
-    a = (inner(grad(u), grad(v)) + lmbda * inner(u, v)) * dx
-    L = inner(f, v) * dx
+    a = inner(grad(u), grad(v))*dx + lmbda * inner(u, v)*dx
+    L = a(v, u_exact)
 
     # Compute solution
-    assemble(a)
-    assemble(L)
     sol = Function(V)
-    solve(a == L, sol, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'})
+    solve(a == L, sol)
 
-    # Analytical solution
-    f.project(cos(x[0]*pi*2)*cos(x[1]*pi*2))
-    return sqrt(assemble(inner(sol - f, sol - f) * dx))
+    # Return the H1 norm
+    return sqrt(assemble(a(sol - u_exact, sol - u_exact)))
 
 
-# Test convergence on Hermite, Bell, and Argyris
+# Test H1 convergence on Hermite, HCT, Bell, and Argyris
 # Morley is omitted since it only can be used on 4th-order problems.
 # It is, somewhat oddly, a suitable C^1 nonconforming element but
 # not a suitable C^0 nonconforming one.
 @pytest.mark.parametrize(('el', 'deg', 'convrate'),
-                         [('Hermite', 3, 3.8),
-                          ('HCT', 3, 3.8),
-                          ('Bell', 5, 4.8),
-                          ('Argyris', 5, 4.8)])
+                         [('Hermite', 3, 3),
+                          ('HCT', 3, 3),
+                          ('HCT', 4, 4),
+                          ('Bell', 5, 4),
+                          ('Argyris', 5, 5),
+                          ('Argyris', 6, 6)])
 @pytest.mark.parametrize("perturb", [False, True], ids=["Regular", "Perturbed"])
 def test_firedrake_helmholtz_scalar_convergence(el, deg, convrate, perturb):
-    diff = np.array([helmholtz(i, el, deg, perturb) for i in range(1, 4)])
+    l = 4
+    diff = np.array([helmholtz(i, el, deg, perturb) for i in range(l, l+2)])
     conv = np.log2(diff[:-1] / diff[1:])
-    assert (np.array(conv) > convrate).all()
+    assert (np.array(conv) > convrate - 0.3).all()

tests/regression/test_projection_zany.py

@@ -2,7 +2,7 @@
 
 on the unit square of a function
 
-  f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi)
+  f(x, y) = LegendreP(degree + 2, x + y)
 
 using elements with nonstandard pullbacks
 """
@@ -13,10 +13,22 @@
 from firedrake import *
 
 
+relative_magnitudes = lambda x: np.array(x)[1:] / np.array(x)[:-1]
+convergence_orders = lambda x: -np.log2(relative_magnitudes(x))
+
+
+def jrc(a, b, n):
+    """Jacobi recurrence coefficients"""
+    an = (2*n+1+a+b)*(2*n+2+a+b) / (2*(n+1)*(n+1+a+b))
+    bn = (a+b)*(a-b)*(2*n+1+a+b) / (2*(n+1)*(n+1+a+b)*(2*n+a+b))
+    cn = (n+a)*(n+b)*(2*n+2+a+b) / ((n+1)*(n+1+a+b)*(2*n+a+b))
+    return an, bn, cn
+
+
 @pytest.fixture
 def hierarchy(request):
     msh = UnitSquareMesh(2, 2)
-    return MeshHierarchy(msh, 2)
+    return MeshHierarchy(msh, 4)
 
 
 def do_projection(mesh, el_type, degree):
@@ -27,28 +39,59 @@ def do_projection(mesh, el_type, degree):
     v = TestFunction(V)
     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))
-
+    m = degree + 2
+    z = sum(x)
+    p0 = 1
+    p1 = z
+    for n in range(1, m):
+        an, bn, cn = jrc(0, 0, n)
+        ptemp = p1
+        p1 = (an * z + bn) * p1 - cn * p0
+        p0 = ptemp
+
+    f = p1 * Constant(np.ones(u.ufl_shape))
     a = inner(u, v) * dx
     L = inner(f, v) * dx
 
     # Compute solution
     fh = Function(V)
-    solve(a == L, fh, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'})
-
+    solve(a == L, fh)
     return sqrt(assemble(inner(fh - f, fh - f) * dx))
 
 
 @pytest.mark.parametrize(('el', 'deg', 'convrate'),
                          [('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)])
+                          ('HCT', 3, 3.7),
+                          ('HCT', 4, 4.8),
+                          ('Hermite', 3, 3.8),
+                          ('Bell', 5, 4.7),
+                          ('Argyris', 5, 5.8),
+                          ('Argyris', 6, 6.7)])
 def test_projection_zany_convergence_2d(hierarchy, el, deg, convrate):
-    diff = np.array([do_projection(m, el, deg) for m in hierarchy])
+    diff = np.array([do_projection(m, el, deg) for m in hierarchy[2:]])
     conv = np.log2(diff[:-1] / diff[1:])
     assert (np.array(conv) > convrate).all()
+
+
+@pytest.mark.parametrize(('element', 'degree'),
+                         [('HCT-red', 3),
+                          ('HCT', 3),
+                          ('HCT', 4),
+                          ('Argyris', 5),
+                          ('Argyris', 6)])
+def test_mass_conditioning(element, degree, hierarchy):
+    mass_cond = []
+    for msh in hierarchy[1:4]:
+        V = FunctionSpace(msh, element, degree)
+        u = TrialFunction(V)
+        v = TestFunction(V)
+        a = inner(u, v)*dx
+        B = assemble(a, mat_type="aij").M.handle
+        A = B.convert("dense").getDenseArray()
+        kappa = np.linalg.cond(A)
+
+        mass_cond.append(kappa)
+
+    assert max(relative_magnitudes(mass_cond)) < 1.1