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

#include <petscdmplex.h>

PetscErrorCode CreateDiscretization(DM dm)
{
  PetscFE        fe;
  PetscInt       dim;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
  /* Define 1 DOF on cell center of each cell */
  ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, "p_", -1, &fe);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe, "p");CHKERRQ(ierr);
  ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
  ierr = DMCreateDS(dm);CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
  DM             dm, dmDist = NULL;
  Vec            Gvec, Nvec;
  PetscScalar   *n;
  PetscInt       dim = 3, lsize, cum_lsize, ii;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;
  ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);
  ierr = CreateDiscretization(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 = DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_TRUE);CHKERRQ(ierr);
  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
  ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);

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

  /* Add entries in the natural vector */
  ierr = VecGetArray(Nvec, &n); CHKERRQ(ierr);
  for (ii = 0; ii < lsize; ++ii) n[ii] = ii + cum_lsize - lsize;
  ierr = VecRestoreArray(Nvec, &n); CHKERRQ(ierr);

  ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Natural vector:\n\n");CHKERRQ(ierr);
  ierr = VecView(Nvec, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

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

  ierr = DMPlexNaturalToGlobalBegin(dm, Nvec, Gvec);CHKERRQ(ierr);
  ierr = DMPlexNaturalToGlobalEnd(dm, Nvec, Gvec);CHKERRQ(ierr);
  ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\nGlobal vector:\n\n");CHKERRQ(ierr);
  ierr = VecView(Gvec, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  {
    PetscReal   vol, centroid[3];
    PetscInt    cStart, cEnd, c, gref;
    PetscMPIInt rank;

    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr);
    ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\nInformation about the mesh:\n");CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    for (c = cStart; c < cEnd; ++c) {
      ierr = DMPlexGetPointGlobal(dm, c, &gref, NULL); CHKERRQ(ierr);
      ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL); CHKERRQ(ierr);
      ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] cell = %02D; (%f, %f, %f); is_local = %d\n", rank, c, centroid[0], centroid[1], centroid[2], gref>=0);CHKERRQ(ierr);
    }
    ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);
  }

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