Integration with Fuse

.github/workflows/build.yml

@@ -87,6 +87,8 @@ jobs:
             --install defcon \
             --install gadopt \
             --install asQ \
+            --package-branch fiat indiamai/integrate_fuse \
+            --package-branch fuse main \
             || (cat firedrake-install.log && /bin/false)
 
       - name: Install test dependencies

firedrake/mesh.py

@@ -37,6 +37,7 @@
 from firedrake.adjoint_utils import MeshGeometryMixin
 from pyadjoint import stop_annotating
 import gem
+from fuse import constructCellComplex
 
 try:
     import netgen
@@ -1216,15 +1217,15 @@ def _ufl_cell(self):
 
         # TODO: this needs to be updated for mixed-cell meshes.
         nfacets = self._comm.allreduce(nfacets, op=MPI.MAX)
-
         # Note that the geometric dimension of the cell is not set here
         # despite it being a property of a UFL cell. It will default to
         # equal the topological dimension.
         # Firedrake mesh topologies, by convention, which specifically
         # represent a mesh topology (as here) have geometric dimension
         # equal their topological dimension. This is reflected in the
         # corresponding UFL mesh.
-        return ufl.Cell(_cells[tdim][nfacets])
+        # return ufl.Cell(_cells[tdim][nfacets])
+        return constructCellComplex(_cells[tdim][nfacets])
 
     @utils.cached_property
     def _ufl_mesh(self):
@@ -1262,6 +1263,7 @@ def cell_closure(self):
 
         cell = self.ufl_cell()
         assert tdim == cell.topological_dimension()
+
         if self.submesh_parent is not None:
             return dmcommon.submesh_create_cell_closure_cell_submesh(plex,
                                                                      self.submesh_parent.topology_dm,
@@ -1273,10 +1275,16 @@ def cell_closure(self):
             entity_per_cell = np.zeros(len(topology), dtype=IntType)
             for d, ents in topology.items():
                 entity_per_cell[d] = len(ents)
-
             return dmcommon.closure_ordering(plex, vertex_numbering,
                                              cell_numbering, entity_per_cell)
-
+        # elif hasattr(cell, "to_fiat"):
+        #       TODO
+        #     topology = cell.to_fiat().topology
+        #     entity_per_cell = np.zeros(len(topology), dtype=IntType)
+        #     for d, ents in topology.items():
+        #         entity_per_cell[d] = len(ents)
+        #     return dmcommon.closure_ordering(plex, vertex_numbering,
+        #                                      cell_numbering, entity_per_cell)
         elif cell.cellname() == "quadrilateral":
             from firedrake_citations import Citations
             Citations().register("Homolya2016")
@@ -1332,7 +1340,6 @@ def merge_ids(x, y, datatype):
             op.Free()
         else:
             unique_markers = None
-
         local_facet_number, facet_cell = \
             dmcommon.facet_numbering(dm, kind, facets,
                                      self._cell_numbering,
@@ -1955,7 +1962,7 @@ def _mark_entity_classes(self):
 
     @utils.cached_property
     def _ufl_cell(self):
-        return ufl.Cell(_cells[0][0])
+        return constructCellComplex(_cells[0][0])
 
     @utils.cached_property
     def _ufl_mesh(self):

pyproject.toml

@@ -39,6 +39,7 @@ dependencies = [
   "pyadjoint-ad @ git+https://github.com/dolfin-adjoint/pyadjoint.git",
   "loopy @ git+https://github.com/firedrakeproject/loopy.git@main",
   "libsupermesh @ git+https://github.com/firedrakeproject/libsupermesh.git",
+  "fuse-element @ git+https://github.com/indiamai/fuse.git",
 ]
 classifiers = [
   "Development Status :: 5 - Production/Stable",

requirements-git.txt

@@ -5,3 +5,4 @@ git+https://github.com/firedrakeproject/loopy.git@main#egg=loopy
 git+https://github.com/firedrakeproject/pytest-mpi.git@main#egg=pytest-mpi
 git+https://github.com/firedrakeproject/petsc.git@firedrake#egg=petsc
 git+https://github.com/firedrakeproject/libsupermesh.git#egg=libsupermesh
+git+https://github.com/indiamai/fuse.git#egg=fuse

tests/firedrake/regression/test_fuse.py

@@ -0,0 +1,24 @@
+from test_helmholtz import helmholtz
+from test_poisson_strong_bcs import run_test
+import pytest
+import numpy as np
+
+
+
+@pytest.mark.parametrize(['params', 'degree', 'quadrilateral'],
+                         [(p, d, q)
+                          for p in [{}, {'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu'}]
+                          for d in (1, 2, 3)
+                          for q in [False, True]])
+def test_poisson_analytic(params, degree, quadrilateral):
+    assert (run_test(2, degree, parameters=params, quadrilateral=False) < 1.e-9)
+
+@pytest.mark.parametrize(['conv_num', 'degree'],
+                         [(p, d)
+                          for p,d in zip([1.8, 2.8, 3.8],[1, 2, 3])])
+def test_helmholtz(conv_num, degree):
+    diff = np.array([helmholtz(i, degree=degree)[0] for i in range(3, 6)])
+    print("l2 error norms:", diff)
+    conv = np.log2(diff[:-1] / diff[1:])
+    print("convergence order:", conv)
+    assert (np.array(conv) > conv_num).all()