static char help[] = "Nothing";

#include<petsc.h>

int main(int argc, char **argv)
{
   PetscErrorCode ierr;
   DM                 dm = NULL, 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]);
      }
   }
   ierr = PetscFree(offsets);CHKERRQ(ierr);
   ierr = PetscFree(adjacency);CHKERRQ(ierr);
   
   // using DMPlexGetAdjacency
   ierr = PetscPrintf(PETSC_COMM_WORLD, "using \"DMPlexGetAdjacency\"\n"); CHKERRQ(ierr);
   for(i=0; i<nCells; ++i)
   {
      PetscInt adjSize = -1, *adj=NULL;
      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;
}