Indiamai/fuse quads

.github/workflows/core.yml

@@ -210,6 +210,10 @@ jobs:
             --no-binary h5py \
             --extra-index-url https://download.pytorch.org/whl/cpu \
             "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]"
+
+          pip install -I "fenics-ufl @ git+https://github.com/firedrakeproject/ufl.git@indiamai/integrate-fuse"
+          pip install -I "firedrake-fiat @ git+https://github.com/firedrakeproject/fiat.git@indiamai/integrate_fuse"
+          pip install -I "fuse-element @ git+https://github.com/firedrakeproject/fuse.git@indiamai/compare-topologies"
 
           firedrake-clean
           pip list

Makefile

@@ -4,7 +4,7 @@ all: modules
 .PHONY: modules
 modules:
 	@echo "    Building extension modules"
-	@python setup.py build_ext --inplace > build.log 2>&1 || cat build.log
+	@python3 setup.py build_ext --inplace > build.log 2>&1 || cat build.log
 
 .PHONY: lint
 lint: srclint actionlint dockerlint
@@ -20,18 +20,18 @@ endif
 .PHONY: srclint
 srclint:
 	@echo "    Linting firedrake"
-	@python -m flake8 $(FLAKE8_FORMAT) firedrake
+	@python3 -m flake8 $(FLAKE8_FORMAT) firedrake
 	@echo "    Linting firedrake scripts"
-	@python -m flake8 $(FLAKE8_FORMAT) firedrake/scripts --filename=*
-	@python -m flake8 $(FLAKE8_FORMAT) scripts --filename=*
+	@python3 -m flake8 $(FLAKE8_FORMAT) firedrake/scripts --filename=*
+	@python3 -m flake8 $(FLAKE8_FORMAT) scripts --filename=*
 	@echo "    Linting firedrake tests"
-	@python -m flake8 $(FLAKE8_FORMAT) tests
+	@python3 -m flake8 $(FLAKE8_FORMAT) tests
 	@echo "    Linting PyOP2"
-	@python -m flake8 $(FLAKE8_FORMAT) pyop2
+	@python3 -m flake8 $(FLAKE8_FORMAT) pyop2
 	@echo "    Linting PyOP2 scripts"
-	@python -m flake8 $(FLAKE8_FORMAT) pyop2/scripts --filename=*
+	@python3 -m flake8 $(FLAKE8_FORMAT) pyop2/scripts --filename=*
 	@echo "    Linting TSFC"
-	@python -m flake8 $(FLAKE8_FORMAT) tsfc
+	@python3 -m flake8 $(FLAKE8_FORMAT) tsfc
 
 .PHONY: actionlint
 actionlint:
@@ -58,7 +58,7 @@ dockerlint:
 .PHONY: clean
 clean:
 	@echo "    Cleaning extension modules"
-	@python setup.py clean > /dev/null 2>&1
+	@python3 setup.py clean > /dev/null 2>&1
 	@echo "    RM firedrake/cython/dmplex.*.so"
 	-@rm -f firedrake/cython/dmplex.so > /dev/null 2>&1
 	@echo "    RM firedrake/cython/dmplex.c"
@@ -88,6 +88,7 @@ clean:
 	@echo "    RM tinyasm/*.so"
 	-@rm -f tinyasm/*.so > /dev/null 2>&1
 
+
 # NOTE: It is recommended to run this command from inside the 'firedrake'
 # Docker image to reduce the likelihood of test failures.
 .PHONY: test_durations

firedrake/cython/dmcommon.pyx

@@ -386,6 +386,21 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm,
         #                         0       2
         #                         |       |
         #                         +---1---+
+        # FUSE                    +---1---+
+        #                         |       |
+        #                         0       2
+        #                         |       |
+        #                         +---3---+
+        #if "fuse" in dm.name:
+        # FUSE rules
+        print("FUSE")
+        #plex_cone_new[0] = plex_cone_old[0]
+        #plex_cone_new[1] = plex_cone_old[3]
+        #plex_cone_new[2] = plex_cone_old[2]
+        #plex_cone_new[3] = plex_cone_old[1]
+        #else:
+        # UFC rules
+        #print("UFC")
         plex_cone_new[0] = plex_cone_old[0]
         plex_cone_new[1] = plex_cone_old[2]
         plex_cone_new[2] = plex_cone_old[1]
@@ -476,7 +491,39 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm,
         #                               8   9
         #                                \ /
         #                                14
-        raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
+        # fuse tet
+        fiat_closure[0] = plex_closure[2 * 13]
+        fiat_closure[1] = plex_closure[2 * 11]
+        fiat_closure[2] = plex_closure[2 * 14]
+        fiat_closure[3] = plex_closure[2 * 12]
+        fiat_closure[4] = plex_closure[2 * 7]
+        fiat_closure[5] = plex_closure[2 * 8]
+        fiat_closure[6] = plex_closure[2 * 10]
+        fiat_closure[7] = plex_closure[2 * 6]
+        fiat_closure[8] = plex_closure[2 * 5]
+        fiat_closure[9] = plex_closure[2 * 9]
+        fiat_closure[10] = plex_closure[2 * 4]
+        fiat_closure[11] = plex_closure[2 * 1]
+        fiat_closure[12] = plex_closure[2 * 3]
+        fiat_closure[13] = plex_closure[2 * 2]
+        fiat_closure[14] = plex_closure[2 * 0]
+        # ufc tet fuse
+        #fiat_closure[0] = plex_closure[2 * 14]
+        #fiat_closure[1] = plex_closure[2 * 13]
+        #fiat_closure[2] = plex_closure[2 * 12]
+        #fiat_closure[3] = plex_closure[2 * 11]
+        #fiat_closure[4] = plex_closure[2 * 5]
+        #fiat_closure[5] = plex_closure[2 * 7]
+        #fiat_closure[6] = plex_closure[2 * 6]
+        #fiat_closure[7] = plex_closure[2 * 8]
+        #fiat_closure[8] = plex_closure[2 * 9]
+        #fiat_closure[9] = plex_closure[2 * 10]
+        #fiat_closure[10] = plex_closure[2 * 2]
+        #fiat_closure[11] = plex_closure[2 * 4]
+        #fiat_closure[12] = plex_closure[2 * 1]
+        #fiat_closure[13] = plex_closure[2 * 3]
+        #fiat_closure[14] = plex_closure[2 * 0]
+        # raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
     elif dm.getCellType(p) == PETSc.DM.PolytopeType.QUADRILATERAL:
         # UFCQuadrilateral:       1---7---3
         #                         |       |
@@ -489,7 +536,22 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm,
         #                         1   0   3
         #                         |       |
         #                         6---2---7
-        raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
+        #
+        # FUSE                    1---5---3
+        #                         |       |
+        #                         4   8   6
+        #                         |       |
+        #                         0---7---2
+        fiat_closure[0] = plex_closure[2 * 6]
+        fiat_closure[1] = plex_closure[2 * 5]
+        fiat_closure[2] = plex_closure[2 * 7]
+        fiat_closure[3] = plex_closure[2 * 8]
+        fiat_closure[4] = plex_closure[2 * 1]
+        fiat_closure[5] = plex_closure[2 * 4]
+        fiat_closure[6] = plex_closure[2 * 3]
+        fiat_closure[7] = plex_closure[2 * 2]
+        fiat_closure[8] = plex_closure[2 * 0]
+        # raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
     elif dm.getCellType(p) == PETSc.DM.PolytopeType.HEXAHEDRON:
         # UFCHexahedron:            3--19---7     3--19---7
         #                         13.       |   13  25  15|
@@ -976,11 +1038,15 @@ cdef inline PetscInt _compute_orientation_simplex(PetscInt *fiat_cone,
 
     CHKERR(PetscMalloc1(coneSize, &cone1))
     CHKERR(PetscMalloc1(coneSize, &inds))
+    #print("plex")
     for k in range(coneSize1):
         cone1[k] = plex_cone[k]
+        #print(plex_cone[k])
     n = 0
+    #print("fiat")
     for e in range(coneSize):
         q = fiat_cone[e]
+        #print(q)
         for k in range(coneSize1):
             if q == cone1[k]:
                 inds[n] = k
@@ -1087,21 +1153,25 @@ cdef inline PetscInt _compute_orientation(PETSc.DM dm,
         const PetscInt *cone = NULL
 
     p = cell_closure[cell, e]
+    #print("p", p)
     if dm.getCellType(p) == PETSc.DM.PolytopeType.POINT:
         return 0
     CHKERR(DMPlexGetConeSize(dm.dm, p, &coneSize))
     CHKERR(DMPlexGetCone(dm.dm, p, &cone))
+    #print(coneSize)
     if (entity_cone_map_offset[e + 1] - entity_cone_map_offset[e]) != coneSize:
         raise RuntimeError("FIAT entity cone size != plex point cone size")
     offset = entity_cone_map_offset[e]
     for i in range(coneSize):
         fiat_cone[i] = cell_closure[cell, entity_cone_map[offset + i]]
     if dm.getCellType(p) == PETSc.DM.PolytopeType.SEGMENT or \
-       dm.getCellType(p) == PETSc.DM.PolytopeType.TRIANGLE or \
-       dm.getCellType(p) == PETSc.DM.PolytopeType.TETRAHEDRON:
+       dm.getCellType(p) == PETSc.DM.PolytopeType.TRIANGLE:
         # UFCInterval      <- PETSc.DM.PolytopeType.SEGMENT
         # UFCTriangle      <- PETSc.DM.PolytopeType.TRIANGLE
+       return _compute_orientation_simplex(fiat_cone, cone, coneSize)
+    elif dm.getCellType(p) == PETSc.DM.PolytopeType.TETRAHEDRON:
         # UFCTetrahedron   <- PETSc.DM.PolytopeType.TETRAHEDRON
+        # _reorder_plex_cone(dm, p, cone, plex_cone)
         return _compute_orientation_simplex(fiat_cone, cone, coneSize)
     elif dm.getCellType(p) == PETSc.DM.PolytopeType.QUADRILATERAL:
         # UFCQuadrilateral <- PETSc.DM.PolytopeType.QUADRILATERAL
@@ -1174,7 +1244,7 @@ def entity_orientations(mesh,
         entity_cone_map[i] = entity_cone_list[i]
     for i in range(len(entity_cone_list_offset)):
         entity_cone_map_offset[i] = entity_cone_list_offset[i]
-    #
+
     dm = mesh.topology_dm
     dim = dm.getDimension()
     numCells = cell_closure.shape[0]
@@ -1511,6 +1581,7 @@ def get_cell_nodes(mesh,
                     perm_offset += ceil_ndofs[i] * num_orientations_c[i]
                 else:
                     # FInAT element must eventually add entity_permutations() method
+
                     if extruded_periodic_1_layer:
                         for j in range(ceil_ndofs[i]):
                             cell_nodes[cell, flat_index[k]] = off + j % offset[flat_index[k]]

firedrake/functionspace.py

@@ -6,6 +6,7 @@
 """
 import itertools
 import ufl
+from ufl.cell import as_cell
 import finat.ufl
 
 from pyop2.utils import flatten
@@ -63,7 +64,7 @@ def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant):
                                      degree=degree, variant=variant)
         # If second element was passed in, use it
         lb = finat.ufl.FiniteElement(vfamily,
-                                     cell=ufl.interval,
+                                     cell=as_cell("interval"),
                                      degree=vdegree, variant=variant)
         # Now make the TensorProductElement
         return finat.ufl.TensorProductElement(la, lb)

firedrake/interpolation.py

@@ -1249,7 +1249,9 @@
     #             weights.append(w)
     # assert len(weights) == num_points
     # but for now we just fix the values to what we know works:
-    if element.degree != 0 or not isinstance(element.cell, FIAT.reference_element.Point):
+
+    # not isinstance(element.cell, FIAT.reference_element.Point)
+    if element.degree != 0 or element.cell.get_dimension() != 0:
         raise NotImplementedError("Cross mesh interpolation only implemented for P0DG on vertex cells.")
     num_points = 1
     weights = [1.]*num_points

firedrake/mesh.py

@@ -39,6 +39,7 @@
 from firedrake.adjoint_utils import MeshGeometryMixin
 from pyadjoint import stop_annotating
 import gem
+from fuse import constructCellComplex
 
 try:
     import netgen
@@ -48,6 +49,7 @@
 # Only for docstring
 import mpi4py  # noqa: F401
 from finat.element_factory import as_fiat_cell
+from finat.ufl import as_cell
 
 
 __all__ = [
@@ -317,6 +319,7 @@ def local_facet_orientation_dat(self):
         local_facet_start = offsets[-3]
         local_facet_end = offsets[-2]
         map_from_cell_to_facet_orientations = self.mesh.entity_orientations[:, local_facet_start:local_facet_end]
+
         # Make output data;
         # this is a map from an exterior/interior facet to the corresponding local facet orientation/orientations.
         # Halo data are required by design, but not actually used.
@@ -762,8 +765,7 @@ def make_cell_node_list(self, global_numbering, entity_dofs, entity_permutations
         :arg entity_permutations: FInAT element entity permutations
         :arg offsets: layer offsets for each entity dof (may be None).
         """
-        return dmcommon.get_cell_nodes(self, global_numbering,
-                                       entity_dofs, entity_permutations, offsets)
+        return dmcommon.get_cell_nodes(self, global_numbering, entity_dofs, entity_permutations, offsets)
 
     def make_dofs_per_plex_entity(self, entity_dofs):
         """Returns the number of DoFs per plex entity for each stratum,
@@ -1171,15 +1173,14 @@ 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 as_cell(_cells[tdim][nfacets])
 
     @utils.cached_property
     def _ufl_mesh(self):
@@ -1221,6 +1222,9 @@ def cell_closure(self):
         vertex_numbering = self._vertex_numbering.createGlobalSection(plex.getPointSF())
 
         cell = self.ufl_cell()
+        if hasattr(cell, "to_fiat"):
+            plex.setName('firedrake_default_topology_fuse')
+
         assert tdim == cell.topological_dimension
         if self.submesh_parent is not None and \
                 not (self.submesh_parent.ufl_cell().cellname == "hexahedron" and cell.cellname == "quadrilateral"):
@@ -1242,34 +1246,36 @@ def cell_closure(self):
                 self.submesh_parent.cell_closure,
                 entity_per_cell,
             )
+        elif hasattr(cell, "to_fiat") and cell.cellname == "tetrahedron":
+            plex.setName('firedrake_default_topology_fuse')
+            #   TODO find better way of branching here
+            topology = cell.to_fiat().topology
+            closureSize = sum([len(ents) for _, ents in topology.items()])
+            return dmcommon.create_cell_closure(plex, cell_numbering, closureSize)
         elif cell.is_simplex:
             topology = FIAT.ufc_cell(cell).get_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":
             petsctools.cite("Homolya2016")
             petsctools.cite("McRae2016")
             # Quadrilateral mesh
             cell_ranks = dmcommon.get_cell_remote_ranks(plex)
-
             facet_orientations = dmcommon.quadrilateral_facet_orientations(
                 plex, vertex_numbering, cell_ranks)
-
             cell_orientations = dmcommon.orientations_facet2cell(
                 plex, vertex_numbering, cell_ranks,
                 facet_orientations, cell_numbering)
-
             dmcommon.exchange_cell_orientations(plex,
                                                 cell_numbering,
                                                 cell_orientations)
 
-            return dmcommon.quadrilateral_closure_ordering(
+            res = dmcommon.quadrilateral_closure_ordering(
                 plex, vertex_numbering, cell_numbering, cell_orientations)
+            return res
         elif cell.cellname == "hexahedron":
             # TODO: Should change and use create_cell_closure() for all cell types.
             topology = FIAT.ufc_cell(cell).get_topology()
@@ -1305,7 +1311,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,
@@ -1708,7 +1713,7 @@ def __init__(self, mesh, layers, periodic=False, name=None):
 
     @utils.cached_property
     def _ufl_cell(self):
-        return ufl.TensorProductCell(self._base_mesh.ufl_cell(), ufl.interval)
+        return ufl.TensorProductCell(self._base_mesh.ufl_cell(), as_cell("interval"))
 
     @utils.cached_property
     def _ufl_mesh(self):
@@ -1931,7 +1936,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):
@@ -3256,9 +3261,9 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', peri
     if extrusion_type == 'radial_hedgehog':
         helement = helement.reconstruct(family="DG", variant="equispaced")
     if periodic:
-        velement = finat.ufl.FiniteElement("DP", ufl.interval, 1, variant="equispaced")
+        velement = finat.ufl.FiniteElement("DP", as_cell("interval"), 1, variant="equispaced")
     else:
-        velement = finat.ufl.FiniteElement("Lagrange", ufl.interval, 1)
+        velement = finat.ufl.FiniteElement("Lagrange", as_cell("interval"), 1)
     element = finat.ufl.TensorProductElement(helement, velement)
 
     if gdim is None:

firedrake/tsfc_interface.py

@@ -217,7 +217,6 @@ def compile_form(form, name, parameters=None, split=True, dont_split=(), diagona
         _ = parameters
         parameters = default_parameters["form_compiler"].copy()
         parameters.update(_)
-
     kernels = []
     numbering = form.terminal_numbering()
     all_meshes = extract_domains(form)

pyproject.toml

@@ -44,6 +44,7 @@ dependencies = [
   "scipy",
   "sympy",
   "islpy>=2025.1.5; sys_platform == 'darwin'",
+  "fuse-element @ git+https://github.com/firedrakeproject/fuse.git",
   "vtk",
 ]
 classifiers = [