Test integral variant for Regge/HHJ (#3839)

* DO NOT MERGE

* Test VOM with point variant

* Update .github/workflows/build.yml

* Test Regge point variant

* add continuity test

* Test H(curl div) element (#3843)

* Test H(curl div) element

* test Hu-Zhang

* checkout ufl master branch

* checkout master branches

tests/regression/test_interpolate_cross_mesh.py

@@ -215,7 +215,7 @@ def parameters(request):
         expected = np.asarray(
             [np.outer(coords[i], coords[i]) for i in range(len(coords))]
         )
-        V_src = FunctionSpace(m_src, "Regge", 2)  # Not point evaluation nodes
+        V_src = FunctionSpace(m_src, "Regge", 2, variant="point")  # Not point evaluation nodes
         V_dest = TensorFunctionSpace(m_dest, "CG", 4)
         V_dest_2 = TensorFunctionSpace(m_dest, "DQ", 2)
     elif request.param == "spheresphere":

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

@@ -0,0 +1,51 @@
+from firedrake import *
+import pytest
+
+
+rg = RandomGenerator(PCG64(seed=1))
+
+
+@pytest.fixture(params=("square", "cube"))
+def mesh(request):
+    if request.param == "square":
+        return UnitSquareMesh(2, 2)
+    elif request.param == "cube":
+        return UnitCubeMesh(2, 2, 2)
+    else:
+        raise ValueError("Unrecognized mesh type")
+
+
+@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)
+
+    n = FacetNormal(mesh)
+
+    space = V.ufl_element().sobolev_space
+    if space == HDivDiv:
+        utrace = dot(n, dot(u, n))
+    elif space == HEin:
+        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))
+    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

tests/vertexonly/test_interpolation_from_parent.py

@@ -102,7 +102,7 @@ def vfs(request, parentmesh):
 
 @pytest.fixture(params=[("CG", 2, TensorFunctionSpace),
                         ("BDM", 2, VectorFunctionSpace),
-                        ("Regge", 2, FunctionSpace)],
+                        ("Regge", 2, lambda *args: FunctionSpace(*args, variant="point"))],
                 ids=lambda x: f"{x[2].__name__}({x[0]}{x[1]})")
 def tfs(request, parentmesh):
     family = request.param[0]