From 8d4700c823a6170ef6b5b495ee720a219f1e896e Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 18 Feb 2019 16:02:13 +0000
Subject: [PATCH 01/43] DROP BEFORE MERGE

---
 Jenkinsfile | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Jenkinsfile b/Jenkinsfile
index 7684a36f23..263ecdf503 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -121,4 +121,3 @@ sudo docker push firedrakeproject/firedrake:latest
     }
   }
 }
-

From e7a3047d3603fcab950e8eb350eb62688e7dc92c Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Tue, 5 Mar 2019 22:15:20 +0000
Subject: [PATCH 02/43] tests for dpc element

---
 tests/regression/test_dpc.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 tests/regression/test_dpc.py

diff --git a/tests/regression/test_dpc.py b/tests/regression/test_dpc.py
new file mode 100644
index 0000000000..783d66e4a2
--- /dev/null
+++ b/tests/regression/test_dpc.py
@@ -0,0 +1,24 @@
+import pytest
+import numpy as np
+from firedrake import *
+
+
+def test_simple():
+    mesh = UnitSquareMesh(10, 10, quadrilateral=True)
+    V = FunctionSpace(mesh, "DPC", 3)
+    x = SpatialCoordinate(mesh)
+    u = project(x[0]**3 + x[1]**3, V)
+    assert (u.dat.data[0]) < 1e-14
+
+
+def test_enrichment():
+    mesh = UnitSquareMesh(10, 10, quadrilateral=True)
+    x = SpatialCoordinate(mesh)
+    dPc = FiniteElement("DPC", "quadrilateral", 3)
+    V = FunctionSpace(mesh, "DPC", 3)
+    W = FunctionSpace(mesh, "CG", 3)
+    u = project(x[0]**3 + x[1]**3, V)
+    exact = Function(W)
+    exact.interpolate(x[0]**3 + x[1]**3)
+    # make sure that these are the same
+    assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14

From d443661712a426fed17938e43b47672bf994cb91 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 29 Mar 2019 12:16:18 +0000
Subject: [PATCH 03/43] modified to test dPc element

---
 tests/regression/test_steady_advection_2D.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/tests/regression/test_steady_advection_2D.py b/tests/regression/test_steady_advection_2D.py
index 2827d9689e..1d574529a5 100644
--- a/tests/regression/test_steady_advection_2D.py
+++ b/tests/regression/test_steady_advection_2D.py
@@ -14,10 +14,12 @@ def mesh(request):
     return UnitSquareMesh(5, 5, quadrilateral=request.param)
 
 
-@pytest.fixture
-def DG0(mesh):
-    return FunctionSpace(mesh, "DG", 0)
-
+@pytest.fixture(params=["DG", "DPC"])
+def DG0(request, mesh):
+    if mesh.ufl_cell() == triangle:
+        return FunctionSpace(mesh, "DG", 0)
+    else:
+        return FunctionSpace(mesh, request.param, 0)
 
 @pytest.fixture
 def DG1(mesh):
@@ -67,7 +69,7 @@ def test_left_to_right(mesh, DG0, W):
     run_left_to_right(mesh, DG0, W)
 
 
-@pytest.mark.parallel
+#@pytest.mark.parallel
 def test_left_to_right_parallel(mesh, DG0, W):
     run_left_to_right(mesh, DG0, W)
 
@@ -104,6 +106,6 @@ def test_up_to_down(mesh, DG1, W):
     run_up_to_down(mesh, DG1, W)
 
 
-@pytest.mark.parallel
+#@pytest.mark.parallel
 def test_up_to_down_parallel(mesh, DG1, W):
     run_up_to_down(mesh, DG1, W)

From 26cbe9118901cb496b06b48bf48ba8f112cc67c1 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 12 Apr 2019 18:51:06 +0100
Subject: [PATCH 04/43] new tests for serendipity

---
 tests/regression/test_dpc.py                  |  24 -----
 tests/regression/test_helmholtz.py            |  17 ++-
 .../regression/test_helmholtz_serendipity.py  | 102 ++++++++++++++++++
 3 files changed, 117 insertions(+), 26 deletions(-)
 delete mode 100644 tests/regression/test_dpc.py
 create mode 100644 tests/regression/test_helmholtz_serendipity.py

diff --git a/tests/regression/test_dpc.py b/tests/regression/test_dpc.py
deleted file mode 100644
index 783d66e4a2..0000000000
--- a/tests/regression/test_dpc.py
+++ /dev/null
@@ -1,24 +0,0 @@
-import pytest
-import numpy as np
-from firedrake import *
-
-
-def test_simple():
-    mesh = UnitSquareMesh(10, 10, quadrilateral=True)
-    V = FunctionSpace(mesh, "DPC", 3)
-    x = SpatialCoordinate(mesh)
-    u = project(x[0]**3 + x[1]**3, V)
-    assert (u.dat.data[0]) < 1e-14
-
-
-def test_enrichment():
-    mesh = UnitSquareMesh(10, 10, quadrilateral=True)
-    x = SpatialCoordinate(mesh)
-    dPc = FiniteElement("DPC", "quadrilateral", 3)
-    V = FunctionSpace(mesh, "DPC", 3)
-    W = FunctionSpace(mesh, "CG", 3)
-    u = project(x[0]**3 + x[1]**3, V)
-    exact = Function(W)
-    exact.interpolate(x[0]**3 + x[1]**3)
-    # make sure that these are the same
-    assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14
diff --git a/tests/regression/test_helmholtz.py b/tests/regression/test_helmholtz.py
index 75f1a692d9..25a49d4fbd 100644
--- a/tests/regression/test_helmholtz.py
+++ b/tests/regression/test_helmholtz.py
@@ -20,11 +20,11 @@
 cwd = abspath(dirname(__file__))
 
 
-def helmholtz(r, quadrilateral=False, degree=2, mesh=None):
+def helmholtz(r, quadrilateral=False, degree=2, mesh=None, element="CG"):
     # Create mesh and define function space
     if mesh is None:
         mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral)
-    V = FunctionSpace(mesh, "CG", degree)
+    V = FunctionSpace(mesh, element, degree)
 
     # Define variational problem
     lmbda = 1
@@ -77,6 +77,19 @@ def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals(testcase, conv
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
 
 
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [((1, (4, 6)), 1.9),
+                          ((2, (3, 6)), 2.9),
+                          ((3, (2, 4)), 3.9),
+                          ((4, (2, 4)), 4.7)])
+def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_dpc(testcase, convrate):
+    degree, (start, end) = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = helmholtz(ii, quadrilateral=True, degree=degree, element="S")[0]
+    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
+
+
 def run_firedrake_helmholtz_on_quadrilateral_mesh_from_file():
     meshfile = join(cwd, "..", "meshes", "unitsquare_unstructured_quadrilaterals.msh")
     assert helmholtz(None, mesh=Mesh(meshfile))[0] <= 0.01
diff --git a/tests/regression/test_helmholtz_serendipity.py b/tests/regression/test_helmholtz_serendipity.py
new file mode 100644
index 0000000000..2d1dafa42d
--- /dev/null
+++ b/tests/regression/test_helmholtz_serendipity.py
@@ -0,0 +1,102 @@
+"""This demo program solves Helmholtz's equation
+
+  - div grad u(x, y) + u(x,y) = f(x, y)
+
+on the unit square with source f given by
+
+  f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi)
+
+and the analytical solution
+
+  u(x, y) = cos(x[0]*2*pi)*cos(x[1]*2*pi)
+"""
+
+from os.path import abspath, dirname, join
+import numpy as np
+import pytest
+
+from firedrake import *
+
+cwd = abspath(dirname(__file__))
+
+
+def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
+    # Create mesh and define function space
+    if mesh is None:
+        mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral)
+    V = FunctionSpace(mesh, "S", degree)
+
+    x, y = SpatialCoordinate(mesh)
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    #uex = cos(pi*x)*cos(pi*y)
+    uex = cos(x*pi*2)*cos(y*pi*2)
+    f = -div(grad(uex)) + uex
+
+    a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=8)
+    L = inner(f, v)*dx(degree=8)
+
+    params = {"snes_type": "ksponly",
+              "ksp_type": "preonly",
+              "pc_type": "lu",
+              "snes_rtol": 1e-16,
+              "snes_atol": 1e-25}
+
+    # Compute solution
+    #assemble(a)
+    #assemble(L)
+    sol = Function(V)
+    solve(a == L, sol, solver_parameters=params)
+
+    # Analytical solution
+    f = Function(V)
+    f.project(cos(x*pi*2)*cos(y*pi*2))
+    return sqrt(assemble(dot(sol - uex, sol - uex) * dx)), sol, uex
+
+
+def run_firedrake_helmholtz():
+    diff = np.array([helmholtz(i)[0] for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > 2.8).all()
+
+
+def test_firedrake_helmholtz_serial():
+    run_firedrake_helmholtz()
+
+
+#@pytest.mark.parallel
+#def test_firedrake_helmholtz_parallel():
+    #run_firedrake_helmholtz()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [((1, (4, 6)), 1.9),
+                          ((2, (3, 6)), 2.9),
+                          ((3, (2, 4)), 3.9),
+                          ((4, (2, 4)), 4.7)])
+def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_s(testcase, convrate):
+    degree, (start, end) = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = helmholtz(ii, quadrilateral=True, degree=degree)[0]
+    print(l2err)
+    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
+
+
+def run_firedrake_helmholtz_on_quadrilateral_mesh_from_file():
+    meshfile = join(cwd, "..", "meshes", "unitsquare_unstructured_quadrilaterals.msh")
+    assert helmholtz(None, mesh=Mesh(meshfile))[0] <= 0.01
+
+
+def test_firedrake_helmholtz_on_quadrilateral_mesh_from_file_serial():
+    run_firedrake_helmholtz_on_quadrilateral_mesh_from_file()
+
+
+#@pytest.mark.parallel
+#def test_firedrake_helmholtz_on_quadrilateral_mesh_from_file_parallel():
+    #run_firedrake_helmholtz_on_quadrilateral_mesh_from_file()

From 04f414528f5643b11f72a2d739d103d528cf9f7f Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sun, 28 Apr 2019 21:32:47 +0100
Subject: [PATCH 05/43] new tests for serendipity element

---
 AUTHORS                                       |  1 +
 .../test_helmholtz_serendipity_3d.py          | 44 +++++++++++++++++++
 .../regression/test_helmholtz_serendipity.py  | 10 ++---
 3 files changed, 50 insertions(+), 5 deletions(-)
 create mode 100644 tests/extrusion/test_helmholtz_serendipity_3d.py

diff --git a/AUTHORS b/AUTHORS
index 50b4a21eec..261853cce2 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -58,3 +58,4 @@ Jemma Shipton
 Tianjiao Sun
 Florian Wechsung
 Fangyi Zhou
+Cyrus Cheng
diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
new file mode 100644
index 0000000000..aba6f689fe
--- /dev/null
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -0,0 +1,44 @@
+"""Tests for scalar Helmholtz convergence on Serendipity elements on extruded meshes"""
+import numpy as np
+import pytest
+
+from firedrake import *
+
+
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [(("S", 1, (4, 6)), 1.9),
+                          (("S", 2, (3, 5)), 2.9),
+                          (("S", 3, (2, 4)), 3.9),
+                          (("S", 4, (2, 4)), 4.8),
+                          (("S", 5, (2, 5)), 5.7),
+                          (("S", 6, (2, 4)), 6.8),
+                          (("S", 7, (2, 4)), 7.7)])
+def test_scalar_convergence(extmesh, testcase, convrate):
+    family, degree, (start, end) = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        mesh = extmesh(2**ii, 2**ii, 2**ii, quadrilateral=True)
+
+        fspace = FunctionSpace(mesh, family, degree, vfamily=family, vdegree=degree)
+
+        u = TrialFunction(fspace)
+        v = TestFunction(fspace)
+
+        x, y, z = SpatialCoordinate(mesh)
+
+        uex = cos(2*np.pi*x)*cos(2*np.pi*y)*cos(2*np.pi*z)
+        f = -div(grad(uex)) + uex
+
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=16)
+        L = inner(f, v)*dx(degree=16)
+
+        params = {"snes_type": "ksponly",
+                  "ksp_type": "preonly",
+                  "pc_type": "lu"}
+
+        sol = Function(fspace)
+        solve(a == L, sol, solver_parameters=params)
+
+        l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
+    print(np.array(np.log2(l2err[:-1]/l2err[1:])))
+    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
diff --git a/tests/regression/test_helmholtz_serendipity.py b/tests/regression/test_helmholtz_serendipity.py
index 2d1dafa42d..baf3165700 100644
--- a/tests/regression/test_helmholtz_serendipity.py
+++ b/tests/regression/test_helmholtz_serendipity.py
@@ -36,8 +36,8 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
     uex = cos(x*pi*2)*cos(y*pi*2)
     f = -div(grad(uex)) + uex
 
-    a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=8)
-    L = inner(f, v)*dx(degree=8)
+    a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=12)
+    L = inner(f, v)*dx(degree=12)
 
     params = {"snes_type": "ksponly",
               "ksp_type": "preonly",
@@ -46,8 +46,6 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
               "snes_atol": 1e-25}
 
     # Compute solution
-    #assemble(a)
-    #assemble(L)
     sol = Function(V)
     solve(a == L, sol, solver_parameters=params)
 
@@ -78,7 +76,9 @@ def test_firedrake_helmholtz_serial():
                          [((1, (4, 6)), 1.9),
                           ((2, (3, 6)), 2.9),
                           ((3, (2, 4)), 3.9),
-                          ((4, (2, 4)), 4.7)])
+                          ((4, (2, 4)), 4.7),
+                          ((5, (2, 4)), 5.7),
+                          ((6, (2, 4)), 6.7)])
 def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_s(testcase, convrate):
     degree, (start, end) = testcase
     l2err = np.zeros(end - start)

From ffaab8746464c154b245e02dec3437230466d636 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sun, 28 Apr 2019 21:57:33 +0100
Subject: [PATCH 06/43] removed print statements

---
 tests/extrusion/test_helmholtz_serendipity_3d.py |  1 -
 tests/regression/test_helmholtz_serendipity.py   | 15 +++++++--------
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
index aba6f689fe..d49c78b042 100644
--- a/tests/extrusion/test_helmholtz_serendipity_3d.py
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -40,5 +40,4 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         solve(a == L, sol, solver_parameters=params)
 
         l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
-    print(np.array(np.log2(l2err[:-1]/l2err[1:])))
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
diff --git a/tests/regression/test_helmholtz_serendipity.py b/tests/regression/test_helmholtz_serendipity.py
index baf3165700..7f0ff82fb6 100644
--- a/tests/regression/test_helmholtz_serendipity.py
+++ b/tests/regression/test_helmholtz_serendipity.py
@@ -32,7 +32,7 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
     u = TrialFunction(V)
     v = TestFunction(V)
 
-    #uex = cos(pi*x)*cos(pi*y)
+    #uex = cos(pi*x)*cos(pi*y)  # Alternative problem
     uex = cos(x*pi*2)*cos(y*pi*2)
     f = -div(grad(uex)) + uex
 
@@ -67,9 +67,9 @@ def test_firedrake_helmholtz_serial():
     run_firedrake_helmholtz()
 
 
-#@pytest.mark.parallel
-#def test_firedrake_helmholtz_parallel():
-    #run_firedrake_helmholtz()
+@pytest.mark.parallel
+def test_firedrake_helmholtz_parallel():
+    run_firedrake_helmholtz()
 
 
 @pytest.mark.parametrize(('testcase', 'convrate'),
@@ -84,7 +84,6 @@ def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_s(testcase, co
     l2err = np.zeros(end - start)
     for ii in [i + start for i in range(len(l2err))]:
         l2err[ii - start] = helmholtz(ii, quadrilateral=True, degree=degree)[0]
-    print(l2err)
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
 
 
@@ -97,6 +96,6 @@ def test_firedrake_helmholtz_on_quadrilateral_mesh_from_file_serial():
     run_firedrake_helmholtz_on_quadrilateral_mesh_from_file()
 
 
-#@pytest.mark.parallel
-#def test_firedrake_helmholtz_on_quadrilateral_mesh_from_file_parallel():
-    #run_firedrake_helmholtz_on_quadrilateral_mesh_from_file()
+@pytest.mark.parallel
+def test_firedrake_helmholtz_on_quadrilateral_mesh_from_file_parallel():
+    run_firedrake_helmholtz_on_quadrilateral_mesh_from_file()

From bd48f1d0fdbb787a05fc5e06ee6a77b12725218d Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 29 Apr 2019 13:17:23 +0100
Subject: [PATCH 07/43] removed parallel comments

---
 tests/regression/test_helmholtz_serendipity.py | 4 +---
 tests/regression/test_steady_advection_2D.py   | 4 ++--
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/tests/regression/test_helmholtz_serendipity.py b/tests/regression/test_helmholtz_serendipity.py
index 7f0ff82fb6..25d53348e9 100644
--- a/tests/regression/test_helmholtz_serendipity.py
+++ b/tests/regression/test_helmholtz_serendipity.py
@@ -41,9 +41,7 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
 
     params = {"snes_type": "ksponly",
               "ksp_type": "preonly",
-              "pc_type": "lu",
-              "snes_rtol": 1e-16,
-              "snes_atol": 1e-25}
+              "pc_type": "lu"}
 
     # Compute solution
     sol = Function(V)
diff --git a/tests/regression/test_steady_advection_2D.py b/tests/regression/test_steady_advection_2D.py
index 1d574529a5..1169c5f3c9 100644
--- a/tests/regression/test_steady_advection_2D.py
+++ b/tests/regression/test_steady_advection_2D.py
@@ -69,7 +69,7 @@ def test_left_to_right(mesh, DG0, W):
     run_left_to_right(mesh, DG0, W)
 
 
-#@pytest.mark.parallel
+@pytest.mark.parallel
 def test_left_to_right_parallel(mesh, DG0, W):
     run_left_to_right(mesh, DG0, W)
 
@@ -106,6 +106,6 @@ def test_up_to_down(mesh, DG1, W):
     run_up_to_down(mesh, DG1, W)
 
 
-#@pytest.mark.parallel
+@pytest.mark.parallel
 def test_up_to_down_parallel(mesh, DG1, W):
     run_up_to_down(mesh, DG1, W)

From 8c9731a182b7bcdbb59f661182020127fe1d755f Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 29 Apr 2019 13:33:42 +0100
Subject: [PATCH 08/43] DROP BEFORE MERGE

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 5685529f05..4eb16f31e8 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -14,7 +14,7 @@ parentdir_prefix = firedrake-
 
 [flake8]
 ignore =
-    E501,F403,F405,E226,E402,E721,E731,E741,W503,F999,
+    E501,F403,F405,E226,E402,E721,E731,E741,W504,F999,
     N801,N802,N803,N806,N813,N814
 exclude = .git,__pycache__,build,.tox,dist,./pylit,firedrake/_version.py
 

From 9321f4b6ff45aa7cafc32f30e82ead542cf4d31d Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 29 Apr 2019 13:47:26 +0100
Subject: [PATCH 09/43] DROP BEFORE MERGE

---
 Jenkinsfile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Jenkinsfile b/Jenkinsfile
index 263ecdf503..c0d4dbb54c 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -26,7 +26,7 @@ pipeline {
         sh 'mkdir tmp'
         dir('tmp') {
           timestamps {
-            sh '../scripts/firedrake-install --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
+            sh '../scripts/firedrake-install --package-branch fiat Cyrus-serendipity-branch --package-branch ufl cyrus-serendipity --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
           }
         }
       }

From 99a76890ef63bed3da9887aeed63bd60b04ecdd4 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 20 May 2019 15:30:27 +0100
Subject: [PATCH 10/43] Test fixes

---
 tests/extrusion/test_helmholtz_serendipity_3d.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
index d49c78b042..92808da5bc 100644
--- a/tests/extrusion/test_helmholtz_serendipity_3d.py
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -10,7 +10,7 @@
                           (("S", 2, (3, 5)), 2.9),
                           (("S", 3, (2, 4)), 3.9),
                           (("S", 4, (2, 4)), 4.8),
-                          (("S", 5, (2, 5)), 5.7),
+                          (("S", 5, (2, 4)), 5.7),
                           (("S", 6, (2, 4)), 6.8),
                           (("S", 7, (2, 4)), 7.7)])
 def test_scalar_convergence(extmesh, testcase, convrate):
@@ -19,7 +19,7 @@ def test_scalar_convergence(extmesh, testcase, convrate):
     for ii in [i + start for i in range(len(l2err))]:
         mesh = extmesh(2**ii, 2**ii, 2**ii, quadrilateral=True)
 
-        fspace = FunctionSpace(mesh, family, degree, vfamily=family, vdegree=degree)
+        fspace = FunctionSpace(mesh, family, degree)
 
         u = TrialFunction(fspace)
         v = TestFunction(fspace)
@@ -29,8 +29,8 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         uex = cos(2*np.pi*x)*cos(2*np.pi*y)*cos(2*np.pi*z)
         f = -div(grad(uex)) + uex
 
-        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=16)
-        L = inner(f, v)*dx(degree=16)
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=3*degree)
+        L = inner(f, v)*dx(degree=3*degree)
 
         params = {"snes_type": "ksponly",
                   "ksp_type": "preonly",

From 12a0c8570776b0341ec4f1488fdfbdc87c13d596 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Mon, 20 May 2019 15:47:55 +0100
Subject: [PATCH 11/43] Correct Jenkinsfile

---
 Jenkinsfile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Jenkinsfile b/Jenkinsfile
index c0d4dbb54c..7875c41912 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -26,7 +26,7 @@ pipeline {
         sh 'mkdir tmp'
         dir('tmp') {
           timestamps {
-            sh '../scripts/firedrake-install --package-branch fiat Cyrus-serendipity-branch --package-branch ufl cyrus-serendipity --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
+            sh '../scripts/firedrake-install --package-branch fiat cyrus-serendipity --package-branch ufl cyrus-serendipity --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
           }
         }
       }

From 14513aec0d0a3633031438a4ca85902fac6cdb22 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Thu, 30 May 2019 16:12:02 +0100
Subject: [PATCH 12/43] new serendipity tests

---
 .../test_helmholtz_serendipity_2d.py          | 42 +++++++++++++++++++
 .../test_helmholtz_serendipity_3d.py          | 20 ++++-----
 2 files changed, 52 insertions(+), 10 deletions(-)
 create mode 100644 tests/extrusion/test_helmholtz_serendipity_2d.py

diff --git a/tests/extrusion/test_helmholtz_serendipity_2d.py b/tests/extrusion/test_helmholtz_serendipity_2d.py
new file mode 100644
index 0000000000..c773af361c
--- /dev/null
+++ b/tests/extrusion/test_helmholtz_serendipity_2d.py
@@ -0,0 +1,42 @@
+"""Tests for scalar Helmholtz convergence on Serendipity elements on extruded meshes"""
+import numpy as np
+import pytest
+
+from firedrake import *
+
+
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [(("S", 1, (4, 6)), 1.9),
+                          (("S", 2, (3, 5)), 2.9),
+                          (("S", 3, (2, 4)), 3.9),
+                          (("S", 4, (2, 4)), 4.8),
+                          (("S", 5, (2, 4)), 5.7)])
+def test_scalar_convergence(extmesh, testcase, convrate):
+    family, degree, (start, end) = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        mesh = UnitIntervalMesh(2**ii)
+        mesh = ExtrudedMesh(mesh, 2**ii)
+
+        fspace = FunctionSpace(mesh, family, degree)
+
+        u = TrialFunction(fspace)
+        v = TestFunction(fspace)
+
+        x, y = SpatialCoordinate(mesh)
+
+        uex = cos(2*np.pi*x)*cos(2*np.pi*y)
+        f = -div(grad(uex)) + uex
+
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=2*degree)
+        L = inner(f, v)*dx(degree=2*degree)
+
+        params = {"snes_type": "ksponly",
+                  "ksp_type": "preonly",
+                  "pc_type": "lu"}
+
+        sol = Function(fspace)
+        solve(a == L, sol, solver_parameters=params)
+
+        l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
+    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
index 92808da5bc..c404217bcc 100644
--- a/tests/extrusion/test_helmholtz_serendipity_3d.py
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -6,13 +6,12 @@
 
 
 @pytest.mark.parametrize(('testcase', 'convrate'),
-                         [(("S", 1, (4, 6)), 1.9),
-                          (("S", 2, (3, 5)), 2.9),
-                          (("S", 3, (2, 4)), 3.9),
-                          (("S", 4, (2, 4)), 4.8),
-                          (("S", 5, (2, 4)), 5.7),
-                          (("S", 6, (2, 4)), 6.8),
-                          (("S", 7, (2, 4)), 7.7)])
+                         [#(("S", 1, (4, 6)), 0.9),
+                          #(("S", 2, (3, 5)), 1.9),
+                          #(("S", 3, (2, 4)), 2.9),
+                          #(("S", 4, (2, 4)), 3.8),
+                          (("S", 5, (2, 4)), 4.7),
+                          (("S", 6, (2, 4)), 5.8)])
 def test_scalar_convergence(extmesh, testcase, convrate):
     family, degree, (start, end) = testcase
     l2err = np.zeros(end - start)
@@ -29,8 +28,8 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         uex = cos(2*np.pi*x)*cos(2*np.pi*y)*cos(2*np.pi*z)
         f = -div(grad(uex)) + uex
 
-        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=3*degree)
-        L = inner(f, v)*dx(degree=3*degree)
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=degree+7)
+        L = inner(f, v)*dx(degree=degree+7)
 
         params = {"snes_type": "ksponly",
                   "ksp_type": "preonly",
@@ -39,5 +38,6 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         sol = Function(fspace)
         solve(a == L, sol, solver_parameters=params)
 
-        l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
+        l2err[ii - start] = errornorm(sol, uex, norm_type="H1")
+        print(l2err)
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()

From 3749bbf9bd0fbe5c6ca95b81661f772dafeb2559 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Thu, 30 May 2019 16:12:30 +0100
Subject: [PATCH 13/43] new dpc tests

---
 tests/extrusion/test_steady_advection_2D_extr.py | 6 +++---
 tests/extrusion/test_steady_advection_3D_extr.py | 9 ++++++---
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/tests/extrusion/test_steady_advection_2D_extr.py b/tests/extrusion/test_steady_advection_2D_extr.py
index 7d2051071e..1eec1688bb 100644
--- a/tests/extrusion/test_steady_advection_2D_extr.py
+++ b/tests/extrusion/test_steady_advection_2D_extr.py
@@ -14,9 +14,9 @@ def mesh():
     return ExtrudedMesh(m, layers=4, layer_height=0.25)
 
 
-@pytest.fixture(scope='module')
-def DG0(mesh):
-    return FunctionSpace(mesh, "DG", 0)
+@pytest.fixture(scope='module', params=["DG", "DPC"])
+def DG0(request, mesh):
+    return FunctionSpace(mesh, request.param, 0)
 
 
 @pytest.fixture(scope='module')
diff --git a/tests/extrusion/test_steady_advection_3D_extr.py b/tests/extrusion/test_steady_advection_3D_extr.py
index 9fe500100a..bd57ab0b70 100644
--- a/tests/extrusion/test_steady_advection_3D_extr.py
+++ b/tests/extrusion/test_steady_advection_3D_extr.py
@@ -15,9 +15,12 @@ def mesh(request):
     return ExtrudedMesh(m, layers=4, layer_height=0.25)
 
 
-@pytest.fixture(scope='module')
-def DG0(mesh):
-    return FunctionSpace(mesh, "DG", 0)
+@pytest.fixture(scope='module', params=["DG", "DPC"])
+def DG0(request, mesh):
+    if mesh.ufl_cell() == triangle:
+        return FunctionSpace(mesh, "DG", 0)
+    else:
+        return FunctionSpace(mesh, request.param, 0)
 
 
 @pytest.fixture(scope='module')

From 4aaaf49691c028199638bd5153a5897957fda748 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Tue, 11 Jun 2019 15:04:26 +0100
Subject: [PATCH 14/43] test for serendipity 3d

---
 .../test_helmholtz_serendipity_3d.py          | 21 ++++++++++---------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
index c404217bcc..9d247c984b 100644
--- a/tests/extrusion/test_helmholtz_serendipity_3d.py
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -6,12 +6,12 @@
 
 
 @pytest.mark.parametrize(('testcase', 'convrate'),
-                         [#(("S", 1, (4, 6)), 0.9),
-                          #(("S", 2, (3, 5)), 1.9),
-                          #(("S", 3, (2, 4)), 2.9),
-                          #(("S", 4, (2, 4)), 3.8),
-                          (("S", 5, (2, 4)), 4.7),
-                          (("S", 6, (2, 4)), 5.8)])
+                         [(("S", 1, (4, 6)), 0.9),
+                          (("S", 2, (3, 5)), 1.8),
+                          (("S", 3, (2, 4)), 2.8),
+                          (("S", 4, (2, 4)), 3.7),
+                          (("S", 5, (3, 5)), 4.6),
+                          (("S", 6, (2, 4)), 5.6)])
 def test_scalar_convergence(extmesh, testcase, convrate):
     family, degree, (start, end) = testcase
     l2err = np.zeros(end - start)
@@ -25,11 +25,11 @@ def test_scalar_convergence(extmesh, testcase, convrate):
 
         x, y, z = SpatialCoordinate(mesh)
 
-        uex = cos(2*np.pi*x)*cos(2*np.pi*y)*cos(2*np.pi*z)
+        uex = cos(8*np.pi*x)*cos(8*np.pi*y)*cos(8*np.pi*z)
         f = -div(grad(uex)) + uex
 
-        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=degree+7)
-        L = inner(f, v)*dx(degree=degree+7)
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=degree+10)
+        L = inner(f, v)*dx(degree=degree+10)
 
         params = {"snes_type": "ksponly",
                   "ksp_type": "preonly",
@@ -38,6 +38,7 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         sol = Function(fspace)
         solve(a == L, sol, solver_parameters=params)
 
-        l2err[ii - start] = errornorm(sol, uex, norm_type="H1")
+        #l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
+        l2err[ii - start] = errornorm(uex, sol, norm_type="H1")
         print(l2err)
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()

From 2f064c775c48d027bc0a212a9969be9c78008cf9 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Tue, 11 Jun 2019 15:06:42 +0100
Subject: [PATCH 15/43] minor fix

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 4eb16f31e8..5685529f05 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -14,7 +14,7 @@ parentdir_prefix = firedrake-
 
 [flake8]
 ignore =
-    E501,F403,F405,E226,E402,E721,E731,E741,W504,F999,
+    E501,F403,F405,E226,E402,E721,E731,E741,W503,F999,
     N801,N802,N803,N806,N813,N814
 exclude = .git,__pycache__,build,.tox,dist,./pylit,firedrake/_version.py
 

From c71eafdf4b1efa61bbce11f9b5c4a3652edacc8d Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Tue, 11 Jun 2019 15:31:06 +0100
Subject: [PATCH 16/43] minor fix

---
 .../test_helmholtz_serendipity_3d.py          | 20 ++++++++-----------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/tests/extrusion/test_helmholtz_serendipity_3d.py b/tests/extrusion/test_helmholtz_serendipity_3d.py
index 9d247c984b..bdf5bd3970 100644
--- a/tests/extrusion/test_helmholtz_serendipity_3d.py
+++ b/tests/extrusion/test_helmholtz_serendipity_3d.py
@@ -6,12 +6,10 @@
 
 
 @pytest.mark.parametrize(('testcase', 'convrate'),
-                         [(("S", 1, (4, 6)), 0.9),
-                          (("S", 2, (3, 5)), 1.8),
-                          (("S", 3, (2, 4)), 2.8),
-                          (("S", 4, (2, 4)), 3.7),
-                          (("S", 5, (3, 5)), 4.6),
-                          (("S", 6, (2, 4)), 5.6)])
+                         [(("S", 1, (4, 6)), 1.9),
+                          (("S", 2, (3, 5)), 2.8),
+                          (("S", 3, (2, 4)), 3.8),
+                          (("S", 4, (2, 4)), 4.7)])
 def test_scalar_convergence(extmesh, testcase, convrate):
     family, degree, (start, end) = testcase
     l2err = np.zeros(end - start)
@@ -25,11 +23,11 @@ def test_scalar_convergence(extmesh, testcase, convrate):
 
         x, y, z = SpatialCoordinate(mesh)
 
-        uex = cos(8*np.pi*x)*cos(8*np.pi*y)*cos(8*np.pi*z)
+        uex = cos(2*np.pi*x)*cos(2*np.pi*y)*cos(2*np.pi*z)
         f = -div(grad(uex)) + uex
 
-        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=degree+10)
-        L = inner(f, v)*dx(degree=degree+10)
+        a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=degree+7)
+        L = inner(f, v)*dx(degree=degree+7)
 
         params = {"snes_type": "ksponly",
                   "ksp_type": "preonly",
@@ -38,7 +36,5 @@ def test_scalar_convergence(extmesh, testcase, convrate):
         sol = Function(fspace)
         solve(a == L, sol, solver_parameters=params)
 
-        #l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
-        l2err[ii - start] = errornorm(uex, sol, norm_type="H1")
-        print(l2err)
+        l2err[ii - start] = sqrt(assemble((sol-uex)*(sol-uex)*dx))
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()

From 1eca4db2fc88068c38af6e5a8ad42323d3e585eb Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Tue, 11 Jun 2019 15:46:49 +0100
Subject: [PATCH 17/43] removed serendipity test

---
 tests/regression/test_helmholtz.py | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/tests/regression/test_helmholtz.py b/tests/regression/test_helmholtz.py
index 25a49d4fbd..cfccda5c08 100644
--- a/tests/regression/test_helmholtz.py
+++ b/tests/regression/test_helmholtz.py
@@ -77,19 +77,6 @@ def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals(testcase, conv
     assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
 
 
-@pytest.mark.parametrize(('testcase', 'convrate'),
-                         [((1, (4, 6)), 1.9),
-                          ((2, (3, 6)), 2.9),
-                          ((3, (2, 4)), 3.9),
-                          ((4, (2, 4)), 4.7)])
-def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_dpc(testcase, convrate):
-    degree, (start, end) = testcase
-    l2err = np.zeros(end - start)
-    for ii in [i + start for i in range(len(l2err))]:
-        l2err[ii - start] = helmholtz(ii, quadrilateral=True, degree=degree, element="S")[0]
-    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
-
-
 def run_firedrake_helmholtz_on_quadrilateral_mesh_from_file():
     meshfile = join(cwd, "..", "meshes", "unitsquare_unstructured_quadrilaterals.msh")
     assert helmholtz(None, mesh=Mesh(meshfile))[0] <= 0.01

From f3bb21856436917aeee816a6147b1aba2f695803 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 21 Jun 2019 23:52:08 +0100
Subject: [PATCH 18/43] new tests

---
 Jenkinsfile~ | 123 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 123 insertions(+)
 create mode 100644 Jenkinsfile~

diff --git a/Jenkinsfile~ b/Jenkinsfile~
new file mode 100644
index 0000000000..c0d4dbb54c
--- /dev/null
+++ b/Jenkinsfile~
@@ -0,0 +1,123 @@
+pipeline {
+  agent {
+    docker {
+      image 'firedrakeproject/firedrake-env:latest'
+      label 'firedrakeproject'
+      args '-v /var/run/docker.sock:/var/run/docker.sock'
+      alwaysPull true
+    }
+  }
+  environment {
+    FIREDRAKE_CI_TESTS = "1"
+    DOCKER_CREDENTIALS = credentials('f52ccab9-5250-4b17-9fb6-c3f1ebdcc986')
+    PETSC_CONFIGURE_OPTIONS = "--with-make-np=11"
+  }
+  stages {
+    stage('Clean') {
+      steps {
+        sh 'git clean -fdx'
+        dir('tmp') {
+          deleteDir()
+        }
+      }
+    }
+    stage('Build') {
+      steps {
+        sh 'mkdir tmp'
+        dir('tmp') {
+          timestamps {
+            sh '../scripts/firedrake-install --package-branch fiat Cyrus-serendipity-branch --package-branch ufl cyrus-serendipity --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)'
+          }
+        }
+      }
+    }
+    stage('Test') {
+      parallel {
+        stage('Test Firedrake') {
+          steps {
+            dir('tmp') {
+              timestamps {
+                sh '''
+. ./firedrake/bin/activate
+python $(which firedrake-clean)
+python -m pip install pytest-cov pytest-xdist
+python -m pip list
+cd firedrake/src/firedrake
+python -m pytest -n 11 --cov firedrake -v tests
+'''
+              }
+            }
+          }
+        }
+        stage('Test pyadjoint'){
+          steps {
+            dir('tmp') {
+              timestamps {
+                sh '''
+. ./firedrake/bin/activate
+cd firedrake/src/pyadjoint; python -m pytest -v tests/firedrake_adjoint
+'''
+              }
+            }
+          }
+        }
+      }
+    }
+    stage('Post-test') {
+      parallel {
+        stage('Test build documentation') {
+          steps {
+            dir('tmp') {
+              timestamps {
+                sh '''
+. ./firedrake/bin/activate
+echo $PATH
+echo $VIRTUAL_ENV
+ls $VIRTUAL_ENV/bin
+cd firedrake/src/firedrake/docs; make html
+'''
+              }
+            }
+          }
+        }
+        stage('Zenodo API canary') {
+          steps {
+            timestamps {
+              sh 'scripts/firedrake-install --test-doi-resolution || (cat firedrake-install.log && /bin/false)'
+            }
+          }
+        }
+        stage('Lint') {
+          steps {
+            dir('tmp') {
+              timestamps {
+                sh '''
+. ./firedrake/bin/activate
+python -m pip install flake8
+cd firedrake/src/firedrake
+make lint
+'''
+              }
+            }
+          }
+        }
+      }
+    }
+    stage('Docker'){
+      when {
+        branch 'master'
+      }
+      steps {
+        sh '''
+sudo docker login -u $DOCKER_CREDENTIALS_USR -p $DOCKER_CREDENTIALS_PSW
+sudo docker build -t firedrakeproject/firedrake-env:latest -f docker/Dockerfile.env .
+sudo docker push firedrakeproject/firedrake-env:latest
+sudo docker build --no-cache --build-arg PETSC_CONFIGURE_OPTIONS -t firedrakeproject/firedrake-vanilla:latest -f docker/Dockerfile.vanilla .
+sudo docker push firedrakeproject/firedrake-vanilla:latest
+sudo docker build --no-cache --build-arg PETSC_CONFIGURE_OPTIONS -t firedrakeproject/firedrake:latest -f docker/Dockerfile.firedrake .
+sudo docker push firedrakeproject/firedrake:latest
+'''
+      }
+    }
+  }
+}

From 6343bb14b3ad0082d3d77e310531a754b041b222 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 21 Jun 2019 23:55:34 +0100
Subject: [PATCH 19/43] new tests

---
 tests/extrusion/test_serendipity_3d_polys.py  | 30 +++++++++++++++++++
 .../test_steady_advection_3D_extr.py          | 16 ++++++----
 tests/regression/test_helmholtz.py            |  4 +--
 .../regression/test_helmholtz_serendipity.py  |  4 +--
 4 files changed, 43 insertions(+), 11 deletions(-)
 create mode 100644 tests/extrusion/test_serendipity_3d_polys.py

diff --git a/tests/extrusion/test_serendipity_3d_polys.py b/tests/extrusion/test_serendipity_3d_polys.py
new file mode 100644
index 0000000000..2c33b81302
--- /dev/null
+++ b/tests/extrusion/test_serendipity_3d_polys.py
@@ -0,0 +1,30 @@
+import numpy as np
+import pytest
+
+from firedrake import *
+
+
+@pytest.mark.parametrize(('testcase'),
+                         [(5), (6)])
+def test_scalar_convergence(testcase):
+    mesh = UnitSquareMesh(2**2, 2**2, quadrilateral=True)
+    mesh = ExtrudedMesh(mesh, 2**2)
+
+    fspace = FunctionSpace(mesh, "S", testcase)
+
+    u = TrialFunction(fspace)
+    v = TestFunction(fspace)
+
+    x, y, z = SpatialCoordinate(mesh)
+
+    uex = x**testcase + y**testcase + x**2*y**3
+    f = uex
+
+    a = inner(u, v)*dx(degree=12)
+    L = inner(f, v)*dx(degree=12)
+
+    sol = Function(fspace)
+    solve(a == L, sol)
+
+    l2err = sqrt(assemble((sol-uex)*(sol-uex)*dx))
+    assert l2err < 1e-6
diff --git a/tests/extrusion/test_steady_advection_3D_extr.py b/tests/extrusion/test_steady_advection_3D_extr.py
index bd57ab0b70..cf8b321d0c 100644
--- a/tests/extrusion/test_steady_advection_3D_extr.py
+++ b/tests/extrusion/test_steady_advection_3D_extr.py
@@ -15,12 +15,16 @@ def mesh(request):
     return ExtrudedMesh(m, layers=4, layer_height=0.25)
 
 
-@pytest.fixture(scope='module', params=["DG", "DPC"])
-def DG0(request, mesh):
-    if mesh.ufl_cell() == triangle:
-        return FunctionSpace(mesh, "DG", 0)
-    else:
-        return FunctionSpace(mesh, request.param, 0)
+@pytest.fixture(scope='module')
+def DG0(mesh):
+    return FunctionSpace(mesh, "DG", 0)
+
+
+@pytest.fixture(scope='module')
+def DPC0():
+    m = UnitSquareMesh(4, 4, quadrilateral=True)
+    mesh = ExtrudedMesh(m, layers=4, layer_height=0.25)
+    return FunctionSpace(mesh, "DPC", 0)
 
 
 @pytest.fixture(scope='module')
diff --git a/tests/regression/test_helmholtz.py b/tests/regression/test_helmholtz.py
index cfccda5c08..75f1a692d9 100644
--- a/tests/regression/test_helmholtz.py
+++ b/tests/regression/test_helmholtz.py
@@ -20,11 +20,11 @@
 cwd = abspath(dirname(__file__))
 
 
-def helmholtz(r, quadrilateral=False, degree=2, mesh=None, element="CG"):
+def helmholtz(r, quadrilateral=False, degree=2, mesh=None):
     # Create mesh and define function space
     if mesh is None:
         mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral)
-    V = FunctionSpace(mesh, element, degree)
+    V = FunctionSpace(mesh, "CG", degree)
 
     # Define variational problem
     lmbda = 1
diff --git a/tests/regression/test_helmholtz_serendipity.py b/tests/regression/test_helmholtz_serendipity.py
index 25d53348e9..77578d3c86 100644
--- a/tests/regression/test_helmholtz_serendipity.py
+++ b/tests/regression/test_helmholtz_serendipity.py
@@ -32,7 +32,6 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
     u = TrialFunction(V)
     v = TestFunction(V)
 
-    #uex = cos(pi*x)*cos(pi*y)  # Alternative problem
     uex = cos(x*pi*2)*cos(y*pi*2)
     f = -div(grad(uex)) + uex
 
@@ -75,8 +74,7 @@ def test_firedrake_helmholtz_parallel():
                           ((2, (3, 6)), 2.9),
                           ((3, (2, 4)), 3.9),
                           ((4, (2, 4)), 4.7),
-                          ((5, (2, 4)), 5.7),
-                          ((6, (2, 4)), 6.7)])
+                          ((5, (2, 4)), 5.7)])
 def test_firedrake_helmholtz_scalar_convergence_on_quadrilaterals_s(testcase, convrate):
     degree, (start, end) = testcase
     l2err = np.zeros(end - start)

From 52001d241d062366bd4215fea7216d6381931e13 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sat, 22 Jun 2019 00:01:53 +0100
Subject: [PATCH 20/43] name change to DG0

---
 tests/extrusion/test_steady_advection_2D_extr.py | 2 +-
 tests/regression/test_steady_advection_2D.py     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/extrusion/test_steady_advection_2D_extr.py b/tests/extrusion/test_steady_advection_2D_extr.py
index 1eec1688bb..96b13113a2 100644
--- a/tests/extrusion/test_steady_advection_2D_extr.py
+++ b/tests/extrusion/test_steady_advection_2D_extr.py
@@ -15,7 +15,7 @@ def mesh():
 
 
 @pytest.fixture(scope='module', params=["DG", "DPC"])
-def DG0(request, mesh):
+def DGDPC0(request, mesh):
     return FunctionSpace(mesh, request.param, 0)
 
 
diff --git a/tests/regression/test_steady_advection_2D.py b/tests/regression/test_steady_advection_2D.py
index 1169c5f3c9..0d3ad198a7 100644
--- a/tests/regression/test_steady_advection_2D.py
+++ b/tests/regression/test_steady_advection_2D.py
@@ -15,7 +15,7 @@ def mesh(request):
 
 
 @pytest.fixture(params=["DG", "DPC"])
-def DG0(request, mesh):
+def DGDPC0(request, mesh):
     if mesh.ufl_cell() == triangle:
         return FunctionSpace(mesh, "DG", 0)
     else:

From aae3da1f5fdc06ab9c4f5184aa08c92a3579a773 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 12 Jul 2019 11:27:54 +0100
Subject: [PATCH 21/43] new tests for bdmc

---
 .../extrusion/test_point_eval_fs_extrusion.py | 20 +++++++++++++++++++
 tests/regression/test_poisson_mixed_no_bcs.py |  8 ++++----
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/tests/extrusion/test_point_eval_fs_extrusion.py b/tests/extrusion/test_point_eval_fs_extrusion.py
index 37a36a6909..9b6854eadf 100644
--- a/tests/extrusion/test_point_eval_fs_extrusion.py
+++ b/tests/extrusion/test_point_eval_fs_extrusion.py
@@ -63,6 +63,26 @@ def test_quad_vector(mesh_quad, family, degree):
     assert np.allclose([1.1, 0.18], f([0.0, 0.9]))
 
 
+@pytest.mark.parametrize(('family', 'degree'),
+                         [('BDMCE', 2)])
+def test_quad_vector(mesh_quad, family, degree):
+    x, y = SpatialCoordinate(mesh_quad)
+    V = FunctionSpace(mesh_quad, family, degree)
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    uex = as_vector([0.2 + y, 0.8*x + 0.2*y])
+
+    a = inner(u, v)*dx(degree=4)
+    L = inner(uex, v)*dx(degree=4)
+
+    sol = Function(V)
+    solve(a == L, sol)
+
+    l2err = errornorm(uex, sol, norm_type="l2")
+    assert l2err < 1e-6
+
+
 @pytest.mark.parametrize(('family', 'degree', 'vfamily', 'vdegree'),
                          [('CG', 3, 'DG', 2),
                           ('DG', 3, 'CG', 2)])
diff --git a/tests/regression/test_poisson_mixed_no_bcs.py b/tests/regression/test_poisson_mixed_no_bcs.py
index 8c369986b8..53804a9938 100644
--- a/tests/regression/test_poisson_mixed_no_bcs.py
+++ b/tests/regression/test_poisson_mixed_no_bcs.py
@@ -35,7 +35,7 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
 
     # Define function spaces and mixed (product) space
     if quadrilateral:
-        BDM = FunctionSpace(mesh, "RTCF", 1)
+        BDM = FunctionSpace(mesh, "BDMCF", 1)
     else:
         BDM = FunctionSpace(mesh, "BDM", 1)
     DG = FunctionSpace(mesh, "DG", 0)
@@ -49,8 +49,8 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
     # Define variational form
-    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
-    L = - f*v*dx
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=10)
+    L = - f*v*dx(degree=10)
 
     # Compute solution
     w = Function(W)
@@ -88,5 +88,5 @@ def test_hdiv_convergence(testcase, convrate):
     start, end = testcase
     l2err = np.zeros(end - start)
     for ii in [i + start for i in range(len(l2err))]:
-        l2err[ii - start] = poisson_mixed(ii, quadrilateral=True)[0]
+        l2err[ii - start]= poisson_mixed(ii, quadrilateral=True)[0]
     assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()

From 821feac685f64c496e23475611ed2a1f42d48d12 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 12 Jul 2019 14:40:49 +0100
Subject: [PATCH 22/43] new tests

---
 .../extrusion/test_point_eval_fs_extrusion.py  |  3 ++-
 tests/regression/test_poisson_mixed_no_bcs.py  |  8 ++++----
 tests/regression/test_steady_advection_2D.py   | 18 +++++++++---------
 3 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/tests/extrusion/test_point_eval_fs_extrusion.py b/tests/extrusion/test_point_eval_fs_extrusion.py
index 9b6854eadf..b4a735214c 100644
--- a/tests/extrusion/test_point_eval_fs_extrusion.py
+++ b/tests/extrusion/test_point_eval_fs_extrusion.py
@@ -64,7 +64,8 @@ def test_quad_vector(mesh_quad, family, degree):
 
 
 @pytest.mark.parametrize(('family', 'degree'),
-                         [('BDMCE', 2)])
+                         [('BDMCE', 2),
+                          ('BDMCF', 2)])
 def test_quad_vector(mesh_quad, family, degree):
     x, y = SpatialCoordinate(mesh_quad)
     V = FunctionSpace(mesh_quad, family, degree)
diff --git a/tests/regression/test_poisson_mixed_no_bcs.py b/tests/regression/test_poisson_mixed_no_bcs.py
index 53804a9938..8c369986b8 100644
--- a/tests/regression/test_poisson_mixed_no_bcs.py
+++ b/tests/regression/test_poisson_mixed_no_bcs.py
@@ -35,7 +35,7 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
 
     # Define function spaces and mixed (product) space
     if quadrilateral:
-        BDM = FunctionSpace(mesh, "BDMCF", 1)
+        BDM = FunctionSpace(mesh, "RTCF", 1)
     else:
         BDM = FunctionSpace(mesh, "BDM", 1)
     DG = FunctionSpace(mesh, "DG", 0)
@@ -49,8 +49,8 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
     # Define variational form
-    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=10)
-    L = - f*v*dx(degree=10)
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
+    L = - f*v*dx
 
     # Compute solution
     w = Function(W)
@@ -88,5 +88,5 @@ def test_hdiv_convergence(testcase, convrate):
     start, end = testcase
     l2err = np.zeros(end - start)
     for ii in [i + start for i in range(len(l2err))]:
-        l2err[ii - start]= poisson_mixed(ii, quadrilateral=True)[0]
+        l2err[ii - start] = poisson_mixed(ii, quadrilateral=True)[0]
     assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
diff --git a/tests/regression/test_steady_advection_2D.py b/tests/regression/test_steady_advection_2D.py
index 0d3ad198a7..1748901a73 100644
--- a/tests/regression/test_steady_advection_2D.py
+++ b/tests/regression/test_steady_advection_2D.py
@@ -34,20 +34,20 @@ def W(mesh):
         return FunctionSpace(mesh, "RTCF", 1)
 
 
-def run_left_to_right(mesh, DG0, W):
+def run_left_to_right(mesh, DGDPC0, W):
     velocity = as_vector((1.0, 0.0))
     u0 = project(velocity, W)
 
     xs = SpatialCoordinate(mesh)
     inflowexpr = conditional(And(xs[1] > 0.25, xs[1] < 0.75), 1.0, 0.5)
-    inflow = Function(DG0)
+    inflow = Function(DGDPC0)
     inflow.interpolate(inflowexpr)
 
     n = FacetNormal(mesh)
     un = 0.5*(dot(u0, n) + abs(dot(u0, n)))
 
-    D = TrialFunction(DG0)
-    phi = TestFunction(DG0)
+    D = TrialFunction(DGDPC0)
+    phi = TestFunction(DGDPC0)
 
     a1 = -D*dot(u0, grad(phi))*dx
     a2 = jump(phi)*(un('+')*D('+') - un('-')*D('-'))*dS
@@ -56,7 +56,7 @@ def run_left_to_right(mesh, DG0, W):
 
     L = -inflow*phi*dot(u0, n)*ds(1)  # inflow at left-hand wall
 
-    out = Function(DG0)
+    out = Function(DGDPC0)
     solve(a == L, out)
 
     # we only use inflow at the left wall, but since the velocity field
@@ -65,13 +65,13 @@ def run_left_to_right(mesh, DG0, W):
     assert max(abs(out.dat.data - inflow.dat.data)) < 1.2e-7
 
 
-def test_left_to_right(mesh, DG0, W):
-    run_left_to_right(mesh, DG0, W)
+def test_left_to_right(mesh, DGDPC0, W):
+    run_left_to_right(mesh, DGDPC0, W)
 
 
 @pytest.mark.parallel
-def test_left_to_right_parallel(mesh, DG0, W):
-    run_left_to_right(mesh, DG0, W)
+def test_left_to_right_parallel(mesh, DGDPC0, W):
+    run_left_to_right(mesh, DGDPC0, W)
 
 
 def run_up_to_down(mesh, DG1, W):

From ea136edb5d25f185fc077a2654ef01b5cadc4564 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sat, 27 Jul 2019 20:55:03 +0100
Subject: [PATCH 23/43] new bdmc tests

---
 .../extrusion/test_point_eval_fs_extrusion.py | 15 +--
 .../test_poisson_mixed_no_bcs_bdmc.py         | 93 ++++++++++++++++++
 .../test_poisson_mixed_strong_bcs_bdmc.py     | 97 +++++++++++++++++++
 3 files changed, 194 insertions(+), 11 deletions(-)
 create mode 100644 tests/regression/test_poisson_mixed_no_bcs_bdmc.py
 create mode 100644 tests/regression/test_poisson_mixed_strong_bcs_bdmc.py

diff --git a/tests/extrusion/test_point_eval_fs_extrusion.py b/tests/extrusion/test_point_eval_fs_extrusion.py
index b4a735214c..501eb61910 100644
--- a/tests/extrusion/test_point_eval_fs_extrusion.py
+++ b/tests/extrusion/test_point_eval_fs_extrusion.py
@@ -69,19 +69,12 @@ def test_quad_vector(mesh_quad, family, degree):
 def test_quad_vector(mesh_quad, family, degree):
     x, y = SpatialCoordinate(mesh_quad)
     V = FunctionSpace(mesh_quad, family, degree)
-    u = TrialFunction(V)
-    v = TestFunction(V)
+    f = Function(V).project(as_vector([0.2 + y, 0.8*x + 0.2*y]))
 
-    uex = as_vector([0.2 + y, 0.8*x + 0.2*y])
+    W = VectorFunctionSpace(mesh_quad, "DG", 1)
+    g = Function(W).interpolate(as_vector([0.2 + y, 0.8*x + 0.2*y]))
 
-    a = inner(u, v)*dx(degree=4)
-    L = inner(uex, v)*dx(degree=4)
-
-    sol = Function(V)
-    solve(a == L, sol)
-
-    l2err = errornorm(uex, sol, norm_type="l2")
-    assert l2err < 1e-6
+    assert sqrt(assemble(dot(g - f, g - f) * dx)) <1e-6
 
 
 @pytest.mark.parametrize(('family', 'degree', 'vfamily', 'vdegree'),
diff --git a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
new file mode 100644
index 0000000000..d0fd6b77c9
--- /dev/null
+++ b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
@@ -0,0 +1,93 @@
+"""Solve the mixed formulation of the Laplacian described in section 2.3.1 of
+Arnold, Falk, Winther 2010 "Finite Element Exterior Calculus: From Hodge Theory
+to Numerical Stability":
+
+    sigma - grad(u) = 0
+         div(sigma) = f
+
+The corresponding weak (variational problem)
+
+    <sigma, tau> + <div(tau), u>   = 0       for all tau
+                   <div(sigma), v> = <f, v>  for all v
+
+is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
+(sigma, tau) and DG (discontinuous Galerkin) elements of degree k - 1
+for (u, v).
+
+No strong boundary conditions are enforced. The forcing function is chosen as
+
+    -2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1]
+
+which reproduces the known analytical solution
+
+    x[0]*(1-x[0])*x[1]*(1-x[1])
+"""
+import pytest
+import numpy as np
+
+from firedrake import *
+
+
+def poisson_mixed(size, parameters={}, quadrilateral=False):
+    # Create mesh
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=quadrilateral)
+    x = SpatialCoordinate(mesh)
+
+    # Define function spaces and mixed (product) space
+    if quadrilateral:
+        BDM = FunctionSpace(mesh, "BDMCF", 1)
+    else:
+        BDM = FunctionSpace(mesh, "BDM", 1)
+    DG = FunctionSpace(mesh, "DG", 0)
+    W = BDM * DG
+
+    # Define trial and test functions
+    sigma, u = TrialFunctions(W)
+    tau, v = TestFunctions(W)
+
+    # Define source function
+    f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
+
+    # Define variational form
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
+    L = -f*v*dx(degree=4)
+
+    # Compute solution
+    w = Function(W)
+    solve(a == L, w, solver_parameters=parameters)
+    sigma, u = w.split()
+
+    # Analytical solution
+    f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1]))
+    return sqrt(assemble(dot(u - f, u - f) * dx)), u, f
+
+
+@pytest.mark.parametrize('parameters',
+                         [{}, {'pc_type': 'fieldsplit',
+                               'pc_fieldsplit_type': 'schur',
+                               'ksp_type': 'preonly',
+                               'pc_fieldsplit_schur_fact_type': 'FULL',
+                               'fieldsplit_0_ksp_type': 'cg',
+                               'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
+                               'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
+                               'fieldsplit_1_ksp_type': 'cg'}])
+def test_poisson_mixed(parameters):
+    """Test second-order convergence of the mixed poisson formulation."""
+    diff = np.array([poisson_mixed(i, parameters)[0] for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > 1.9).all()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [((3, 6), 1.9)])
+def test_hdiv_convergence(testcase, convrate):
+    """Test second-order convergence of the mixed poisson formulation
+    on quadrilaterals with H(div) elements."""
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = poisson_mixed(ii, quadrilateral=True)[0]
+    print(l2err)
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
diff --git a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
new file mode 100644
index 0000000000..72dab2a484
--- /dev/null
+++ b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
@@ -0,0 +1,97 @@
+"""Solve the mixed formulation of the Laplacian on the unit square
+
+    sigma - grad(u) = 0
+         div(sigma) = f
+
+The corresponding weak (variational problem)
+
+    <sigma, tau> + <div(tau), u>   = 0       for all tau
+                   <div(sigma), v> = <f, v>  for all v
+
+is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
+(sigma, tau) and DG (discontinuous Galerkin) elements of degree k - 1
+for (u, v).
+
+The boundary conditions on the left and right are enforced strongly as
+
+    dot(sigma, n) = 0
+
+which corresponds to a Neumann condition du/dn = 0.
+
+The top is fixed to 42 with a Dirichlet boundary condition, which enters
+the weak formulation of the right hand side as
+
+    42*dot(tau, n)*ds
+"""
+import pytest
+from firedrake import *
+
+
+def poisson_mixed(size, parameters={}, quadrilateral=False):
+    # Create mesh
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=quadrilateral)
+    x = SpatialCoordinate(mesh)
+
+    # Define function spaces and mixed (product) space
+    BDM = FunctionSpace(mesh, "BDM" if not quadrilateral else "BDMCF", 1)
+    DG = FunctionSpace(mesh, "DG", 0)
+    W = BDM * DG
+
+    # Define trial and test functions
+    sigma, u = TrialFunctions(W)
+    tau, v = TestFunctions(W)
+
+    # Define source function
+    f = Function(DG).assign(0)
+
+    # Define variational form
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
+    n = FacetNormal(mesh)
+    L = -f*v*dx + 42*dot(tau, n)*ds(4)
+
+    # Apply dot(sigma, n) == 0 on left and right boundaries strongly
+    # (corresponding to Neumann condition du/dn = 0)
+    bcs = DirichletBC(W.sub(0), Constant((0, 0)), (1, 2))
+    # Compute solution
+    w = Function(W)
+    solve(a == L, w, bcs=bcs, solver_parameters=parameters)
+    sigma, u = w.split()
+
+    # Analytical solution
+    f.interpolate(42*x[1])
+    return sqrt(assemble(dot(u - f, u - f) * dx))
+
+
+@pytest.mark.parametrize('quadrilateral', [False, True])
+@pytest.mark.parametrize('parameters',
+                         [{},
+                          {'ksp_type': 'fgmres',
+                           'pc_type': 'fieldsplit',
+                           'pc_fieldsplit_type': 'schur',
+                           'pc_fieldsplit_schur_fact_type': 'diag',
+                           'fieldsplit_0_ksp_type': 'preonly',
+                           'fieldsplit_0_pc_type': 'icc',
+                           'fieldsplit_1_ksp_type': 'cg',
+                           'fieldsplit_1_pc_type': 'none'}])
+def test_poisson_mixed(parameters, quadrilateral):
+    assert poisson_mixed(3, parameters, quadrilateral=quadrilateral) < 2e-5
+
+
+@pytest.mark.parallel(nprocs=3)
+def test_poisson_mixed_parallel_fieldsplit():
+    x = poisson_mixed(3, parameters={'ksp_type': 'fgmres',
+                                     'pc_type': 'fieldsplit',
+                                     'pc_fieldsplit_type': 'schur',
+                                     'fieldsplit_schur_fact_type': 'diag',
+                                     'fieldsplit_0_ksp_type': 'preonly',
+                                     'fieldsplit_1_ksp_type': 'cg',
+                                     'fieldsplit_0_pc_type': 'bjacobi',
+                                     'fieldsplit_0_sub_pc_type': 'icc',
+                                     'fieldsplit_1_pc_type': 'none'})
+    assert x < 2e-5
+
+
+@pytest.mark.parallel(nprocs=3)
+@pytest.mark.parametrize('quadrilateral', [False, True])
+def test_poisson_mixed_parallel(quadrilateral):
+    assert poisson_mixed(3, quadrilateral=quadrilateral) < 2e-5

From 73541c57693804de82b08cbaa097fc9a703ccf88 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sun, 28 Jul 2019 03:25:04 +0100
Subject: [PATCH 24/43] bdmcf test

---
 tests/regression/test_bdmcf.py | 47 ++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 tests/regression/test_bdmcf.py

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
new file mode 100644
index 0000000000..ceec3756fd
--- /dev/null
+++ b/tests/regression/test_bdmcf.py
@@ -0,0 +1,47 @@
+from firedrake import *
+import matplotlib.pyplot as plt
+
+mesh = UnitSquareMesh(1, 1, quadrilateral=True)
+x = SpatialCoordinate(mesh)
+
+BDM = FunctionSpace(mesh, "BDMCF", 3)
+DG = FunctionSpace(mesh, "DG", 2)
+W = BDM * DG
+
+# Define trial and test functions
+sigma, u = TrialFunctions(W)
+tau, v = TestFunctions(W)
+
+# Define source function
+f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
+
+# Define variational form
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
+L = -f*v*dx(degree=4)
+
+# Compute solution
+w = Function(W)
+solve(a == L, w, solver_parameters={'pc_type': 'fieldsplit',
+      'pc_fieldsplit_type': 'schur',
+      'ksp_type': 'preonly',
+      'pc_fieldsplit_schur_fact_type': 'FULL',
+      'fieldsplit_0_ksp_type': 'cg',
+      'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
+      'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
+      'fieldsplit_1_ksp_type': 'cg'})
+# solve(a == L, w)
+sigma, u = w.split()
+
+# Analytical solution
+f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1]))
+
+l2err = errornorm(f, u, norm_type="l2")
+print(l2err)
+
+plt.figure(1)
+plot(u)
+
+plt.figure(2)
+plot(f)
+
+plt.show()

From 885b256cbfc1088a44764e2ff0cc13d9b51344fc Mon Sep 17 00:00:00 2001
From: David Ham <david.ham@imperial.ac.uk>
Date: Thu, 8 Aug 2019 14:55:12 +0100
Subject: [PATCH 25/43] Correct functionspace and quadrature degree

---
 tests/regression/test_bdmcf.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index ceec3756fd..646b710db9 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -5,7 +5,7 @@
 x = SpatialCoordinate(mesh)
 
 BDM = FunctionSpace(mesh, "BDMCF", 3)
-DG = FunctionSpace(mesh, "DG", 2)
+DG = FunctionSpace(mesh, "DP", 2)
 W = BDM * DG
 
 # Define trial and test functions
@@ -16,8 +16,8 @@
 f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
 # Define variational form
-a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
-L = -f*v*dx(degree=4)
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=6)
+L = -f*v*dx(degree=6)
 
 # Compute solution
 w = Function(W)

From d96768c2b7fdc896784a640cc2fc3fa52ae4cfe1 Mon Sep 17 00:00:00 2001
From: David Ham <david.ham@imperial.ac.uk>
Date: Thu, 8 Aug 2019 17:00:06 +0100
Subject: [PATCH 26/43] passing bdmce/f test

---
 tests/regression/test_bdmc.py | 46 +++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 tests/regression/test_bdmc.py

diff --git a/tests/regression/test_bdmc.py b/tests/regression/test_bdmc.py
new file mode 100644
index 0000000000..2a3cbb56e2
--- /dev/null
+++ b/tests/regression/test_bdmc.py
@@ -0,0 +1,46 @@
+import numpy as np
+from firedrake import *
+import pytest
+
+
+def project_bdmc(size, degree, family):
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
+    x = SpatialCoordinate(mesh)
+
+    fs = FunctionSpace(mesh, family, degree)
+
+    f = Function(fs)
+    v = TestFunction(fs)
+
+    expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
+                      cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
+
+    solve(inner(f-expr,v) * dx(degree=8) == 0, f)
+
+    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8))) 
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmcf(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCF")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmce(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCE")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()

From 1ec3acc54142f73f14eeddfb2bd7af205cd89894 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Thu, 8 Aug 2019 21:51:25 +0100
Subject: [PATCH 27/43] minor updates

---
 tests/regression/test_bdmcf.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index ceec3756fd..f49a896da2 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -1,7 +1,7 @@
 from firedrake import *
 import matplotlib.pyplot as plt
 
-mesh = UnitSquareMesh(1, 1, quadrilateral=True)
+mesh = UnitSquareMesh(32, 32, quadrilateral=True)
 x = SpatialCoordinate(mesh)
 
 BDM = FunctionSpace(mesh, "BDMCF", 3)
@@ -16,8 +16,8 @@
 f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
 # Define variational form
-a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
-L = -f*v*dx(degree=4)
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=12)
+L = -f*v*dx(degree=12)
 
 # Compute solution
 w = Function(W)

From e27bded03b6195cb48049b3dd6a6d3f023bf8bf7 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 03:03:29 +0100
Subject: [PATCH 28/43] minor changes

---
 tests/regression/test_bdmcf.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index f1f398a4f7..d2f5889064 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -1,11 +1,11 @@
 from firedrake import *
 import matplotlib.pyplot as plt
 
-mesh = UnitSquareMesh(32, 32, quadrilateral=True)
+mesh = UnitSquareMesh(16, 16, quadrilateral=True)
 x = SpatialCoordinate(mesh)
 
-BDM = FunctionSpace(mesh, "BDMCF", 3)
-DG = FunctionSpace(mesh, "DP", 2)
+BDM = FunctionSpace(mesh, "BDMCF", 1)
+DG = FunctionSpace(mesh, "DPC", 0)
 W = BDM * DG
 
 # Define trial and test functions
@@ -16,8 +16,8 @@
 f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
 # Define variational form
-a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=6)
-L = -f*v*dx(degree=6)
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=10)
+L = -f*v*dx(degree=10)
 
 # Compute solution
 w = Function(W)

From 204891d4655fecd0913f89293dce4408eb22b1f8 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 10:38:16 +0100
Subject: [PATCH 29/43] correct functionspace name to DGDPC0

---
 .../test_steady_advection_2D_extr.py          | 36 +++++++++----------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/tests/extrusion/test_steady_advection_2D_extr.py b/tests/extrusion/test_steady_advection_2D_extr.py
index 96b13113a2..63f1cac198 100644
--- a/tests/extrusion/test_steady_advection_2D_extr.py
+++ b/tests/extrusion/test_steady_advection_2D_extr.py
@@ -38,21 +38,21 @@ def W(mesh):
     return FunctionSpace(mesh, W0+W1)
 
 
-def run_left_to_right(mesh, DG0, W):
+def run_left_to_right(mesh, DGDPC0, W):
     velocity = as_vector([1.0, 0.0])
     u0 = project(velocity, W)
 
     xs = SpatialCoordinate(mesh)
 
     inflowexpr = conditional(And(xs[1] > 0.25, xs[1] < 0.75), 1.0, 0.5)
-    inflow = Function(DG0)
+    inflow = Function(DGDPC0)
     inflow.interpolate(inflowexpr)
 
     n = FacetNormal(mesh)
     un = 0.5*(dot(u0, n) + abs(dot(u0, n)))
 
-    D = TrialFunction(DG0)
-    phi = TestFunction(DG0)
+    D = TrialFunction(DGDPC0)
+    phi = TestFunction(DGDPC0)
 
     a1 = -D*dot(u0, grad(phi))*dx
     a2 = jump(phi)*(un('+')*D('+') - un('-')*D('-'))*dS_v
@@ -61,7 +61,7 @@ def run_left_to_right(mesh, DG0, W):
 
     L = -inflow*phi*dot(u0, n)*ds_v(1)  # inflow at left-hand wall
 
-    out = Function(DG0)
+    out = Function(DGDPC0)
     solve(a == L, out)
 
     # we only use inflow at the left wall, but since the velocity field
@@ -70,13 +70,13 @@ def run_left_to_right(mesh, DG0, W):
     assert max(abs(out.dat.data - inflow.dat.data)) < 1e-14
 
 
-def test_left_to_right(mesh, DG0, W):
-    run_left_to_right(mesh, DG0, W)
+def test_left_to_right(mesh, DGDPC0, W):
+    run_left_to_right(mesh, DGDPC0, W)
 
 
 @pytest.mark.parallel
-def test_left_to_right_parallel(mesh, DG0, W):
-    run_left_to_right(mesh, DG0, W)
+def test_left_to_right_parallel(mesh, DGDPC0, W):
+    run_left_to_right(mesh, DGDPC0, W)
 
 
 def run_right_to_left(mesh, DG1, W):
@@ -116,20 +116,20 @@ def test_right_to_left_parallel(mesh, DG1, W):
     run_right_to_left(mesh, DG1, W)
 
 
-def run_bottom_to_top(mesh, DG0, W):
+def run_bottom_to_top(mesh, DGDPC0, W):
     velocity = as_vector([0.0, 1.0])
     u0 = project(velocity, W)
 
     xs = SpatialCoordinate(mesh)
     inflowexpr = conditional(And(xs[0] > 0.25, xs[0] < 0.75), 1.0, 0.5)
-    inflow = Function(DG0)
+    inflow = Function(DGDPC0)
     inflow.interpolate(inflowexpr)
 
     n = FacetNormal(mesh)
     un = 0.5*(dot(u0, n) + abs(dot(u0, n)))
 
-    D = TrialFunction(DG0)
-    phi = TestFunction(DG0)
+    D = TrialFunction(DGDPC0)
+    phi = TestFunction(DGDPC0)
 
     a1 = -D*dot(u0, grad(phi))*dx
     a2 = jump(phi)*(un('+')*D('+') - un('-')*D('-'))*dS_h
@@ -138,19 +138,19 @@ def run_bottom_to_top(mesh, DG0, W):
 
     L = -inflow*phi*dot(u0, n)*ds_b  # inflow at bottom wall
 
-    out = Function(DG0)
+    out = Function(DGDPC0)
     solve(a == L, out)
 
     assert max(abs(out.dat.data - inflow.dat.data)) < 1e-14
 
 
-def test_bottom_to_top(mesh, DG0, W):
-    run_bottom_to_top(mesh, DG0, W)
+def test_bottom_to_top(mesh, DGDPC0, W):
+    run_bottom_to_top(mesh, DGDPC0, W)
 
 
 @pytest.mark.parallel
-def test_bottom_to_top_parallel(mesh, DG0, W):
-    run_bottom_to_top(mesh, DG0, W)
+def test_bottom_to_top_parallel(mesh, DGDPC0, W):
+    run_bottom_to_top(mesh, DGDPC0, W)
 
 
 def run_top_to_bottom(mesh, DG1, W):

From 13844df800f1a297925773723b636ce1519eb319 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 10:54:17 +0100
Subject: [PATCH 30/43] rename DG to DPC

---
 tests/regression/test_bdmcf.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index d2f5889064..1974ba6a8c 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -5,15 +5,15 @@
 x = SpatialCoordinate(mesh)
 
 BDM = FunctionSpace(mesh, "BDMCF", 1)
-DG = FunctionSpace(mesh, "DPC", 0)
-W = BDM * DG
+DPC = FunctionSpace(mesh, "DPC", 0)
+W = BDM * DPC
 
 # Define trial and test functions
 sigma, u = TrialFunctions(W)
 tau, v = TestFunctions(W)
 
 # Define source function
-f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
+f = Function(DPC).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
 # Define variational form
 a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=10)

From eb5d8e80c42c3b19feb8b6a737a6979ce47ddcd4 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 14:45:11 +0100
Subject: [PATCH 31/43] remove bdmc tests

---
 tests/extrusion/test_point_eval_fs_extrusion.py | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/tests/extrusion/test_point_eval_fs_extrusion.py b/tests/extrusion/test_point_eval_fs_extrusion.py
index 501eb61910..37a36a6909 100644
--- a/tests/extrusion/test_point_eval_fs_extrusion.py
+++ b/tests/extrusion/test_point_eval_fs_extrusion.py
@@ -63,20 +63,6 @@ def test_quad_vector(mesh_quad, family, degree):
     assert np.allclose([1.1, 0.18], f([0.0, 0.9]))
 
 
-@pytest.mark.parametrize(('family', 'degree'),
-                         [('BDMCE', 2),
-                          ('BDMCF', 2)])
-def test_quad_vector(mesh_quad, family, degree):
-    x, y = SpatialCoordinate(mesh_quad)
-    V = FunctionSpace(mesh_quad, family, degree)
-    f = Function(V).project(as_vector([0.2 + y, 0.8*x + 0.2*y]))
-
-    W = VectorFunctionSpace(mesh_quad, "DG", 1)
-    g = Function(W).interpolate(as_vector([0.2 + y, 0.8*x + 0.2*y]))
-
-    assert sqrt(assemble(dot(g - f, g - f) * dx)) <1e-6
-
-
 @pytest.mark.parametrize(('family', 'degree', 'vfamily', 'vdegree'),
                          [('CG', 3, 'DG', 2),
                           ('DG', 3, 'CG', 2)])

From ac1f55e666a909de54b1d22a3e7c166587588b94 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 14:45:35 +0100
Subject: [PATCH 32/43] test BDMCF1 DPC0 case

---
 .../test_poisson_mixed_no_bcs_bdmc.py         | 38 +++++++++----------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
index d0fd6b77c9..abf889bdd1 100644
--- a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
+++ b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
@@ -38,7 +38,7 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
         BDM = FunctionSpace(mesh, "BDMCF", 1)
     else:
         BDM = FunctionSpace(mesh, "BDM", 1)
-    DG = FunctionSpace(mesh, "DG", 0)
+    DG = FunctionSpace(mesh, "DPC", 0)
     W = BDM * DG
 
     # Define trial and test functions
@@ -49,8 +49,8 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
     # Define variational form
-    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
-    L = -f*v*dx(degree=4)
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=6)
+    L = -f*v*dx(degree=6)
 
     # Compute solution
     w = Function(W)
@@ -62,22 +62,22 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     return sqrt(assemble(dot(u - f, u - f) * dx)), u, f
 
 
-@pytest.mark.parametrize('parameters',
-                         [{}, {'pc_type': 'fieldsplit',
-                               'pc_fieldsplit_type': 'schur',
-                               'ksp_type': 'preonly',
-                               'pc_fieldsplit_schur_fact_type': 'FULL',
-                               'fieldsplit_0_ksp_type': 'cg',
-                               'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
-                               'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
-                               'fieldsplit_1_ksp_type': 'cg'}])
-def test_poisson_mixed(parameters):
-    """Test second-order convergence of the mixed poisson formulation."""
-    diff = np.array([poisson_mixed(i, parameters)[0] for i in range(3, 6)])
-    print("l2 error norms:", diff)
-    conv = np.log2(diff[:-1] / diff[1:])
-    print("convergence order:", conv)
-    assert (np.array(conv) > 1.9).all()
+# @pytest.mark.parametrize('parameters',
+#                          [{}, {'pc_type': 'fieldsplit',
+#                                'pc_fieldsplit_type': 'schur',
+#                                'ksp_type': 'preonly',
+#                                'pc_fieldsplit_schur_fact_type': 'FULL',
+#                                'fieldsplit_0_ksp_type': 'cg',
+#                                'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
+#                                'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
+#                                'fieldsplit_1_ksp_type': 'cg'}])
+# def test_poisson_mixed(parameters):
+#     """Test second-order convergence of the mixed poisson formulation."""
+#     diff = np.array([poisson_mixed(i, parameters)[0] for i in range(3, 6)])
+#     print("l2 error norms:", diff)
+#     conv = np.log2(diff[:-1] / diff[1:])
+#     print("convergence order:", conv)
+#     assert (np.array(conv) > 1.9).all()
 
 
 @pytest.mark.parametrize(('testcase', 'convrate'),

From 86b2c8016ed587b4f5497ced3176b7b95382fb33 Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 14:45:54 +0100
Subject: [PATCH 33/43] change DG to DPC

---
 tests/regression/test_poisson_mixed_strong_bcs_bdmc.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
index 72dab2a484..aaf4e27881 100644
--- a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
+++ b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
@@ -34,15 +34,15 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
 
     # Define function spaces and mixed (product) space
     BDM = FunctionSpace(mesh, "BDM" if not quadrilateral else "BDMCF", 1)
-    DG = FunctionSpace(mesh, "DG", 0)
-    W = BDM * DG
+    DPC = FunctionSpace(mesh, "DPC", 0)
+    W = BDM * DPC
 
     # Define trial and test functions
     sigma, u = TrialFunctions(W)
     tau, v = TestFunctions(W)
 
     # Define source function
-    f = Function(DG).assign(0)
+    f = Function(DPC).assign(0)
 
     # Define variational form
     a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)

From 06e9b34db80557744743718ba0aee25f084d202f Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Fri, 9 Aug 2019 14:50:11 +0100
Subject: [PATCH 34/43] change DG to DPC

---
 .../test_poisson_mixed_strong_bcs_bdmc.py        | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
index aaf4e27881..c4bfee31df 100644
--- a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
+++ b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
@@ -27,13 +27,13 @@
 from firedrake import *
 
 
-def poisson_mixed(size, parameters={}, quadrilateral=False):
+def poisson_mixed(size, parameters={}):
     # Create mesh
-    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=quadrilateral)
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
     x = SpatialCoordinate(mesh)
 
     # Define function spaces and mixed (product) space
-    BDM = FunctionSpace(mesh, "BDM" if not quadrilateral else "BDMCF", 1)
+    BDM = FunctionSpace(mesh, "BDMCF", 1)
     DPC = FunctionSpace(mesh, "DPC", 0)
     W = BDM * DPC
 
@@ -62,7 +62,6 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     return sqrt(assemble(dot(u - f, u - f) * dx))
 
 
-@pytest.mark.parametrize('quadrilateral', [False, True])
 @pytest.mark.parametrize('parameters',
                          [{},
                           {'ksp_type': 'fgmres',
@@ -73,8 +72,8 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
                            'fieldsplit_0_pc_type': 'icc',
                            'fieldsplit_1_ksp_type': 'cg',
                            'fieldsplit_1_pc_type': 'none'}])
-def test_poisson_mixed(parameters, quadrilateral):
-    assert poisson_mixed(3, parameters, quadrilateral=quadrilateral) < 2e-5
+def test_poisson_mixed(parameters):
+    assert poisson_mixed(3, parameters) < 2e-5
 
 
 @pytest.mark.parallel(nprocs=3)
@@ -92,6 +91,5 @@ def test_poisson_mixed_parallel_fieldsplit():
 
 
 @pytest.mark.parallel(nprocs=3)
-@pytest.mark.parametrize('quadrilateral', [False, True])
-def test_poisson_mixed_parallel(quadrilateral):
-    assert poisson_mixed(3, quadrilateral=quadrilateral) < 2e-5
+def test_poisson_mixed_parallel():
+    assert poisson_mixed(3) < 2e-5

From a7ac0c017efcfb35742d092a544c38a8c1c2e75f Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Thu, 15 Aug 2019 13:09:00 +0100
Subject: [PATCH 35/43] test degree 2 bdmc

---
 tests/regression/test_bdmcf.py                |  4 +-
 .../test_poisson_mixed_no_bcs_bdmc.py         | 38 ++++---------------
 2 files changed, 10 insertions(+), 32 deletions(-)

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index 1974ba6a8c..1e0904c502 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -4,8 +4,8 @@
 mesh = UnitSquareMesh(16, 16, quadrilateral=True)
 x = SpatialCoordinate(mesh)
 
-BDM = FunctionSpace(mesh, "BDMCF", 1)
-DPC = FunctionSpace(mesh, "DPC", 0)
+BDM = FunctionSpace(mesh, "BDMCF", 2)
+DPC = FunctionSpace(mesh, "DPC", 1)
 W = BDM * DPC
 
 # Define trial and test functions
diff --git a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
index abf889bdd1..d3172a8e06 100644
--- a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
+++ b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
@@ -11,7 +11,7 @@
                    <div(sigma), v> = <f, v>  for all v
 
 is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
-(sigma, tau) and DG (discontinuous Galerkin) elements of degree k - 1
+(sigma, tau) and DPC (discontinuous Galerkin) elements of degree k - 1
 for (u, v).
 
 No strong boundary conditions are enforced. The forcing function is chosen as
@@ -28,25 +28,22 @@
 from firedrake import *
 
 
-def poisson_mixed(size, parameters={}, quadrilateral=False):
+def poisson_mixed(size, parameters={}):
     # Create mesh
-    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=quadrilateral)
+    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
     x = SpatialCoordinate(mesh)
 
     # Define function spaces and mixed (product) space
-    if quadrilateral:
-        BDM = FunctionSpace(mesh, "BDMCF", 1)
-    else:
-        BDM = FunctionSpace(mesh, "BDM", 1)
-    DG = FunctionSpace(mesh, "DPC", 0)
-    W = BDM * DG
+    BDM = FunctionSpace(mesh, "BDMCF", 2)
+    DPC = FunctionSpace(mesh, "DPC", 1)
+    W = BDM * DPC
 
     # Define trial and test functions
     sigma, u = TrialFunctions(W)
     tau, v = TestFunctions(W)
 
     # Define source function
-    f = Function(DG).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
+    f = Function(DPC).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
 
     # Define variational form
     a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=6)
@@ -62,24 +59,6 @@ def poisson_mixed(size, parameters={}, quadrilateral=False):
     return sqrt(assemble(dot(u - f, u - f) * dx)), u, f
 
 
-# @pytest.mark.parametrize('parameters',
-#                          [{}, {'pc_type': 'fieldsplit',
-#                                'pc_fieldsplit_type': 'schur',
-#                                'ksp_type': 'preonly',
-#                                'pc_fieldsplit_schur_fact_type': 'FULL',
-#                                'fieldsplit_0_ksp_type': 'cg',
-#                                'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
-#                                'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
-#                                'fieldsplit_1_ksp_type': 'cg'}])
-# def test_poisson_mixed(parameters):
-#     """Test second-order convergence of the mixed poisson formulation."""
-#     diff = np.array([poisson_mixed(i, parameters)[0] for i in range(3, 6)])
-#     print("l2 error norms:", diff)
-#     conv = np.log2(diff[:-1] / diff[1:])
-#     print("convergence order:", conv)
-#     assert (np.array(conv) > 1.9).all()
-
-
 @pytest.mark.parametrize(('testcase', 'convrate'),
                          [((3, 6), 1.9)])
 def test_hdiv_convergence(testcase, convrate):
@@ -88,6 +67,5 @@ def test_hdiv_convergence(testcase, convrate):
     start, end = testcase
     l2err = np.zeros(end - start)
     for ii in [i + start for i in range(len(l2err))]:
-        l2err[ii - start] = poisson_mixed(ii, quadrilateral=True)[0]
-    print(l2err)
+        l2err[ii - start] = poisson_mixed(ii)[0]
     assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()

From 65bdea680ab9e3a66a8a84ef1c391c09f12de98b Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Thu, 15 Aug 2019 13:24:07 +0100
Subject: [PATCH 36/43] extrusion test for bdmc

---
 tests/extrusion/test_bdmc_extr.py | 46 +++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 tests/extrusion/test_bdmc_extr.py

diff --git a/tests/extrusion/test_bdmc_extr.py b/tests/extrusion/test_bdmc_extr.py
new file mode 100644
index 0000000000..a446c87b13
--- /dev/null
+++ b/tests/extrusion/test_bdmc_extr.py
@@ -0,0 +1,46 @@
+import numpy as np
+from firedrake import *
+import pytest
+
+
+def project_bdmc(size, degree, family):
+    mesh = ExtrudedMesh(UnitIntervalMesh(2 ** size), layers=2 ** size)
+    x = SpatialCoordinate(mesh)
+
+    fs = FunctionSpace(mesh, family, degree)
+
+    f = Function(fs)
+    v = TestFunction(fs)
+
+    expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
+                      cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
+
+    solve(inner(f-expr,v) * dx(degree=8) == 0, f)
+
+    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8)))
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmcf(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCF")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
+
+
+@pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
+                         [((3, 6), 1.9, 1),
+                          ((3, 6), 2.9, 2),
+                          ((3, 6), 3.9, 3),
+                          ((3, 6), 4.9, 4)])
+def test_bdmce(testcase, convrate, degree):
+    start, end = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = project_bdmc(ii, degree, "BDMCE")
+    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()

From 2704ce5f88c9b3b24ecd7863b75f1c68d6fda6bf Mon Sep 17 00:00:00 2001
From: cyruscycheng21 <cyruscycheng21@gmail.com>
Date: Sat, 24 Aug 2019 17:29:47 +0100
Subject: [PATCH 37/43] delete some bdmc tests

---
 .../test_poisson_mixed_no_bcs_bdmc.py         | 71 --------------
 .../test_poisson_mixed_strong_bcs_bdmc.py     | 95 -------------------
 2 files changed, 166 deletions(-)
 delete mode 100644 tests/regression/test_poisson_mixed_no_bcs_bdmc.py
 delete mode 100644 tests/regression/test_poisson_mixed_strong_bcs_bdmc.py

diff --git a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py b/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
deleted file mode 100644
index d3172a8e06..0000000000
--- a/tests/regression/test_poisson_mixed_no_bcs_bdmc.py
+++ /dev/null
@@ -1,71 +0,0 @@
-"""Solve the mixed formulation of the Laplacian described in section 2.3.1 of
-Arnold, Falk, Winther 2010 "Finite Element Exterior Calculus: From Hodge Theory
-to Numerical Stability":
-
-    sigma - grad(u) = 0
-         div(sigma) = f
-
-The corresponding weak (variational problem)
-
-    <sigma, tau> + <div(tau), u>   = 0       for all tau
-                   <div(sigma), v> = <f, v>  for all v
-
-is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
-(sigma, tau) and DPC (discontinuous Galerkin) elements of degree k - 1
-for (u, v).
-
-No strong boundary conditions are enforced. The forcing function is chosen as
-
-    -2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1]
-
-which reproduces the known analytical solution
-
-    x[0]*(1-x[0])*x[1]*(1-x[1])
-"""
-import pytest
-import numpy as np
-
-from firedrake import *
-
-
-def poisson_mixed(size, parameters={}):
-    # Create mesh
-    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
-    x = SpatialCoordinate(mesh)
-
-    # Define function spaces and mixed (product) space
-    BDM = FunctionSpace(mesh, "BDMCF", 2)
-    DPC = FunctionSpace(mesh, "DPC", 1)
-    W = BDM * DPC
-
-    # Define trial and test functions
-    sigma, u = TrialFunctions(W)
-    tau, v = TestFunctions(W)
-
-    # Define source function
-    f = Function(DPC).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
-
-    # Define variational form
-    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=6)
-    L = -f*v*dx(degree=6)
-
-    # Compute solution
-    w = Function(W)
-    solve(a == L, w, solver_parameters=parameters)
-    sigma, u = w.split()
-
-    # Analytical solution
-    f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1]))
-    return sqrt(assemble(dot(u - f, u - f) * dx)), u, f
-
-
-@pytest.mark.parametrize(('testcase', 'convrate'),
-                         [((3, 6), 1.9)])
-def test_hdiv_convergence(testcase, convrate):
-    """Test second-order convergence of the mixed poisson formulation
-    on quadrilaterals with H(div) elements."""
-    start, end = testcase
-    l2err = np.zeros(end - start)
-    for ii in [i + start for i in range(len(l2err))]:
-        l2err[ii - start] = poisson_mixed(ii)[0]
-    assert (np.log2(l2err[:-1] / l2err[1:]) > convrate).all()
diff --git a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py b/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
deleted file mode 100644
index c4bfee31df..0000000000
--- a/tests/regression/test_poisson_mixed_strong_bcs_bdmc.py
+++ /dev/null
@@ -1,95 +0,0 @@
-"""Solve the mixed formulation of the Laplacian on the unit square
-
-    sigma - grad(u) = 0
-         div(sigma) = f
-
-The corresponding weak (variational problem)
-
-    <sigma, tau> + <div(tau), u>   = 0       for all tau
-                   <div(sigma), v> = <f, v>  for all v
-
-is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
-(sigma, tau) and DG (discontinuous Galerkin) elements of degree k - 1
-for (u, v).
-
-The boundary conditions on the left and right are enforced strongly as
-
-    dot(sigma, n) = 0
-
-which corresponds to a Neumann condition du/dn = 0.
-
-The top is fixed to 42 with a Dirichlet boundary condition, which enters
-the weak formulation of the right hand side as
-
-    42*dot(tau, n)*ds
-"""
-import pytest
-from firedrake import *
-
-
-def poisson_mixed(size, parameters={}):
-    # Create mesh
-    mesh = UnitSquareMesh(2 ** size, 2 ** size, quadrilateral=True)
-    x = SpatialCoordinate(mesh)
-
-    # Define function spaces and mixed (product) space
-    BDM = FunctionSpace(mesh, "BDMCF", 1)
-    DPC = FunctionSpace(mesh, "DPC", 0)
-    W = BDM * DPC
-
-    # Define trial and test functions
-    sigma, u = TrialFunctions(W)
-    tau, v = TestFunctions(W)
-
-    # Define source function
-    f = Function(DPC).assign(0)
-
-    # Define variational form
-    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=4)
-    n = FacetNormal(mesh)
-    L = -f*v*dx + 42*dot(tau, n)*ds(4)
-
-    # Apply dot(sigma, n) == 0 on left and right boundaries strongly
-    # (corresponding to Neumann condition du/dn = 0)
-    bcs = DirichletBC(W.sub(0), Constant((0, 0)), (1, 2))
-    # Compute solution
-    w = Function(W)
-    solve(a == L, w, bcs=bcs, solver_parameters=parameters)
-    sigma, u = w.split()
-
-    # Analytical solution
-    f.interpolate(42*x[1])
-    return sqrt(assemble(dot(u - f, u - f) * dx))
-
-
-@pytest.mark.parametrize('parameters',
-                         [{},
-                          {'ksp_type': 'fgmres',
-                           'pc_type': 'fieldsplit',
-                           'pc_fieldsplit_type': 'schur',
-                           'pc_fieldsplit_schur_fact_type': 'diag',
-                           'fieldsplit_0_ksp_type': 'preonly',
-                           'fieldsplit_0_pc_type': 'icc',
-                           'fieldsplit_1_ksp_type': 'cg',
-                           'fieldsplit_1_pc_type': 'none'}])
-def test_poisson_mixed(parameters):
-    assert poisson_mixed(3, parameters) < 2e-5
-
-
-@pytest.mark.parallel(nprocs=3)
-def test_poisson_mixed_parallel_fieldsplit():
-    x = poisson_mixed(3, parameters={'ksp_type': 'fgmres',
-                                     'pc_type': 'fieldsplit',
-                                     'pc_fieldsplit_type': 'schur',
-                                     'fieldsplit_schur_fact_type': 'diag',
-                                     'fieldsplit_0_ksp_type': 'preonly',
-                                     'fieldsplit_1_ksp_type': 'cg',
-                                     'fieldsplit_0_pc_type': 'bjacobi',
-                                     'fieldsplit_0_sub_pc_type': 'icc',
-                                     'fieldsplit_1_pc_type': 'none'})
-    assert x < 2e-5
-
-
-@pytest.mark.parallel(nprocs=3)
-def test_poisson_mixed_parallel():
-    assert poisson_mixed(3) < 2e-5

From 6a76ce8671bee062c09a7972e49e9dcbf8670ff7 Mon Sep 17 00:00:00 2001
From: Patrick Farrell <patrick@pefarrell.org>
Date: Wed, 22 Jul 2020 10:25:36 +0100
Subject: [PATCH 38/43] Get lint passing

---
 tests/extrusion/test_bdmc_extr.py            | 2 +-
 tests/regression/test_bdmc.py                | 4 ++--
 tests/regression/test_bdmcf.py               | 6 ++++--
 tests/regression/test_steady_advection_2D.py | 1 +
 4 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/tests/extrusion/test_bdmc_extr.py b/tests/extrusion/test_bdmc_extr.py
index a446c87b13..8091eb808a 100644
--- a/tests/extrusion/test_bdmc_extr.py
+++ b/tests/extrusion/test_bdmc_extr.py
@@ -15,7 +15,7 @@ def project_bdmc(size, degree, family):
     expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
                       cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
 
-    solve(inner(f-expr,v) * dx(degree=8) == 0, f)
+    solve(inner(f-expr, v) * dx(degree=8) == 0, f)
 
     return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8)))
 
diff --git a/tests/regression/test_bdmc.py b/tests/regression/test_bdmc.py
index 2a3cbb56e2..6cd881feb1 100644
--- a/tests/regression/test_bdmc.py
+++ b/tests/regression/test_bdmc.py
@@ -15,9 +15,9 @@ def project_bdmc(size, degree, family):
     expr = as_vector((sin(2 * pi * x[0]) * sin(2 * pi * x[1]),
                       cos(2 * pi * x[0]) * cos(2 * pi * x[1])))
 
-    solve(inner(f-expr,v) * dx(degree=8) == 0, f)
+    solve(inner(f-expr, v) * dx(degree=8) == 0, f)
 
-    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8))) 
+    return np.sqrt(assemble(inner(f-expr, f-expr) * dx(degree=8)))
 
 
 @pytest.mark.parametrize(('testcase', 'convrate', 'degree'),
diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
index 1e0904c502..9ec2887d38 100644
--- a/tests/regression/test_bdmcf.py
+++ b/tests/regression/test_bdmcf.py
@@ -21,14 +21,16 @@
 
 # Compute solution
 w = Function(W)
-solve(a == L, w, solver_parameters={'pc_type': 'fieldsplit',
+sp = {'pc_type': 'fieldsplit',
       'pc_fieldsplit_type': 'schur',
       'ksp_type': 'preonly',
       'pc_fieldsplit_schur_fact_type': 'FULL',
       'fieldsplit_0_ksp_type': 'cg',
       'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
       'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
-      'fieldsplit_1_ksp_type': 'cg'})
+      'fieldsplit_1_ksp_type': 'cg'}
+
+solve(a == L, w, solver_parameters=sp)
 # solve(a == L, w)
 sigma, u = w.split()
 
diff --git a/tests/regression/test_steady_advection_2D.py b/tests/regression/test_steady_advection_2D.py
index 8fd13ca596..f37f54ae15 100644
--- a/tests/regression/test_steady_advection_2D.py
+++ b/tests/regression/test_steady_advection_2D.py
@@ -21,6 +21,7 @@ def DGDPC0(request, mesh):
     else:
         return FunctionSpace(mesh, request.param, 0)
 
+
 @pytest.fixture(scope='module', params=["DG", "DPC"])
 def DGDPC1(request, mesh):
     if mesh.ufl_cell() == triangle:

From b8ec6ad0c1beeddfad2cce48d79d6bae570b626d Mon Sep 17 00:00:00 2001
From: Patrick Farrell <patrick@pefarrell.org>
Date: Wed, 22 Jul 2020 10:27:24 +0100
Subject: [PATCH 39/43] Remove dud test

---
 tests/regression/test_bdmcf.py | 49 ----------------------------------
 1 file changed, 49 deletions(-)
 delete mode 100644 tests/regression/test_bdmcf.py

diff --git a/tests/regression/test_bdmcf.py b/tests/regression/test_bdmcf.py
deleted file mode 100644
index 9ec2887d38..0000000000
--- a/tests/regression/test_bdmcf.py
+++ /dev/null
@@ -1,49 +0,0 @@
-from firedrake import *
-import matplotlib.pyplot as plt
-
-mesh = UnitSquareMesh(16, 16, quadrilateral=True)
-x = SpatialCoordinate(mesh)
-
-BDM = FunctionSpace(mesh, "BDMCF", 2)
-DPC = FunctionSpace(mesh, "DPC", 1)
-W = BDM * DPC
-
-# Define trial and test functions
-sigma, u = TrialFunctions(W)
-tau, v = TestFunctions(W)
-
-# Define source function
-f = Function(DPC).interpolate(-2*(x[0]-1)*x[0] - 2*(x[1]-1)*x[1])
-
-# Define variational form
-a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=10)
-L = -f*v*dx(degree=10)
-
-# Compute solution
-w = Function(W)
-sp = {'pc_type': 'fieldsplit',
-      'pc_fieldsplit_type': 'schur',
-      'ksp_type': 'preonly',
-      'pc_fieldsplit_schur_fact_type': 'FULL',
-      'fieldsplit_0_ksp_type': 'cg',
-      'fieldsplit_0_pc_factor_shift_type': 'INBLOCKS',
-      'fieldsplit_1_pc_factor_shift_type': 'INBLOCKS',
-      'fieldsplit_1_ksp_type': 'cg'}
-
-solve(a == L, w, solver_parameters=sp)
-# solve(a == L, w)
-sigma, u = w.split()
-
-# Analytical solution
-f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1]))
-
-l2err = errornorm(f, u, norm_type="l2")
-print(l2err)
-
-plt.figure(1)
-plot(u)
-
-plt.figure(2)
-plot(f)
-
-plt.show()

From dfc0bf8a9f17fa32af6c8d760c88c52e4879aafa Mon Sep 17 00:00:00 2001
From: Patrick Farrell <patrick@pefarrell.org>
Date: Wed, 22 Jul 2020 11:05:32 +0100
Subject: [PATCH 40/43] Initial effort at Riesz map tests for BDMC[E/F]. Not
 converging for lowest order.

---
 tests/regression/test_bdmc_riesz_map.py | 89 +++++++++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 tests/regression/test_bdmc_riesz_map.py

diff --git a/tests/regression/test_bdmc_riesz_map.py b/tests/regression/test_bdmc_riesz_map.py
new file mode 100644
index 0000000000..ffe6089fc4
--- /dev/null
+++ b/tests/regression/test_bdmc_riesz_map.py
@@ -0,0 +1,89 @@
+from firedrake import *
+import pytest
+import numpy
+
+
+@pytest.fixture(scope='module', params=["div", "curl"])
+def problem(request):
+    return request.param
+
+
+# FIXME: MMS does not converge for degree 1
+# @pytest.fixture(scope='module', params=[1, 2, 3])
+@pytest.fixture(scope='module', params=[2, 3])
+def degree(request):
+    return request.param
+
+
+@pytest.fixture(scope='module', params=[2])
+def dim(request):
+    return request.param
+
+
+sp = {"snes_type": "ksponly",
+      "ksp_type": "preonly",
+      "pc_type": "lu",
+      "pc_factor_mat_solver_type": "mumps",
+      "mat_mumps_icntl_14": 200}
+
+
+def error(N, problem, degree, dim):
+    mesh = UnitSquareMesh(N, N, quadrilateral=True)
+    if dim == 3:
+        mesh = ExtrudedMesh(mesh, N)
+
+    if problem == "div":
+        op = div
+        family = "BDMCF"
+    elif problem == "curl":
+        op = curl
+        family = "BDMCE"
+    else:
+        raise ValueError
+
+    V = FunctionSpace(mesh, family, degree)
+    u = Function(V)
+    v = TestFunction(V)
+
+    k = 1
+    if dim == 2:
+        (x, y) = SpatialCoordinate(mesh)
+        u_ex = as_vector([sin(2*pi*x) * cos(2*pi*y),
+                          x * (1-x) * y * (1-y)])
+    else:
+        (x, y, z) = SpatialCoordinate(mesh)
+        u_ex = as_vector([sin(2*pi*k*x) * cos(2*pi*k*y) * sin(4*pi*z),
+                          x * (1-x) * z * (1-z),
+                          y * (1-y) * z * (1-z)])
+
+    if problem == "div":
+        f = u_ex - grad(div(u_ex))
+    else:
+        f = u_ex + curl(curl(u_ex))
+
+    F = inner(op(u), op(v))*dx + inner(u, v)*dx - inner(f, v)*dx
+    bc = DirichletBC(V, project(u_ex, V, solver_parameters=sp), "on_boundary")
+
+    solve(F == 0, u, bc, solver_parameters=sp)
+
+    err = errornorm(u_ex, u, "L2")
+    return err
+
+
+def test_bdmc_riesz_map(problem, degree, dim):
+
+    if dim == 2:
+        Ns = [10, 20, 40]
+    else:
+        Ns = [2, 4, 8]
+
+    errors = []
+    for N in Ns:
+        errors.append(error(N, problem, degree, dim))
+
+    convergence_orders = lambda x: numpy.log2(numpy.array(x)[:-1] / numpy.array(x)[1:])
+    conv = convergence_orders(errors)
+    print("errors: ", errors)
+    print("convergence order: ", conv)
+
+    assert numpy.allclose(conv, degree, atol=0.1)

From 18814f77618568b87042d1444fda8b6acd7eb1fa Mon Sep 17 00:00:00 2001
From: Patrick Farrell <patrick@pefarrell.org>
Date: Wed, 22 Jul 2020 11:10:36 +0100
Subject: [PATCH 41/43] Simplify test; AAF/AAE elements not implemented, so no
 need to test them

---
 tests/regression/test_bdmc_riesz_map.py | 33 ++++++-------------------
 1 file changed, 7 insertions(+), 26 deletions(-)

diff --git a/tests/regression/test_bdmc_riesz_map.py b/tests/regression/test_bdmc_riesz_map.py
index ffe6089fc4..157a4e50de 100644
--- a/tests/regression/test_bdmc_riesz_map.py
+++ b/tests/regression/test_bdmc_riesz_map.py
@@ -15,11 +15,6 @@ def degree(request):
     return request.param
 
 
-@pytest.fixture(scope='module', params=[2])
-def dim(request):
-    return request.param
-
-
 sp = {"snes_type": "ksponly",
       "ksp_type": "preonly",
       "pc_type": "lu",
@@ -27,10 +22,8 @@ def dim(request):
       "mat_mumps_icntl_14": 200}
 
 
-def error(N, problem, degree, dim):
+def error(N, problem, degree):
     mesh = UnitSquareMesh(N, N, quadrilateral=True)
-    if dim == 3:
-        mesh = ExtrudedMesh(mesh, N)
 
     if problem == "div":
         op = div
@@ -45,16 +38,9 @@ def error(N, problem, degree, dim):
     u = Function(V)
     v = TestFunction(V)
 
-    k = 1
-    if dim == 2:
-        (x, y) = SpatialCoordinate(mesh)
-        u_ex = as_vector([sin(2*pi*x) * cos(2*pi*y),
-                          x * (1-x) * y * (1-y)])
-    else:
-        (x, y, z) = SpatialCoordinate(mesh)
-        u_ex = as_vector([sin(2*pi*k*x) * cos(2*pi*k*y) * sin(4*pi*z),
-                          x * (1-x) * z * (1-z),
-                          y * (1-y) * z * (1-z)])
+    (x, y) = SpatialCoordinate(mesh)
+    u_ex = as_vector([sin(2*pi*x) * cos(2*pi*y),
+                      x * (1-x) * y * (1-y)])
 
     if problem == "div":
         f = u_ex - grad(div(u_ex))
@@ -70,16 +56,11 @@ def error(N, problem, degree, dim):
     return err
 
 
-def test_bdmc_riesz_map(problem, degree, dim):
-
-    if dim == 2:
-        Ns = [10, 20, 40]
-    else:
-        Ns = [2, 4, 8]
+def test_bdmc_riesz_map(problem, degree):
 
     errors = []
-    for N in Ns:
-        errors.append(error(N, problem, degree, dim))
+    for N in [10, 20, 40]:
+        errors.append(error(N, problem, degree))
 
     convergence_orders = lambda x: numpy.log2(numpy.array(x)[:-1] / numpy.array(x)[1:])
     conv = convergence_orders(errors)

From eac78b8b94e31480337bcdbfd34029a17bf875cb Mon Sep 17 00:00:00 2001
From: Tom Bendall <t.bendall15@imperial.ac.uk>
Date: Tue, 26 Jan 2021 19:16:15 +0000
Subject: [PATCH 42/43] update riesz map test for BDMCE/F elements

---
 tests/regression/test_bdmc_riesz_map.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/regression/test_bdmc_riesz_map.py b/tests/regression/test_bdmc_riesz_map.py
index 157a4e50de..7b8c17809a 100644
--- a/tests/regression/test_bdmc_riesz_map.py
+++ b/tests/regression/test_bdmc_riesz_map.py
@@ -8,9 +8,7 @@ def problem(request):
     return request.param
 
 
-# FIXME: MMS does not converge for degree 1
-# @pytest.fixture(scope='module', params=[1, 2, 3])
-@pytest.fixture(scope='module', params=[2, 3])
+@pytest.fixture(scope='module', params=[1, 2, 3])
 def degree(request):
     return request.param
 
@@ -67,4 +65,6 @@ def test_bdmc_riesz_map(problem, degree):
     print("errors: ", errors)
     print("convergence order: ", conv)
 
-    assert numpy.allclose(conv, degree, atol=0.1)
+    tol = 0.11
+
+    assert (conv > (degree + 1) - tol).all()

From f0fff2051e55660e253c4cfbd572954a9f118065 Mon Sep 17 00:00:00 2001
From: Jack Betteridge <J.Betteridge@imperial.ac.uk>
Date: Wed, 20 Sep 2023 16:37:09 +0100
Subject: [PATCH 43/43] Revert package branches

---
 .github/workflows/build.yml | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 0b6c7c4b5c..291634e4f6 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -74,9 +74,6 @@ jobs:
             --torch \
             --netgen \
             --slepc \
-            --package-branch fiat bdmc \
-            --package-branch tsfc bdmc \
-            --package-branch FInAT bdmc \
             --install thetis \
             --install gusto \
             --install icepack \