test Hu-Zhang

tests/regression/test_stress_els.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=["conforming", "nonconforming", "macro", "high-order"])
 def stress_element(request):
     if request.param == "conforming":
         return FiniteElement("AWc", triangle, 3)
     elif request.param == "nonconforming":
         return FiniteElement("AWnc", triangle, 2)
     elif request.param == "macro":
         return FiniteElement("JM", triangle, 1)
+    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,21 +131,23 @@ def epsilon(u):
         l2_sigma.append(error_sigma)
         l2_div_sigma.append(error_div_sigma)
 
+    sdegree = V[0].ufl_element().degree()
+    vdegree = V[1].ufl_element().degree()
     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
+        rates = (vdegree+1, sdegree, sdegree-1)
     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
+        rates = (vdegree+1, sdegree-1, sdegree)
     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
+        rates = (vdegree+1, sdegree, sdegree)
+    elif stress_element.family().startswith("Hu-Zhang"):
+        rates = (vdegree+1, sdegree+1, sdegree)
     else:
         raise ValueError("Don't know what the convergence should be")
 
+    assert min(convergence_orders(l2_u)) > rates[0] * 0.9
+    assert min(convergence_orders(l2_sigma)) > rates[1] * 0.9
+    assert min(convergence_orders(l2_div_sigma)) > rates[2] * 0.9
+
 
 def test_mass_conditioning(stress_element, mesh_hierarchy):
     mass_cond = []