<div dir="ltr"><div dir="ltr">On Mon, May 18, 2020 at 2:49 AM Shashwat Tiwari <<a href="mailto:shaswat121994@gmail.com">shaswat121994@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div>Hi, <br></div><div>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:</div><div><br></div><div>int main(int argc, char **argv)<br>{<br> PetscErrorCode ierr;<br> DM dm, dmDist = NULL;<br> PetscBool interpolate = PETSC_TRUE, status;<br> PetscBool useCone = PETSC_TRUE, useClosure = PETSC_TRUE;<br> PetscInt overlap = 1;<br> PetscMPIInt rank, size;<br> char filename[PETSC_MAX_PATH_LEN];<br><br> ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);<br> MPI_Comm_rank(PETSC_COMM_WORLD, &rank);<br> MPI_Comm_size(PETSC_COMM_WORLD, &size);<br><br> // get mesh filename<br> ierr = PetscOptionsGetString(NULL, NULL, "-mesh", filename, PETSC_MAX_PATH_LEN, &status);<br> if(status) // gmsh file provided by user<br> {<br> char file[PETSC_MAX_PATH_LEN];<br> ierr = PetscStrcpy(file, filename); CHKERRQ(ierr);<br> ierr = PetscSNPrintf(filename, sizeof filename,"./%s", file); CHKERRQ(ierr);<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading gmsh %s ... ", file); CHKERRQ(ierr);<br> ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, filename, interpolate, &dm); CHKERRQ(ierr);<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);<br> }<br><br> // distribute mesh over processes;<br> ierr = DMSetBasicAdjacency(dm, useCone, useClosure); CHKERRQ(ierr);<br> ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);<br> if(dmDist) <br> {<br> ierr = DMDestroy(&dm); CHKERRQ(ierr);<br> dm = dmDist;<br> }<br> // print mesh information<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "overlap: %d, " <br> "distributed among %d processors\n",<br> overlap, size); CHKERRQ(ierr);<br> <br> // construct ghost cells<br> PetscInt nGhost; // number of ghost cells<br> DM dmG; // DM with ghost cells<br> ierr = DMPlexConstructGhostCells(dm, NULL, &nGhost, &dmG); CHKERRQ(ierr);<br> if(dmG) <br> {<br> ierr = DMDestroy(&dm); CHKERRQ(ierr);<br> dm = dmG;<br> } <br><br> ierr = DMSetUp(dm); CHKERRQ(ierr);<br> ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD); <br><br> // testing cell-cell adjacency<br> ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure); CHKERRQ(ierr);<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "useCone = %d, useClosure = %d\n", useCone, useClosure); CHKERRQ(ierr);<br> <br> // using DMPlexCreateNeighborCSR<br> PetscInt numVertices, *offsets, *adjacency;<br> ierr = DMPlexCreateNeighborCSR(dm, 0, &numVertices, &offsets, &adjacency); CHKERRQ(ierr);<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "using \"DMPlexCreateNeighborCSR\"\n"); CHKERRQ(ierr);<br> ierr = PetscPrintf(PETSC_COMM_WORLD, "numVertices = %d\n", numVertices); CHKERRQ(ierr);<br> PetscInt i, j, nCells = 32;<br> for(i=0; i<nCells; ++i)<br> {<br> PetscPrintf(PETSC_COMM_WORLD, "cell number = %d\n", i);<br> for(j=offsets[i]; j<offsets[i+1]; ++j)<br> {<br> PetscPrintf(PETSC_COMM_WORLD, "nbr cell = %d\n", adjacency[j]);<br> }<br> } <br> <br> // using DMPlexGetAdjacency<br> PetscInt adjSize = 10, *adj=NULL; <br> ierr = PetscPrintf(PETSC_COMM_WORLD, "using \"DMPlexGetAdjacency\"\n"); CHKERRQ(ierr);<br> for(i=0; i<nCells; ++i)<br> {<br> ierr = DMPlexGetAdjacency(dm, i, &adjSize, &adj); CHKERRQ(ierr);<br> ierr = PetscIntView(adjSize, adj, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);<br> ierr = PetscFree(adj); CHKERRQ(ierr);<br> }<br><br> // cleanup<br> ierr = DMDestroy(&dm); CHKERRQ(ierr);<br> ierr = PetscFinalize(); CHKERRQ(ierr);<br> return ierr;<br>}</div><div><br></div><div>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</div></div></div></blockquote><div><br></div><div>Yes, by definition this is using the dual mesh. We could make another function that used the adjacency information, but here</div><div>I was following the input prescription for the mesh partitioners.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div> while "DMPlexGetAdjacency" gives the following error:</div></div></div></blockquote><div><br></div><div>This does look like a bug in the preallocation for this case. I will check it out.</div><div><br></div><div> THanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div>[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>[0]PETSC ERROR: Petsc has generated inconsistent data<br>[0]PETSC ERROR: Invalid mesh exceeded adjacency allocation (10)<br>[0]PETSC ERROR: See <a href="https://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">https://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>[0]PETSC ERROR: Petsc Release Version 3.13.0, unknown <br>[0]PETSC ERROR: ./esue_test on a arch-linux2-c-debug named shashwat by shashwat Mon May 18 12:04:27 2020<br>[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-mpich --download-fblaslapack --with-debugging=1 --download-triangle<br>[0]PETSC ERROR: #1 DMPlexGetAdjacency_Transitive_Internal() line 173 in /home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c<br>[0]PETSC ERROR: #2 DMPlexGetAdjacency_Internal() line 223 in /home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c<br>[0]PETSC ERROR: #3 DMPlexGetAdjacency() line 300 in /home/shashwat/local/petsc/src/dm/impls/plex/plexdistribute.c<br>[0]PETSC ERROR: #4 main() line 583 in esue_test.c<br>[0]PETSC ERROR: PETSc Option Table entries:<br>[0]PETSC ERROR: -mesh test.msh<br>[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<br>application called MPI_Abort(MPI_COMM_SELF, 583077) - process 0</div><div><br></div><div>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.</div><div><br></div><div>Regards,</div><div>Shashwat<br></div></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>