static char help[] = "Define a simple field over the mesh\n\n";

#include <petscdmplex.h>

int main(int argc, char **argv)
{
  DM             dm, dmDist = NULL;
  PetscInt N = 2, dim = 3;
  PetscInt size, myrank;
  PetscInt c, cStart, cEnd;
  Vec            coordinates;
  PetscSection   coordSection;
  PetscScalar   *coords;
  PetscErrorCode ierr;

  const PetscInt  faces[3] = {N,N,N  };
  const PetscReal lower[3] = {0.0,0.0,0.0};
  const PetscReal upper[3] = {1.0,1.0,1.0};

  ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;


  ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dim,PETSC_FALSE,faces,lower,upper,
                             NULL,PETSC_TRUE,&dm); CHKERRQ(ierr);

  ierr = DMSetUseNatural(dm,PETSC_TRUE); CHKERRQ(ierr);
  ierr = DMPlexDistribute(dm, 1, NULL, &dmDist);CHKERRQ(ierr);
  if (dmDist) {ierr = DMDestroy(&dm);CHKERRQ(ierr); dm = dmDist;}
  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);

  PetscSection sec;
  PetscInt pStart, pEnd;

  // Define 1 DOF on cell center of the mesh
  ierr = PetscSectionCreate(PETSC_COMM_WORLD, &sec);
  ierr = PetscSectionSetNumFields(sec,1);
  ierr = PetscSectionSetFieldName(sec,0,"p");
  ierr = PetscSectionSetFieldComponents(sec,0,1);
  ierr = DMPlexGetChart(dm,&pStart,&pEnd);
  ierr = PetscSectionSetChart(sec,pStart,pEnd);
  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
  for (PetscInt c=cStart; c<cEnd; c++) {
    ierr = PetscSectionSetFieldDof(sec,c,0,1);
    ierr = PetscSectionSetDof(sec,c,1);
  }
  ierr = PetscSectionSetUp(sec); CHKERRQ(ierr);
  ierr = DMSetSection(dm,sec); CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
  ierr = DMSetBasicAdjacency(dm,PETSC_TRUE,PETSC_TRUE); CHKERRQ(ierr);

  MPI_Comm_size(PETSC_COMM_WORLD, &size);
  MPI_Comm_rank(PETSC_COMM_WORLD, &myrank);

  Vec Gvec,Nvec;
  PetscInt lsize, cum_lsize;
  PetscScalar *n, *g;

  // Create the natural vector
  ierr = DMCreateGlobalVector(dm,&Nvec);
  ierr = VecGetLocalSize(Nvec,&lsize);
  MPI_Scan(&lsize,&cum_lsize,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);

  // Add entries in the natural vector
  ierr = VecGetArray(Nvec, &n); CHKERRQ(ierr);
  for (PetscInt ii=0; ii<lsize; ii++) n[ii] = ii + cum_lsize - lsize;
  ierr = VecRestoreArray(Nvec, &n); CHKERRQ(ierr);
  
  if (myrank==0) printf("Natural vector:\n\n");
  ierr = VecView(Nvec,PETSC_VIEWER_STDOUT_WORLD);

  // Create a global vector
  ierr = DMCreateGlobalVector(dm,&Gvec);
  ierr = VecGetLocalSize(Gvec,&lsize);

  // Do Nat-to-Global
  ierr = DMPlexNaturalToGlobalBegin(dm,Nvec,Gvec);
  ierr = DMPlexNaturalToGlobalEnd(  dm,Nvec,Gvec);
  if (myrank==0) printf("\nGlobal vector:\n\n");
  ierr = VecView(Gvec,PETSC_VIEWER_STDOUT_WORLD);

  ierr = DMGetCoordinateSection(dm, &coordSection); CHKERRQ(ierr);
  ierr = DMGetCoordinatesLocal (dm, &coordinates); CHKERRQ(ierr);
  ierr = VecGetArray(coordinates,&coords); CHKERRQ(ierr);
  ierr = VecGetArray(Gvec, &g); CHKERRQ(ierr);

  PetscReal V, X[3], Normal[3];
  PetscInt gref, junkInt;

  if (myrank==0) printf("\nInformation about the mesh:\n");
  for(PetscInt rank=0;rank<size;rank++) {
    if (rank==myrank) {
      printf("\nRank = %d\n",myrank);
      for (c=cStart;c<cEnd;c++) {
        ierr = DMPlexGetPointGlobal(dm,c,&gref,&junkInt); CHKERRQ(ierr);
        ierr = DMPlexComputeCellGeometryFVM(dm,c,&V,X,Normal); CHKERRQ(ierr);
        printf("local_id = %02d; (%f, %f, %f); is_local = %d\n", c, X[0],X[1],X[2],(gref>=0));
      }
    }
    MPI_Barrier(PETSC_COMM_WORLD);
  }

  ierr = VecRestoreArray(coordinates,&coords); CHKERRQ(ierr);
  ierr = VecRestoreArray(Gvec, &g); CHKERRQ(ierr);


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