#include "petsc.h"
/* We need access to the exodus API to format the file */
#include <exodusII.h>

int
main(int argc, char * argv[])
{
    MPI_Comm comm;
    PetscReal lower[2] = { 0, 0 };
    PetscReal upper[2] = { 1, 1 };
    PetscInt faces[2] = { 1, 1 };
    DMBoundaryType periodicity[2] = { DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED };
    DM dm;
    PetscInt nc[1] = { 1 };
    PetscInt n_dofs[3] = { 1, 0, 0 };
    PetscSection s;
    Vec x;
    PetscViewer viewer;
    PetscErrorCode ierr;

    PetscInitialize(&argc, &argv, NULL, NULL);
    comm = PETSC_COMM_WORLD;

    DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, lower, upper, periodicity, PETSC_TRUE, &dm);
    DMSetUp(dm);

    DMSetNumFields(dm, 1);
    DMPlexCreateSection(dm, NULL, nc, n_dofs, 0, NULL, NULL, NULL, NULL, &s);
    DMSetLocalSection(dm, s);
    DMCreateLabel(dm, "Cell Sets");
    DMLabel cs_label;
    DMGetLabel(dm, "Cell Sets", &cs_label);
    DMLabelAddStratum(cs_label, 10);
    PetscInt idxs[] = { 0 };
    IS is;
    ISCreateGeneral(comm, 1, idxs, PETSC_COPY_VALUES, &is);
    DMLabelSetStratumIS(cs_label, 10, is);
    ISDestroy(&is);

    DMCreateGlobalVector(dm, &x);
    VecSetValue(x, 0, 2, INSERT_VALUES);
    VecSetValue(x, 1, 3, INSERT_VALUES);
    VecSetValue(x, 2, 5, INSERT_VALUES);
    VecSetValue(x, 3, 8, INSERT_VALUES);
    ierr = PetscObjectSetName((PetscObject) x, "U");
    CHKERRQ(ierr);

    /* creating the exo viewer and saving the geometry (the DM) */
    PetscViewerCreate(comm, &viewer);
    PetscViewerSetType(viewer, PETSCVIEWEREXODUSII);
    PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
    PetscViewerFileSetName(viewer, "out.exo");
    /* We are using first order nodal elements so order = 1 */
    ierr = PetscViewerExodusIISetOrder(viewer, 1);
    CHKERRQ(ierr);
    DMView(dm, viewer);

    /* "Formatting" the exodus file, i.e. allocating space for the data, time steps etc */
    int exoid = -1;
    char * nodalVarName[1];
    nodalVarName[0] = (char *) "U";
    PetscInt numNodalVar = 1;
    PetscInt step = 0;
    PetscReal time = 0;

    ierr = PetscViewerExodusIIGetId(viewer, &exoid);
    CHKERRQ(ierr);
    PetscStackCallStandard(ex_put_variable_param, (exoid, EX_NODAL, numNodalVar));
    PetscStackCallStandard(ex_put_variable_names, (exoid, EX_NODAL, numNodalVar, nodalVarName));

    /* setting the time step at which we are outputting in the exo file */
    ierr = DMSetOutputSequenceNumber(dm, step, time);
    CHKERRQ(ierr);
    VecView(x, viewer);

    PetscViewerDestroy(&viewer);
    ierr = VecDestroy(&x);
    CHKERRQ(ierr);
    ierr = DMDestroy(&dm);
    CHKERRQ(ierr);

    PetscFinalize();
}