Hi Matt I don’t think the problem is within Petsc - rather somewhere in my code. When I dump the DMPlex using DMView (ascii-info–detail) the ghost mapping seems to be setup correctly.
Is there a better way to determine if a local point is a ghost point? The way I iterate the DMPlex is like this: void iterateDMPlex(DM dm){ Vec coordinates; DMGetCoordinatesLocal(dm, &coordinates); PetscSection defaultSection; DMGetDefaultSection(dm, &defaultSection); PetscSection coordSection; DMGetCoordinateSection(dm, &coordSection); PetscScalar *coords; VecGetArray(coordinates, &coords); DM cdm; DMGetCoordinateDM(dm, &cdm); // iterate (local) mesh PetscInt cellsFrom, cellsTo; std::string s = ""; DMPlexGetHeightStratum(dm, 0, &cellsFrom, &cellsTo); for (PetscInt i=cellsFrom;i<cellsTo;i++) { PetscInt edgesSize; const PetscInt *edgeIndices; DMPlexGetConeSize(dm, i, &edgesSize); DMPlexGetCone(dm, i, &edgeIndices); s = s + "Element: "+std::to_string(i)+"\n"; for (int edgeId = 0;edgeId <edgesSize;edgeId ++){ // ignore edge orientation PetscInt points = edgeIndices[edgeId]; PetscInt edgePoint = edgeIndices[edgeId]; s = s + "\tEdge: "+std::to_string(edgePoint)+"\n"; PetscInt vertexSize; const PetscInt *vertexIndices; DMPlexGetConeSize(dm, edgePoint, &vertexSize); DMPlexGetCone(dm, edgePoint, &vertexIndices); for (int j=0;j<vertexSize;j++){ s = s + "\t\tVertex: "+std::to_string(vertexIndices[j]); s = s + " coords "; PetscScalar* values; VecGetValuesSection(coordinates, coordSection,vertexIndices[j],&values); int dim = 2; for (int k=0;k<dim;k++){ double coordinate = values[k]; s = s +std::to_string(coordinate)+" "; } s = s + (isGhost(cdm, vertexIndices[j])?" ghost":""); s = s + "\n"; } } } VecRestoreArray(coordinates, &coords); int rank; MPI_Comm_rank (PETSC_COMM_WORLD, &rank); // get current process id PetscSynchronizedPrintf(PETSC_COMM_WORLD,”dmplex iteration rank %d \n %s\n",rank, s.c_str()); PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT); } From: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>> Date: Monday 30 November 2015 at 14:08 To: Morten Nobel-Jørgensen <m...@mek.dtu.dk<mailto:m...@mek.dtu.dk>> Cc: "petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>" <petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>> Subject: Re: [petsc-users] DMPlex: Ghost points after DMRefine On Mon, Nov 30, 2015 at 7:01 AM, Morten Nobel-Jørgensen <m...@mek.dtu.dk<mailto:m...@mek.dtu.dk>> wrote: I have a very simple unstructured mesh composed of two triangles (four vertices) with one shared edge using a DMPlex: /|\ / | \ \ | / \|/ After distributing this mesh to two processes, each process owns a triangle. However one process owns tree vertices, while the last vertex is owned by the other process. The problem occurs after uniformly refining the dm. The mesh now looks like this: /|\ /\|/\ \/|\/ \|/ The new center vertex is now not listed as a ghost vertex but instead exists as two individual points. Is there any way that this new center vertex could be created as a ghost vertex during refinement? This could be a bug with the l2g mapping. I do not recreate it when refining, only the SF defining the mapping. Here is an experiment: do not retrieve the mapping until after the refinement. Do you get what you want? If so, I can easily fix this by destroying the map when I refine. Thanks, Matt Kind regards, Morten Ps. Here are some code snippets for getting global point index and test of point is a ghost point: int localToGlobal(DM dm, PetscInt point){ const PetscInt* array; ISLocalToGlobalMapping ltogm; DMGetLocalToGlobalMapping(dm,<ogm); ISLocalToGlobalMappingGetIndices(ltogm, &array); PetscInt res = array[point]; if (res < 0){ // if ghost res = -res +1; } return res; } bool isGhost(DM dm, PetscInt point){ const PetscInt* array; ISLocalToGlobalMapping ltogm; DMGetLocalToGlobalMapping(dm,<ogm); ISLocalToGlobalMappingGetIndices(ltogm, &array); return array[point]<0; } -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener