[petsc-users] Cell-Cell adjacency in DMPlex

[Shashwat Tiwari]

Hi,
I am writing a second order finite volume scheme on 2D unstructured grids
using DMPlex. For this, I need a vertex neighbor stencil for the cells
(i.e. stencil of the cells which share atleast one vertex with the given
cell). I am trying to achieve this using two of the Petsc functions,
"DMPlexCreateNeighborCSR" and "DMPlexGetAdjacency" and using
"DMSetBasicAdjacency"  function to set whether to use closure or not in
order to get face neighbors and vertex neighbors. Here is a simple code I
have written to test these functions out:

int main(int argc, char **argv)
{
   PetscErrorCode ierr;
   DM                 dm, dmDist = NULL;
   PetscBool      interpolate = PETSC_TRUE, status;
   PetscBool      useCone = PETSC_TRUE, useClosure = PETSC_TRUE;
   PetscInt         overlap = 1;
   PetscMPIInt   rank, size;
   char               filename[PETSC_MAX_PATH_LEN];

   ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
   MPI_Comm_size(PETSC_COMM_WORLD, &size);

   // get mesh filename
   ierr = PetscOptionsGetString(NULL, NULL, "-mesh", filename,
PETSC_MAX_PATH_LEN, &status);
   if(status) // gmsh file provided by user
   {
      char file[PETSC_MAX_PATH_LEN];
      ierr = PetscStrcpy(file, filename); CHKERRQ(ierr);
      ierr = PetscSNPrintf(filename, sizeof filename,"./%s", file);
CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading gmsh %s ... ", file);
CHKERRQ(ierr);
      ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, filename,
interpolate, &dm); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);
   }

   // distribute mesh over processes;
   ierr = DMSetBasicAdjacency(dm, useCone, useClosure); CHKERRQ(ierr);
   ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);
   if(dmDist)
   {
      ierr = DMDestroy(&dm); CHKERRQ(ierr);
      dm = dmDist;
   }
   // print mesh information
   ierr = PetscPrintf(PETSC_COMM_WORLD, "overlap: %d, "
                                        "distributed among %d processors\n",
                                         overlap, size); CHKERRQ(ierr);

   // construct ghost cells
   PetscInt nGhost; // number of ghost cells
   DM dmG; // DM with ghost cells
   ierr = DMPlexConstructGhostCells(dm, NULL, &nGhost, &dmG); CHKERRQ(ierr);
   if(dmG)
   {
      ierr = DMDestroy(&dm); CHKERRQ(ierr);
      dm = dmG;
   }

   ierr = DMSetUp(dm); CHKERRQ(ierr);
   ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);

   // testing cell-cell adjacency
   ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "useCone = %d, useClosure = %d\n",
useCone, useClosure); CHKERRQ(ierr);

   // using DMPlexCreateNeighborCSR
   PetscInt numVertices, *offsets, *adjacency;
   ierr = DMPlexCreateNeighborCSR(dm, 0, &numVertices, &offsets,
&adjacency); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "using
\"DMPlexCreateNeighborCSR\"\n"); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "numVertices = %d\n", numVertices);
CHKERRQ(ierr);
   PetscInt i, j, nCells = 32;
   for(i=0; i<nCells; ++i)
   {
      PetscPrintf(PETSC_COMM_WORLD, "cell number = %d\n", i);
      for(j=offsets[i]; j<offsets[i+1]; ++j)
      {
         PetscPrintf(PETSC_COMM_WORLD, "nbr cell = %d\n", adjacency[j]);
      }
   }

   // using DMPlexGetAdjacency
   PetscInt adjSize = 10, *adj=NULL;
   ierr = PetscPrintf(PETSC_COMM_WORLD, "using \"DMPlexGetAdjacency\"\n");
CHKERRQ(ierr);
   for(i=0; i<nCells; ++i)
   {
      ierr = DMPlexGetAdjacency(dm, i, &adjSize, &adj); CHKERRQ(ierr);
      ierr = PetscIntView(adjSize, adj, PETSC_VIEWER_STDOUT_WORLD);
CHKERRQ(ierr);
      ierr = PetscFree(adj); CHKERRQ(ierr);
   }

   // cleanup
   ierr = DMDestroy(&dm); CHKERRQ(ierr);
   ierr = PetscFinalize(); CHKERRQ(ierr);
   return ierr;
}

I am reading the mesh from a gmsh file which I have attached. When I set
the "useClosure" to PETSC_FALSE, both the functions give the expected
result, i.e. the list of all the neighbors sharing a face with the given
cell. But, when I set it to PETSC_TRUE, the "DMPlexCreateNeighborCSR" still
gives the face neighbors while "DMPlexGetAdjacency" gives the following
error:

[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Invalid mesh exceeded adjacency allocation (10)
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.13.0, unknown
[0]PETSC ERROR: ./esue_test on a arch-linux2-c-debug named shashwat by
shashwat Mon May 18 12:04:27 2020
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
--with-fc=gfortran --download-mpich --download-fblaslapack
--with-debugging=1 --download-triangle
[0]PETSC ERROR: #1 DMPlexGetAdjacency_Transitive_Internal() line 173 in
/home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c
[0]PETSC ERROR: #2 DMPlexGetAdjacency_Internal() line 223 in
/home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c
[0]PETSC ERROR: #3 DMPlexGetAdjacency() line 300 in
/home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c
[0]PETSC ERROR: #4 main() line 583 in esue_test.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -mesh test.msh
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-maint at mcs.anl.gov----------
application called MPI_Abort(MPI_COMM_SELF, 583077) - process 0

Kindly look into it and let me know what mistake I might be making, and
also, if there is some other way in Petsc to get the Vertex Neighbors for a
given cell.

Regards,
Shashwat
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200518/968e3283/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.geo
Type: application/octet-stream
Size: 374 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200518/968e3283/attachment.obj>