diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 1a695ce141..4e1dfc4cce 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -341,6 +341,7 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm, comparison with the associated FIAT cones. The reordered DMPlex cones are later used to construct cell_closure and to compute orientations. + """ if dm.getCellType(p) == PETSc.DM.PolytopeType.POINT: raise RuntimeError(f"POINT has no cone") @@ -392,16 +393,19 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm, raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}") -cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm, - PetscInt p, - PetscInt *plex_closure, - PetscInt *fiat_closure): +cdef inline PetscInt _reorder_plex_closure( + PETSc.DM dm, + PetscInt p, + PetscInt *plex_closure, + PetscInt *fiat_closure, + PetscInt *orientations, +): """Reorder DMPlex closure for FIAT closure. :arg dm: The DMPlex object :arg p: The plex point :arg plex_closure: The original DMPlex closure - :arg fiat_closure: The reorderd closure (output) + :arg fiat_closure: The reordered closure (output) This function defines fixed rules to reorder DMPlex closures for FIAT closures. @@ -411,6 +415,8 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm, Indeed the same FIAT closure can be obtained merely using _reorder_plex_cone() and ensuring that the cell orientation is 0. """ + cdef PetscDMPolytopeType ct = dm.getCellType(p) + if dm.getCellType(p) == PETSc.DM.PolytopeType.POINT: raise RuntimeError(f"POINT has no cone") elif dm.getCellType(p) == PETSc.DM.PolytopeType.SEGMENT: @@ -418,20 +424,47 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm, # # PETSc.DM.PolytopeType. 1---0---2 # SEGMENT: - raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}") + fiat_closure[0] = plex_closure[2*0] + orientations[0] = plex_closure[2*0+1] + + fiat_closure[1] = plex_closure[2*1] + orientations[1] = plex_closure[2*1+1] elif dm.getCellType(p) == PETSc.DM.PolytopeType.TRIANGLE: # UFCTriangle: 1 - # | \ + # ↑ ↘ # 5 3 - # | 6 \ - # 0---4---2 + # ↑ 6 ↘ + # 0 → 4 → 2 # - # PETSc.DM.PolytopeType. 4 - # TRIANGLE: | \ - # 3 1 - # | 0 \ - # 6---2---5 - raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}") + # PETSc.DM.PolytopeType. 5 + # TRIANGLE: ↑ ↘ + # 1 2 + # ↑ 0 ↘ + # 4 ← 3 ← 6 + + # TODO: could straightforwardly-ish convert these directly to FIAT orientations + fiat_closure[0] = plex_closure[2*4] + orientations[0] = plex_closure[2*4+1] + + fiat_closure[1] = plex_closure[2*5] + orientations[1] = plex_closure[2*5+1] + + fiat_closure[2] = plex_closure[2*6] + orientations[2] = plex_closure[2*6+1] + + fiat_closure[3] = plex_closure[2*2] + orientations[3] = plex_closure[2*2+1] + + # this edge is flipped + fiat_closure[4] = plex_closure[2*3] + orientations[4] = DMPolytopeTypeComposeOrientation(ct, plex_closure[2*3 + 1], -1) + + fiat_closure[5] = plex_closure[2 * 1] + orientations[5] = plex_closure[2 * 1+1] + + fiat_closure[6] = plex_closure[2 * 0] + orientations[6] = plex_closure[2 * 0+1] + elif dm.getCellType(p) == PETSc.DM.PolytopeType.TETRAHEDRON: # UFCTetrahedron: 0---9---1---9---0 # \ 12 / \ 13 / @@ -523,6 +556,61 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm, raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}") +def cell_orientation_from_vertex_arrangement( + cell_type, + vertex_ordering +): + cdef: + PetscDMPolytopeType ct + PetscInt nvert, narr + PetscInt *arr + + ct = cell_type + nvert = DMPolytopeTypeGetNumVertices(ct) + narr = DMPolytopeTypeGetNumArrangments(ct) + shift = narr // 2 + for o in range(-shift, shift): + arr = DMPolytopeTypeGetVertexArrangment(ct, o) + + # compare to vertex_ordering + equal = True + for i in range(nvert): + if arr[i] != vertex_ordering[i]: + equal = False + break + if equal: + return o + + raise ValueError(f"No valid orientation found for vertex arrangement {vertex_ordering}") + + +def _ordering(iterable): + iterable = list(iterable) + return tuple(iterable.index(p) for p in sorted(iterable)) + + +def orient_mesh(PETSc.DM dm, vertex_numbering): + cdef: + PetscDMPolytopeType ct + PetscInt o, nverts + + pstart, pend = dm.getChart() + for p in range(pstart, pend): + cell_type = dm.getCellType(p) + + # vertices have just a single orientation + if cell_type == PETSc.DM.PolytopeType.POINT: + continue + + ct = cell_type + nverts = DMPolytopeTypeGetNumVertices(ct) + verts = dm.getTransitiveClosure(p)[0][-nverts:] + verts_renum = [vertex_numbering.getOffset(v) for v in verts] + vertex_ordering = _ordering(verts_renum) + o = cell_orientation_from_vertex_arrangement(cell_type, vertex_ordering) + DMPlexOrientPoint(dm.dm, p, o) + + @cython.boundscheck(False) @cython.wraparound(False) def create_cell_closure(PETSc.DM dm, @@ -538,7 +626,7 @@ def create_cell_closure(PETSc.DM dm, PetscInt c, cStart, cEnd, cell, i PetscInt closureSize = _closureSize, closureSize1 PetscInt *closure = NULL - PetscInt *fiat_closure = NULL + PetscInt *fiat_closure = NULL, *fiat_orientations = NULL np.ndarray[PetscInt, ndim=2, mode="c"] cell_closure get_height_stratum(dm.dm, 0, &cStart, &cEnd) @@ -550,16 +638,20 @@ def create_cell_closure(PETSc.DM dm, raise RuntimeError(f"point 0 and point {c} have different cell types") restore_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure) cell_closure = np.empty((cEnd - cStart, closureSize), dtype=IntType) + orientations = np.empty((cEnd - cStart, closureSize), dtype=IntType) CHKERR(PetscMalloc1(closureSize, &fiat_closure)) + CHKERR(PetscMalloc1(closureSize, &fiat_orientations)) for c in range(cStart, cEnd): CHKERR(PetscSectionGetOffset(cell_numbering.sec, c, &cell)) get_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure) - _reorder_plex_closure(dm, c, closure, fiat_closure) + _reorder_plex_closure(dm, c, closure, fiat_closure, fiat_orientations) restore_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure) for i in range(closureSize): cell_closure[cell, i] = fiat_closure[i] + orientations[cell, i] = fiat_orientations[i] CHKERR(PetscFree(fiat_closure)) - return cell_closure + CHKERR(PetscFree(fiat_orientations)) + return cell_closure, orientations @cython.boundscheck(False) diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index b64a5620f5..26184f8028 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -26,6 +26,28 @@ cdef extern from "petsc.h": PETSC_COMPLEX, PETSC_DATATYPE_UNKNOWN +cdef extern from * nogil: + ctypedef enum PetscDMPolytopeType "DMPolytopeType": + DM_POLYTOPE_POINT + DM_POLYTOPE_SEGMENT + DM_POLYTOPE_POINT_PRISM_TENSOR + DM_POLYTOPE_TRIANGLE + DM_POLYTOPE_QUADRILATERAL + DM_POLYTOPE_SEG_PRISM_TENSOR + DM_POLYTOPE_TETRAHEDRON + DM_POLYTOPE_HEXAHEDRON + DM_POLYTOPE_TRI_PRISM + DM_POLYTOPE_TRI_PRISM_TENSOR + DM_POLYTOPE_QUAD_PRISM_TENSOR + DM_POLYTOPE_PYRAMID + DM_POLYTOPE_FV_GHOST + DM_POLYTOPE_INTERIOR_GHOST + DM_POLYTOPE_UNKNOWN + DM_POLYTOPE_UNKNOWN_CELL + DM_POLYTOPE_UNKNOWN_FACE + DM_NUM_POLYTOPES + + cdef extern from "petscsys.h" nogil: int PetscMalloc1(PetscInt,void*) int PetscMalloc2(PetscInt,void*,PetscInt,void*) @@ -54,6 +76,8 @@ cdef extern from "petscdmplex.h" nogil: int DMPlexCreatePointNumbering(PETSc.PetscDM,PETSc.PetscIS*) int DMPlexLabelComplete(PETSc.PetscDM, PETSc.PetscDMLabel) + int DMPlexOrientPoint(PETSc.PetscDM,PetscInt,PetscInt) + cdef extern from "petscdmlabel.h" nogil: struct _n_DMLabel ctypedef _n_DMLabel* DMLabel "DMLabel" @@ -72,6 +96,12 @@ cdef extern from "petscdm.h" nogil: int DMGetLabel(PETSc.PetscDM,char[],DMLabel*) int DMGetPointSF(PETSc.PetscDM,PETSc.PetscSF*) + PetscInt DMPolytopeTypeGetNumVertices(PetscDMPolytopeType) + PetscInt DMPolytopeTypeComposeOrientation(PetscDMPolytopeType,PetscInt,PetscInt) + PetscInt DMPolytopeTypeGetNumArrangments(PetscDMPolytopeType) + PetscInt* DMPolytopeTypeGetArrangment(PetscDMPolytopeType,PetscInt) + PetscInt* DMPolytopeTypeGetVertexArrangment(PetscDMPolytopeType,PetscInt) + cdef extern from "petscdmswarm.h" nogil: int DMSwarmGetLocalSize(PETSc.PetscDM,PetscInt*) int DMSwarmGetCellDM(PETSc.PetscDM, PETSc.PetscDM*) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 9ac82b2e19..da917fa074 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -553,6 +553,7 @@ def callback(self): self.sfXC = sfXB.compose(self.sfBC) if self.sfBC else self.sfXB dmcommon.label_facets(self.topology_dm) dmcommon.complete_facet_labels(self.topology_dm) + # TODO: Allow users to set distribution name if they want to save # conceptually the same mesh but with different distributions, # e.g., those generated by different partitioners. @@ -589,6 +590,20 @@ def callback(self): entity_dofs[-2] = 1 facet_numbering = self.create_section(entity_dofs) self._facet_ordering = dmcommon.get_facet_ordering(self.topology_dm, facet_numbering) + + # Try to give the mesh a globally consistent orientation. This will set + # entity_orientations to all zeros, meaning that no transformations are + # required when packing. + # NOTE: I believe that the quadrilateral numbering algorithm has + # effectively been reproduced in DMPlexOrient which claims to spit + # out an error if the mesh is non-orientable. It appears that all of + # our test meshes *are* orientable at the moment so this is hard to + # verify. + # The difference between having an orientable and a non-orientable mesh + # is that an orientable mesh can be coerced into having a closure ordering + # with a (FIAT) orientation of all zeros. + dmcommon.orient_mesh(self.topology_dm, self._vertex_numbering) + self._callback = callback self.name = name # Set/Generate names to be used when checkpointing. @@ -1089,55 +1104,25 @@ def cell_closure(self): Each row contains ordered cell entities for a cell, one row per cell. """ - plex = self.topology_dm - tdim = plex.getDimension() + return self._closure_and_orientations[0] - # Cell numbering and global vertex numbering - cell_numbering = self._cell_numbering - vertex_numbering = self._vertex_numbering.createGlobalSection(plex.getPointSF()) + @utils.cached_property + def entity_orientations(self): + return self._closure_and_orientations[1] + @utils.cached_property + def _closure_and_orientations(self): cell = self.ufl_cell() - assert tdim == cell.topological_dimension() - if 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": + # can I replace this with calling DMPlexOrient? + if cell.cellname() == "quadrilateral": from firedrake_citations import Citations Citations().register("Homolya2016") Citations().register("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( - plex, vertex_numbering, cell_numbering, cell_orientations) - elif cell.cellname() == "hexahedron": - # TODO: Should change and use create_cell_closure() for all cell types. - topology = FIAT.ufc_cell(cell).get_topology() - closureSize = sum([len(ents) for _, ents in topology.items()]) - return dmcommon.create_cell_closure(plex, cell_numbering, closureSize) - else: - raise NotImplementedError("Cell type '%s' not supported." % cell) - @utils.cached_property - def entity_orientations(self): - return dmcommon.entity_orientations(self, self.cell_closure) + topology = FIAT.ufc_cell(cell).get_topology() + closureSize = sum([len(ents) for _, ents in topology.items()]) + return dmcommon.create_cell_closure(self.topology_dm, self._cell_numbering, closureSize) @PETSc.Log.EventDecorator() def _facets(self, kind):