static char help[] = "";

#include <petscdmplex.h>

const int dof = 1;

#undef __FUNCT__
#define __FUNCT__ "SetupDOFs"
/* Assign dofs (3 for each node) */
PetscErrorCode SetupDOFs(DM dm) {
    PetscSection   s;
    PetscInt       pStart, pEnd, cStart, cEnd, c, fStart, fEnd, f, vStart, vEnd, v, eStart, eEnd, e;
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);CHKERRQ(ierr);
    ierr = PetscSectionSetNumFields(s, 1);CHKERRQ(ierr);
    ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 2, &eStart, &eEnd);CHKERRQ(ierr);
    ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
    ierr = PetscSectionSetChart(s, pStart, pEnd);CHKERRQ(ierr);

    for (v = vStart; v < vEnd; ++v) {
        ierr = PetscSectionSetDof(s, v, dof);CHKERRQ(ierr);
        ierr = PetscSectionSetFieldDof(s, v, 0, dof);CHKERRQ(ierr);
    }
    ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
    ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
    ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "AssignSomeValues"
PetscErrorCode AssignSomeValues(Vec V){
  PetscInt       vectorSize, rStart, rEnd, r;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = VecGetLocalSize(V, &vectorSize);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(V, &rStart, &rEnd);CHKERRQ(ierr);
  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Vector size %d\n",vectorSize);CHKERRQ(ierr);
  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);

  for (r = rStart; r < rEnd; ++r) {
    const PetscScalar value = r;

    ierr = VecSetValues(V, 1, &r, &value, INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(V);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(V);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv) {
    DM             dm, dist;
    Mat            mat;
    Vec            V;
    PetscInt       dim      = 3;
    PetscInt       cells[3] = {1, 1, 2};
    PetscInt       overlap  = 0;
    PetscMPIInt    rank;
    PetscViewer    viewer;
    const char    *fn = "output.vtk";
    PetscErrorCode ierr;

    ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
    ierr = DMPlexCreateHexBoxMesh(PETSC_COMM_WORLD, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &dm); CHKERRQ(ierr);
    ierr = DMPlexDistribute(dm, overlap, NULL, &dist); CHKERRQ(ierr);
    if (dist) {
      ierr = DMDestroy(&dm);CHKERRQ(ierr);
      dm   = dist;
    }
    ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr);
    ierr = SetupDOFs(dm);CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(dm, &V);CHKERRQ(ierr);
    ierr = AssignSomeValues(V);CHKERRQ(ierr);

    ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,fn,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
    ierr = VecView(V,viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

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