Test H(curl div) element (#3843)

* Test H(curl div) element

* test Hu-Zhang

.github/workflows/build.yml

@@ -86,6 +86,7 @@ jobs:
             --package-branch fiat pbrubeck/h1curl \
             --package-branch finat pbrubeck/hhj \
             --package-branch tsfc pbrubeck/hhj \
+            --package-branch ufl pbrubeck/hcurldiv \
             || (cat firedrake-install.log && /bin/false)
       - name: Install test dependencies
         run: |

tests/regression/test_stress_els.py → tests/regression/test_stress_elements.py

@@ -7,19 +7,21 @@
 convergence_orders = lambda x: -np.log2(relative_magnitudes(x))
 
 
-@pytest.fixture(scope='module', params=["conforming", "nonconforming", "macro"])
+@pytest.fixture(params=["macro", "nonconforming", "conforming", "high-order"])
 def stress_element(request):
-    if request.param == "conforming":
-        return FiniteElement("AWc", triangle, 3)
+    if request.param == "macro":
+        return FiniteElement("JM", triangle, 1)
     elif request.param == "nonconforming":
         return FiniteElement("AWnc", triangle, 2)
-    elif request.param == "macro":
-        return FiniteElement("JM", triangle, 1)
+    elif request.param == "conforming":
+        return FiniteElement("AWc", triangle, 3)
+    elif request.param == "high-order":
+        return FiniteElement("HZ", triangle, 3)
     else:
         raise ValueError("Unknown family")
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture
 def mesh_hierarchy(request):
     N_base = 2
     mesh = UnitSquareMesh(N_base, N_base)
@@ -64,7 +66,10 @@ def epsilon(u):
     l2_u = []
     l2_sigma = []
     l2_div_sigma = []
-    displacement_element = VectorElement("DG", cell=mesh.ufl_cell(), degree=1, variant="integral")
+    sdegree = stress_element.degree()
+    sfamily = stress_element.family()
+    vdegree = sdegree-1 if sfamily == "Hu-Zhang" else 1
+    displacement_element = VectorElement("DG", cell=mesh.ufl_cell(), degree=vdegree, variant="integral")
     element = MixedElement([stress_element, displacement_element])
     for msh in mesh_hierarchy[1:]:
         x, y = SpatialCoordinate(msh)
@@ -126,20 +131,18 @@ def epsilon(u):
         l2_sigma.append(error_sigma)
         l2_div_sigma.append(error_div_sigma)
 
-    if stress_element.family().startswith("Conforming"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 2.9
-        assert min(convergence_orders(l2_div_sigma)) > 1.9
-    elif stress_element.family().startswith("Nonconforming"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 1
-        assert min(convergence_orders(l2_div_sigma)) > 1.9
-    elif stress_element.family().startswith("Johnson-Mercier"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 1
-        assert min(convergence_orders(l2_div_sigma)) > 0.9
+    sdegree = V[0].ufl_element().degree()
+    vdegree = V[1].ufl_element().degree()
+    if stress_element.family().startswith("Nonconforming"):
+        expected_rates = (vdegree+1, sdegree-1, sdegree)
+    elif stress_element.family().startswith("Conforming"):
+        expected_rates = (vdegree+1, sdegree, sdegree-1)
     else:
-        raise ValueError("Don't know what the convergence should be")
+        expected_rates = (vdegree+1, sdegree+1, sdegree)
+
+    assert min(convergence_orders(l2_u[1:])) > expected_rates[0] * 0.9
+    assert min(convergence_orders(l2_sigma[1:])) > expected_rates[1] * 0.9
+    assert min(convergence_orders(l2_div_sigma[1:])) > expected_rates[2] * 0.9
 
 
 def test_mass_conditioning(stress_element, mesh_hierarchy):

tests/regression/test_tensor_elements.py

@@ -15,25 +15,37 @@ def mesh(request):
         raise ValueError("Unrecognized mesh type")
 
 
-@pytest.mark.parametrize("degree", (0, 1, 2))
-@pytest.mark.parametrize("family", ("Regge", "HHJ"))
+@pytest.mark.parametrize("family, degree", [
+                         ("Regge", 0),
+                         ("Regge", 1),
+                         ("Regge", 2),
+                         ("HHJ", 0),
+                         ("HHJ", 1),
+                         ("HHJ", 2),
+                         ("GLS", 1),
+                         ("GLS", 2),
+                         ])
 def test_tensor_continuity(mesh, family, degree):
     V = FunctionSpace(mesh, family, degree)
     u = rg.beta(V, 1.0, 2.0)
 
-    dim = mesh.topological_dimension()
     n = FacetNormal(mesh)
 
     space = V.ufl_element().sobolev_space
     if space == HDivDiv:
         utrace = dot(n, dot(u, n))
     elif space == HEin:
-        if dim == 2:
-            t = dot(as_matrix([[0, -1], [1, 0]]), n)
-            utrace = dot(t, dot(u, t))
+        if mesh.topological_dimension() == 2:
+            t = perp(n)
         else:
             t = as_matrix([[0, n[2], -n[1]], [-n[2], 0, n[0]], [n[1], -n[0], 0]])
-            utrace = dot(t, dot(u, t.T))
+        utrace = dot(t, dot(u, t))
+    else:
+        if mesh.topological_dimension() == 2:
+            t = perp(n)
+            utrace = dot(t, dot(u, n))
+        else:
+            utrace = cross(n, dot(u, n))
 
     ujump = utrace('+') - utrace('-')
     assert assemble(inner(ujump, ujump)*dS) < 1E-10