[petsc-users] Cell-Cell adjacency in DMPlex
Matthew Knepley
knepley at gmail.com
Tue May 19 08:12:45 CDT 2020
On Mon, May 18, 2020 at 2:49 AM Shashwat Tiwari <shaswat121994 at gmail.com>
wrote:
> 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.
>
You are specifying the size of the array, and its too small. I am attaching
a working version.
Thanks,
Matt
> Regards,
> Shashwat
>
--
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
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200519/f9603714/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: tester.c
Type: application/octet-stream
Size: 3470 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200519/f9603714/attachment-0001.obj>
More information about the petsc-users
mailing list