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);
    PetscFunctionReturn(0);
}

void AssignSomeValues(Vec V){
    PetscInt vectorSize;
    VecGetSize(V, &vectorSize);
    PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Vector size %d\n",vectorSize);
    PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);

    for (PetscInt i=0;i<vectorSize;i++){
        const PetscScalar value = i;
        VecSetValues(V,1,&i,&value, INSERT_VALUES);
    }
    VecAssemblyBegin(V);
    VecAssemblyEnd(V);
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv) {
    DM             dm, dist;
    Mat            mat;
    PetscInt       dim         = 3;
    PetscInt       cells[]     = {1, 1, 2};
    PetscInt       overlap     = 1;
    PetscErrorCode ierr;
    PetscViewer view;

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

    int rank;
    MPI_Comm_rank (PETSC_COMM_WORLD, &rank);	// get current process id

    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 = SetupDOFs(dm);CHKERRQ(ierr);

    Vec V;
    DMCreateGlobalVector(dm, &V);
    AssignSomeValues(V);

    PetscViewer viewer;
    const char* fn = "output.vtk";
    PetscViewerVTKOpen(PETSC_COMM_WORLD,fn,FILE_MODE_WRITE,&viewer);
    VecView(V,viewer);
    PetscViewerDestroy(&viewer);

    PetscFinalize();

    return 0;
}