[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