From 3053928436bfd0f6528b1d22ff8a1a2526c7ff42 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 19 Nov 2024 23:51:35 +0000 Subject: [PATCH] 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 --- .../regression/test_interpolate_cross_mesh.py | 2 +- ..._stress_els.py => test_stress_elements.py} | 43 ++++++++-------- tests/regression/test_tensor_elements.py | 51 +++++++++++++++++++ .../test_interpolation_from_parent.py | 2 +- 4 files changed, 76 insertions(+), 22 deletions(-) rename tests/regression/{test_stress_els.py => test_stress_elements.py} (81%) create mode 100644 tests/regression/test_tensor_elements.py diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py index 94252c3df6..ef9fa28769 100644 --- a/tests/regression/test_interpolate_cross_mesh.py +++ b/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": diff --git a/tests/regression/test_stress_els.py b/tests/regression/test_stress_elements.py similarity index 81% rename from tests/regression/test_stress_els.py rename to tests/regression/test_stress_elements.py index 16ca97bc58..35cd8bb521 100755 --- a/tests/regression/test_stress_els.py +++ b/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): diff --git a/tests/regression/test_tensor_elements.py b/tests/regression/test_tensor_elements.py new file mode 100644 index 0000000000..6c0a6f99ff --- /dev/null +++ b/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 diff --git a/tests/vertexonly/test_interpolation_from_parent.py b/tests/vertexonly/test_interpolation_from_parent.py index 7016da12c6..1ac9a0a393 100644 --- a/tests/vertexonly/test_interpolation_from_parent.py +++ b/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]