change quad orderin based on given cell, add cg3 to fuse tests

Author: indiamai

firedrake/cython/dmcommon.pyx

@@ -368,15 +368,18 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm,
         #                         0       2
         #                         |       |
         #                         +---3---+
-        plex_cone_new[0] = plex_cone_old[0]
-        plex_cone_new[1] = plex_cone_old[2]
-        plex_cone_new[2] = plex_cone_old[1]
-        plex_cone_new[3] = plex_cone_old[3]
-        # print("SWITCHING")
-        #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]
+        if "fuse" in dm.name:
+            # FUSE rules
+            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
+            plex_cone_new[0] = plex_cone_old[0]
+            plex_cone_new[1] = plex_cone_old[2]
+            plex_cone_new[2] = plex_cone_old[1]
+            plex_cone_new[3] = plex_cone_old[3]
     elif dm.getCellType(p) == PETSc.DM.PolytopeType.HEXAHEDRON:
         # UFCHexahedron:            +-------+     +-------+
         #                          /.       |    /   5   /|
@@ -1180,7 +1183,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]
@@ -1517,6 +1520,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/mesh.py

@@ -335,6 +335,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.
@@ -1264,21 +1265,19 @@ def cell_closure(self):
         assert tdim == cell.topological_dimension()
         # plex.viewFromOptions("-dm_view")
         closure, _ = plex.getTransitiveClosure(0)
+        if hasattr(cell, "to_fiat"):
+            plex.setName('firedrake_default_topology_fuse')
+
         if self.submesh_parent is not None:
             return dmcommon.submesh_create_cell_closure_cell_submesh(plex,
                                                                      self.submesh_parent.topology_dm,
                                                                      cell_numbering,
                                                                      self.submesh_parent._cell_numbering,
                                                                      self.submesh_parent.cell_closure)
-        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 hasattr(cell, "to_fiat"):
-        #     #   TODO
+        #     breakpoint()
+        #     plex.setName('firedrake_default_topology_fuse')
+        #     #   TODO find better way of branching here
         #     topology = cell.to_fiat().topology
         #     # entity_per_cell = np.zeros(len(topology), dtype=IntType)
         #     # for d, ents in topology.items():
@@ -1287,7 +1286,15 @@ def cell_closure(self):
         #     #  cell_numbering, entity_per_cell)
         #     topology = FIAT.ufc_cell(cell).get_topology()
         #     closureSize = sum([len(ents) for _, ents in topology.items()])
+        #     breakpoint()
         #     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":
             from firedrake_citations import Citations
             Citations().register("Homolya2016")

tests/firedrake/regression/test_fuse.py

@@ -1,5 +1,6 @@
 from test_helmholtz import helmholtz
 from test_poisson_strong_bcs import run_test
+from fuse.cells import firedrake_triangle
 import pytest
 import numpy as np
 
@@ -15,9 +16,9 @@ def test_poisson_analytic(params, degree, quadrilateral):
 
 @pytest.mark.parametrize(['conv_num', 'degree'],
                          [(p, d)
-                          for p, d in zip([1.8, 2.8], [1, 2])])
-                        #   3.8, 3
-def test_helmholtz(conv_num, degree):
+                          for p, d in zip([1.8, 2.8, 3.8], [1, 2, 3])])
+def test_helmholtz(mocker, conv_num, degree):
+    # mocker.patch('firedrake.mesh.as_cell', return_value=firedrake_triangle().to_ufl("triangle"))
     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:])