Normals +Env vars

conftest.py

@@ -1,9 +1,14 @@
 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(
         "--run-cleared",
         action="store_true",
         default=False,
         help="Run tests that require a cleared cache",
-    )
+    )

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,32 @@ def reconstruct(self, **kwargs):
 
 
 def constructCellComplex(name):
+    # import ufl
+    # return ufl.Cell(name)
+
+    import os
+    fuse = os.environ.get("FIREDRAKE_USE_FUSE", "False")
     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:

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=[

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"))])
+@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)
@@ -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()