diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index c0acee0dfe..4a4a8922fd 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -142,9 +142,10 @@ jobs: --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ https://gitlab.com/petsc/petsc.git else - git clone --depth 1 https://gitlab.com/petsc/petsc.git + git clone https://gitlab.com/petsc/petsc.git fi cd petsc + git checkout ksagiyam/fix_submesh_coordinates python3 ../firedrake-repo/scripts/firedrake-configure \ --arch ${{ matrix.arch }} --show-petsc-configure-options | \ xargs -L1 ./configure --with-make-np=8 --download-slepc diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index b67c0c21ec..696a92e805 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -1913,6 +1913,7 @@ def _set_dg_coordinates(PETSc.DM dm, PetscScalar *dg_coords const PetscScalar *firedrake_dg_coords PetscInt n, gdim, cStart, cEnd, c, offset, firedrake_offset, i, j, coord_size, ndof + PETSc.FE coord_fe gdim = firedrake_dg_coord_vec.getBlockSize() coord_dm = dm.getCoordinateDM() @@ -1945,6 +1946,7 @@ def _set_dg_coordinates(PETSc.DM dm, dg_coord_vec, True) dm.setCellCoordinatesLocal(dg_coord_vec) + dm.createCoordinateSpace(1, True, False) @cython.boundscheck(False) @@ -3691,7 +3693,8 @@ def create_halo_exchange_sf(PETSc.DM dm): @cython.wraparound(False) def submesh_create(PETSc.DM dm, label_name, - PetscInt label_value): + PetscInt label_value, + PetscBool ignore_label_halo): """Create submesh. Parameters @@ -3702,6 +3705,8 @@ def submesh_create(PETSc.DM dm, Name of the label label_value : int Value in the label + ignore_label_halo : bool + If labeled points in the halo are ignored. """ cdef: @@ -3710,7 +3715,7 @@ def submesh_create(PETSc.DM dm, PETSc.SF ownership_transfer_sf = PETSc.SF() label = dm.getLabel(label_name) - CHKERR(DMPlexFilter(dm.dm, label.dmlabel, label_value, PETSC_FALSE, PETSC_TRUE, &ownership_transfer_sf.sf, &subdm.dm)) + CHKERR(DMPlexFilter(dm.dm, label.dmlabel, label_value, ignore_label_halo, PETSC_TRUE, &ownership_transfer_sf.sf, &subdm.dm)) submesh_update_facet_labels(dm, subdm) submesh_correct_entity_classes(dm, subdm, ownership_transfer_sf) return subdm @@ -3923,3 +3928,36 @@ def submesh_create_cell_closure_cell_submesh(PETSc.DM subdm, CHKERR(PetscFree(subpoint_indices_inv)) CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) return subcell_closure + + +@cython.boundscheck(False) +@cython.wraparound(False) +def get_dm_cell_types(PETSc.DM dm): + """Return all cell types in the mesh. + + Parameters + ---------- + dm : PETSc.DM + The parent dm. + + Returns + ------- + tuple + Tuple of all cell types in the mesh. + + """ + cdef: + PetscInt cStart, cEnd, c + np.ndarray found, found_all + PetscDMPolytopeType celltype + + cStart, cEnd = dm.getHeightStratum(0) + found = np.zeros((DM_NUM_POLYTOPES, ), dtype=IntType) + found_all = np.zeros((DM_NUM_POLYTOPES, ), dtype=IntType) + for c in range(cStart, cEnd): + CHKERR(DMPlexGetCellType(dm.dm, c, &celltype)) + found[celltype] = 1 + dm.comm.tompi4py().Allreduce(found, found_all, op=MPI.MAX) + return tuple( + polytope_type_enum for polytope_type_enum, found in enumerate(found_all) if found + ) diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index 750a3c3714..a94de11e56 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -29,6 +29,27 @@ cdef extern from "petscsys.h" nogil: int PetscFree2(void*,void*) int PetscSortIntWithArray(PetscInt,PetscInt[],PetscInt[]) +cdef extern from "petscdmtypes.h" 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 "petscdmplex.h" nogil: int DMPlexGetHeightStratum(PETSc.PetscDM,PetscInt,PetscInt*,PetscInt*) int DMPlexGetDepthStratum(PETSc.PetscDM,PetscInt,PetscInt*,PetscInt*) @@ -56,6 +77,9 @@ cdef extern from "petscdmplex.h" nogil: int DMPlexGetSubpointMap(PETSc.PetscDM,PETSc.PetscDMLabel*) int DMPlexSetSubpointMap(PETSc.PetscDM,PETSc.PetscDMLabel) + int DMPlexSetCellType(PETSc.PetscDM,PetscInt,PetscDMPolytopeType) + int DMPlexGetCellType(PETSc.PetscDM,PetscInt,PetscDMPolytopeType*) + cdef extern from "petscdmlabel.h" nogil: struct _n_DMLabel ctypedef _n_DMLabel* DMLabel "DMLabel" diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 0be58d6d95..28ff0af726 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -742,6 +742,12 @@ def ufl_mesh(self): """ return self._ufl_mesh + @property + @abc.abstractmethod + def dm_cell_types(self): + """All ``DM.PolytopeType``s of cells in the mesh.""" + pass + @property @abc.abstractmethod def cell_closure(self): @@ -1247,6 +1253,11 @@ def _renumber_entities(self, reorder): reordering = None return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, reordering) + @property + def dm_cell_types(self): + """All ``DM.PolytopeType``s of cells in the mesh.""" + return dmcommon.get_dm_cell_types(self.topology_dm) + @utils.cached_property def cell_closure(self): """2D array of ordered cell closures @@ -1744,6 +1755,11 @@ def _ufl_mesh(self): cell = self._ufl_cell return ufl.Mesh(finat.ufl.VectorElement("Lagrange", cell, 1, dim=cell.topological_dimension())) + @property + def dm_cell_types(self): + """All ``DM.PolytopeType``s of cells in the mesh.""" + raise NotImplementedError("Notimplemented for ExtrudedMeshTopology") + @utils.cached_property def cell_closure(self): """2D array of ordered cell closures @@ -1982,6 +1998,11 @@ def _renumber_entities(self, reorder): else: return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, None) + @property + def dm_cell_types(self): + """All ``DM.PolytopeType``s of cells in the mesh.""" + return (PETSc.DM.PolytopeType.POINT,) + @utils.cached_property # TODO: Recalculate if mesh moves def cell_closure(self): """2D array of ordered cell closures @@ -4596,7 +4617,7 @@ def SubDomainData(geometric_expr): return op2.Subset(m.cell_set, indices) -def Submesh(mesh, subdim, subdomain_id, label_name=None, name=None): +def Submesh(mesh, subdim, subdomain_id, label_name=None, ignore_label_halo=False, name=None): """Construct a submesh from a given mesh. Parameters @@ -4609,6 +4630,8 @@ def Submesh(mesh, subdim, subdomain_id, label_name=None, name=None): Subdomain ID representing the submesh. label_name : str Name of the label to search ``subdomain_id`` in. + ignore_label_halo : bool + If labeled points in the halo are ignored. name : str Name of the submesh. @@ -4662,7 +4685,7 @@ def Submesh(mesh, subdim, subdomain_id, label_name=None, name=None): if label_name is None: label_name = dmcommon.CELL_SETS_LABEL name = name or _generate_default_submesh_name(mesh.name) - subplex = dmcommon.submesh_create(plex, label_name, subdomain_id) + subplex = dmcommon.submesh_create(plex, label_name, subdomain_id, ignore_label_halo) subplex.setName(_generate_default_mesh_topology_name(name)) if subplex.getDimension() != subdim: raise RuntimeError(f"Found subplex dim ({subplex.getDimension()}) != expected ({subdim})") diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index 446e042b35..9900a739e9 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -131,36 +131,32 @@ def MeshHierarchy(mesh, refinement_levels, raise RuntimeError("Cannot create a NetgenHierarchy from a mesh that has not been generated by\ Netgen.") + mesh_with_overlap = mesh + mesh.init() + dm_cell_type, = mesh.dm_cell_types + tdim = mesh.topology_dm.getDimension() + # Virtually "invert" addOverlap. + # -- This is algorithmically guaranteed. + mesh = firedrake.Submesh(mesh_with_overlap, tdim, dm_cell_type, label_name="celltype", ignore_label_halo=True) cdm = mesh.topology_dm + cdm.removeLabel("pyop2_core") + cdm.removeLabel("pyop2_owned") + cdm.removeLabel("pyop2_ghost") cdm.setRefinementUniform(True) dms = [] if mesh.comm.size > 1 and mesh._grown_halos: raise RuntimeError("Cannot refine parallel overlapped meshes " "(make sure the MeshHierarchy is built immediately after the Mesh)") - parameters = {} - if distribution_parameters is not None: - parameters.update(distribution_parameters) - else: - parameters.update(mesh._distribution_parameters) - - parameters["partition"] = False - distribution_parameters = parameters - if callbacks is not None: before, after = callbacks else: before = after = lambda dm, i: None - for i in range(refinement_levels*refinements_per_level): if i % refinements_per_level == 0: before(cdm, i) rdm = cdm.refine() if i % refinements_per_level == 0: after(rdm, i) - rdm.removeLabel("pyop2_core") - rdm.removeLabel("pyop2_owned") - rdm.removeLabel("pyop2_ghost") - dms.append(rdm) cdm = rdm # Fix up coords if refining embedded circle or sphere @@ -172,20 +168,38 @@ def MeshHierarchy(mesh, refinement_levels, coords = cdm.getCoordinatesLocal().array.reshape(-1, mesh.geometric_dimension()) scale = mesh._radius / np.linalg.norm(coords, axis=1).reshape(-1, 1) coords *= scale - - meshes = [mesh] + [mesh_builder(dm, dim=mesh.geometric_dimension(), - distribution_parameters=distribution_parameters, - reorder=reorder, comm=mesh.comm) - for dm in dms] - - lgmaps = [] - for i, m in enumerate(meshes): - no = impl.create_lgmap(m.topology_dm) + lgmaps_without_overlap = [ + impl.create_lgmap(dm) for dm in [mesh.topology_dm] + dms + ] + parameters = {} + if distribution_parameters is not None: + parameters.update(distribution_parameters) + else: + parameters.update(mesh_with_overlap._distribution_parameters) + parameters["partition"] = False + meshes = [mesh_with_overlap] + [ + mesh_builder( + dm, + dim=mesh_with_overlap.geometric_dimension(), + distribution_parameters=parameters, + reorder=reorder, + comm=mesh.comm, + ) + for dm in dms + ] + for m in meshes: m.init() - o = impl.create_lgmap(m.topology_dm) + #distribution_parameters_noop={ + # "partition": False, + # "overlap_type": (firedrake.mesh.DistributedMeshOverlapType.NONE, 0), + #} + lgmaps_with_overlap = [] + for i, m in enumerate(meshes): + lgmaps_with_overlap.append(impl.create_lgmap(m.topology_dm)) m.topology_dm.setRefineLevel(i) - lgmaps.append((no, o)) - + lgmaps = [ + (no, o) for no, o in zip(lgmaps_without_overlap, lgmaps_with_overlap) + ] coarse_to_fine_cells = [] fine_to_coarse_cells = [None] for (coarse, fine), (clgmaps, flgmaps) in zip(zip(meshes[:-1], meshes[1:]), diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index 99d687449e..2590682c8d 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -734,6 +734,7 @@ def TensorRectangleMesh( plex = mesh.plex_from_cell_list( 2, cells, coords, comm, mesh._generate_default_mesh_topology_name(name) ) + plex.createCoordinateSpace(1, False, False) # mark boundary facets plex.createLabel(dmcommon.FACE_SETS_LABEL) diff --git a/tests/firedrake/multigrid/test_basics.py b/tests/firedrake/multigrid/test_basics.py index ec80c1ee4f..d4161ebca1 100644 --- a/tests/firedrake/multigrid/test_basics.py +++ b/tests/firedrake/multigrid/test_basics.py @@ -60,13 +60,3 @@ def test_refine_square_ncell_parallel(): # Should be fewer than 4 times the number of coarse cells due to # halo shrinking. assert mh[1].num_cells() < 4 * mh[0].num_cells() - - -@pytest.mark.parallel(nprocs=2) -def test_refining_overlapped_mesh_fails_parallel(): - m = UnitSquareMesh(4, 4) - - m.init() - - with pytest.raises(RuntimeError): - MeshHierarchy(m, 1) diff --git a/tests/firedrake/multigrid/test_grid_transfer.py b/tests/firedrake/multigrid/test_grid_transfer.py index 40b21bb3b2..0168e710b1 100644 --- a/tests/firedrake/multigrid/test_grid_transfer.py +++ b/tests/firedrake/multigrid/test_grid_transfer.py @@ -269,7 +269,8 @@ def test_grid_transfer_deformed(deformed_hierarchy, deformed_transfer_type): run_prolongation(deformed_hierarchy, vector, space, degrees) -@pytest.fixture(params=["interval", "triangle", "quadrilateral", "tetrahedron"], scope="module") +#@pytest.fixture(params=["interval", "triangle", "quadrilateral", "tetrahedron"], scope="module") +@pytest.fixture(params=["interval"], scope="module") def periodic_cell(request): return request.param @@ -313,7 +314,6 @@ def exact_primal_periodic(mesh, vector, degree): return expr -@pytest.mark.parallel(nprocs=3) def test_grid_transfer_periodic(periodic_hierarchy, periodic_space): degrees = [4] vector = False