From 2dc873ccc1e623679e83da0842a26ac510c4dd8f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 13 Jun 2025 12:15:46 +0100 Subject: [PATCH 1/4] fix normals and allow fuse use to be turned on/off with env vars --- fuse/cells.py | 38 ++++++++++++++++++++++++++++++++++---- pyproject.toml | 4 ++++ 2 files changed, 38 insertions(+), 4 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 14d4f4d..b5ee47a 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -638,6 +638,21 @@ def basis_vectors(self, return_coords=True, entity=None): basis_vecs.append((v, v_0)) return basis_vecs + def compute_normal(self, facet_id): + # computes normal to a facet of codimension 1 + facet = self.d_entities(self.dimension - 1)[facet_id] + entityBasis = np.array(facet.basis_vectors()) + cellEntityBasis = np.array(self.basis_vectors(entity=facet)) + basis = np.matmul(entityBasis, cellEntityBasis) + + if facet.dimension == 1: + result = np.matmul(basis, np.array([[0, -1], [1, 0]])) + elif facet.dimension == 2: + result = np.cross(basis[0], basis[1]) + else: + raise NotImplementedError("Normals in higher dimensions") + return result.squeeze() + def to_tikz(self, show=True, scale=3): tikz_commands = [] if show: @@ -956,6 +971,9 @@ def get_facet_element(self): dimension = self.get_spatial_dimension() return self.construct_subelement(dimension - 1) + def compute_normal(self, facet_id): + return self.fe_cell.compute_normal(facet_id) + class CellComplexToFiatTensorProduct(FiatTensorProductCell): """ @@ -1033,6 +1051,9 @@ def get_facet_element(self): def get_dimension(self): return self.get_spatial_dimension() + def compute_normal(self, facet_id): + return self.fe_cell.compute_normal(facet_id) + class CellComplexToUFL(Cell): """ @@ -1091,24 +1112,33 @@ def reconstruct(self, **kwargs): def constructCellComplex(name): + # import ufl + # return ufl.Cell(name) + + import os + fuse = os.environ.get("FIREDRAKE_USE_FUSE", "False") + # breakpoint() if name == "vertex": return Point(0).to_ufl(name) elif name == "interval": return Point(1, [Point(0), Point(0)], vertex_num=2).to_ufl(name) elif name == "triangle": + if fuse == "False": + return ufc_triangle().to_ufl(name) return polygon(3).to_ufl(name) - # return ufc_triangle().to_ufl(name) elif name == "quadrilateral": interval = Point(1, [Point(0), Point(0)], vertex_num=2) return TensorProductPoint(interval, interval).flatten().to_ufl(name) # return ufc_quad().to_ufl(name) # return polygon(4).to_ufl(name) elif name == "tetrahedron": - # return ufc_tetrahedron().to_ufl(name) + if fuse == "False": + return ufc_tetrahedron().to_ufl(name) return make_tetrahedron().to_ufl(name) elif name == "hexahedron": - import warnings - warnings.warn("Hexahedron unimplemented in Fuse") + if fuse == "True": + import warnings + warnings.warn("Hexahedron unimplemented in Fuse") import ufl return ufl.Cell(name) elif "*" in name: diff --git a/pyproject.toml b/pyproject.toml index 542ced8..604c8ac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,7 @@ docs = [ ] dev = [ "pytest", + "pytest-env", "pytest-mock", "coverage", "flake8", @@ -48,6 +49,9 @@ dev = [ testpaths = [ "tests", ] +env = [ + "FIREDRAKE_USE_FUSE=True" +] [tool.coverage.run] include=[ From 81e1626df01f2f64a4d0671d53f862ae3cb12171 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 17 Sep 2025 16:18:59 +0100 Subject: [PATCH 2/4] scalar convergence test --- test/test_convert_to_fiat.py | 44 ++++++++++++++++++++++++++++-------- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index a08a74a..8ee4837 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -551,7 +551,7 @@ def test_non_tensor_quad(): create_cg1_quad() -def project(U, mesh, func): +def project_func(U, mesh, func): f = assemble(interpolate(func, U)) out = Function(U) @@ -581,10 +581,10 @@ def test_project(elem_gen, elem_code, deg): mesh = UnitTriangleMesh() U = FunctionSpace(mesh, elem_code, deg) - assert np.allclose(project(U, mesh, Constant(1)), 0, rtol=1e-5) + assert np.allclose(project_func(U, mesh, Constant(1)), 0, rtol=1e-5) U = FunctionSpace(mesh, elem.to_ufl()) - assert np.allclose(project(U, mesh, Constant(1)), 0, rtol=1e-5) + assert np.allclose(project_func(U, mesh, Constant(1)), 0, rtol=1e-5) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_dg1_tet, "DG", 1)]) @@ -595,13 +595,13 @@ def test_project_3d(elem_gen, elem_code, deg): mesh = UnitCubeMesh(3, 3, 3) U = FunctionSpace(mesh, elem_code, deg) - assert np.allclose(project(U, mesh, Constant(1)), 0, rtol=1e-5) + assert np.allclose(project_func(U, mesh, Constant(1)), 0, rtol=1e-5) U = FunctionSpace(mesh, elem.to_ufl()) - assert np.allclose(project(U, mesh, Constant(1)), 0, rtol=1e-5) + assert np.allclose(project_func(U, mesh, Constant(1)), 0, rtol=1e-5) - -@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [pytest.param(create_dg1_tet, "DG", 1, 0.8, marks=pytest.mark.xfail(reason="DG on tets - check test written correctly"))]) +# marks=pytest.mark.xfail(reason="DG on tets - check test written correctly") +@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(create_dg1_tet, "CG", 3, 0.8)]) def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate): cell = make_tetrahedron() elem = elem_gen(cell) @@ -615,11 +615,11 @@ def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate): x = SpatialCoordinate(mesh) V = FunctionSpace(mesh, elem_code, deg) - res1 = project(V, mesh, function(x)) + res1 = project_func(V, mesh, function(x)) diff2[i - 1] = res1 V2 = FunctionSpace(mesh, elem.to_ufl()) - res2 = project(V2, mesh, function(x)) + res2 = project_func(V2, mesh, function(x)) diff[i - 1] = res2 assert np.allclose(res1, res2) @@ -635,3 +635,29 @@ def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate): assert (np.array(conv1) > conv_rate).all() assert (np.array(conv2) > conv_rate).all() + + +# @pytest.mark.parametrize(('testcase', 'convrate'), +# [(("CG", 1), 1.5), (("CG", 2), 2.6), +# (("DG", 0), 0.9), (("DG", 1), 1.7)]) +def test_scalar_convergence(): + convrate = 1.7 + elem = create_dg1_tet(make_tetrahedron()) + l2err = np.zeros(2) + for ii in range(len(l2err)): + mesh = UnitCubeMesh(2 ** ii, 2 ** ii, 2 ** ii) + + fspace = FunctionSpace(mesh, elem.to_ufl()) + exactfspace = FunctionSpace(mesh, "Lagrange", 3) + + u = TrialFunction(fspace) + v = TestFunction(fspace) + + x, y, z = SpatialCoordinate(mesh) + expr = x*x*y*z + exact = project(expr, exactfspace) + + out = Function(fspace) + solve(inner(u, v)*dx == inner(exact, v)*dx, out) + l2err[ii] = sqrt(assemble(inner((out-exact), (out-exact))*dx)) + assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all() From 989973ce4284fb6b463c8919d01333d4019fe49a Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 25 Sep 2025 16:50:41 +0100 Subject: [PATCH 3/4] add env variable to conftest --- conftest.py | 7 ++++++- fuse/cells.py | 1 - 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/conftest.py b/conftest.py index c8c4e00..fb65b29 100644 --- a/conftest.py +++ b/conftest.py @@ -1,4 +1,9 @@ import pytest +import os + +@pytest.fixture(scope="session", autouse=True) +def set_env(): + os.environ["FIREDRAKE_USE_FUSE"] = "True" def pytest_addoption(parser): parser.addoption( @@ -6,4 +11,4 @@ def pytest_addoption(parser): action="store_true", default=False, help="Run tests that require a cleared cache", - ) \ No newline at end of file + ) diff --git a/fuse/cells.py b/fuse/cells.py index b5ee47a..1f60b7b 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -1117,7 +1117,6 @@ def constructCellComplex(name): import os fuse = os.environ.get("FIREDRAKE_USE_FUSE", "False") - # breakpoint() if name == "vertex": return Point(0).to_ufl(name) elif name == "interval": From 1eb04d425f4fd5ec83eed5bbb4ed7b0d5f5044a1 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 25 Sep 2025 16:54:12 +0100 Subject: [PATCH 4/4] add xfail and lint --- test/test_convert_to_fiat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 8ee4837..1081559 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -600,8 +600,8 @@ def test_project_3d(elem_gen, elem_code, deg): U = FunctionSpace(mesh, elem.to_ufl()) assert np.allclose(project_func(U, mesh, Constant(1)), 0, rtol=1e-5) -# marks=pytest.mark.xfail(reason="DG on tets - check test written correctly") -@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(create_dg1_tet, "CG", 3, 0.8)]) + +@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [pytest.param(create_dg1_tet, "CG", 3, 0.8, marks=pytest.mark.xfail(reason="DG on tets - check test written correctly"))]) def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate): cell = make_tetrahedron() elem = elem_gen(cell)