static char help[] = "Second Order TVD Finite Volume Example.\n";

#include <petscdmplex.h>
#include <petscds.h>
#include <petscts.h>
#include <petscsf.h> /* For SplitFaces() */

#define DIM 2                   /* Geometric dimension */
#define ALEN(a) (sizeof(a)/sizeof((a)[0]))

static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer);
PetscErrorCode SetInitialCondition(DM dm, Vec X);
PetscErrorCode CalcConVarGradient(DM dm, Vec X, const PetscInt gradStart);
PetscErrorCode CalcLimiter(DM dm);

    
PetscScalar Dot(const PetscScalar *x,const PetscScalar *y, const PetscInt numComp);
PetscReal Norm(const PetscScalar *x, const PetscInt numComp); 




#undef __FUNCT__
#define __FUNCT__ "Dot"
PetscReal Dot(const PetscScalar *x,const PetscScalar *y, const PetscInt numComp) 
{ 
    int i=0;
    PetscReal temp = 0;
    for(; i < numComp; i++)
    {
        temp += x[i] * y[i];
    }
    return temp;
}

#undef __FUNCT__
#define __FUNCT__ "Norm"
PetscReal Norm(const PetscScalar *x, const PetscInt numComp) 
{ 
    return PetscSqrtReal(PetscAbsScalar(Dot(x,x,numComp)));
}


#undef __FUNCT__
#define __FUNCT__ "OutputVTK"
static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);CHKERRQ(ierr);

  ierr = PetscViewerSetType(*viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
  ierr = PetscViewerPushFormat(*viewer, PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
  /*ierr = PetscViewerPushFormat(*viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);*/
  ierr = PetscViewerFileSetName(*viewer, filename);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SetInitialCondition"
PetscErrorCode SetInitialCondition(DM dm, Vec X)
{
    DM                 dmCell;
    Vec                cellgeom;
    const PetscScalar *cgeom;
    PetscScalar       *x;
    PetscInt           cStart, cEnd, cEndInterior, c;
    PetscErrorCode     ierr;

    PetscFunctionBeginUser;
    ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);

    ierr = DMPlexTSGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecGetArray(X, &x);CHKERRQ(ierr);
    for (c = cStart; c < cEndInterior; ++c) {
        PetscScalar    *xc;

        ierr = DMPlexPointLocalRef(dm,c,x,&xc);CHKERRQ(ierr);

        /* conVar */
        xc[0] = 1.0;
        xc[1] = 2.0;
        xc[2] = 3.0;
        xc[3] = 4.0;
        for(int i = 4; i < 12; i++) xc[i] = 1.0; /* gradient */
        for(int i = 12; i < 16; i++) xc[i] = 1.0; /* limiter */
        for(int i = 16; i < 20; i++) xc[i] = 0.0; /* flux */
    }
    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArray(X, &x);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcConVarGradient"
PetscErrorCode CalcConVarGradient(DM dm, Vec X, const PetscInt gradStart)
{
    DM                 dmCell;
    Vec                cellgeom;
    const PetscScalar *cgeom;
    PetscScalar       *x;
    PetscInt           cStart, cEnd, cEndInterior, c;
    PetscErrorCode     ierr;

    PetscInt i, j;
    PetscFunctionBeginUser;
    ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);

    ierr = DMPlexTSGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecGetArray(X, &x);CHKERRQ(ierr);
    for (c = cStart; c < cEndInterior; ++c) 
    {
        PetscInt adjSize, adjSizeInterior, *adj=NULL;
        PetscScalar *xc;
        const PetscScalar *xcAdj;
        const PetscFVCellGeom *cg, *cgAdj;
        PetscScalar *arrayConVar, *arrayDeltaX, *arrayDeltaY;

        ierr = DMPlexGetAdjacency(dm, c, &adjSize, &adj); CHKERRQ(ierr);

        /*PetscSynchronizedPrintf(PETSC_COMM_WORLD, "adjSize = %d\n", adjSize);*/
        adjSizeInterior = -1;
        for(i = 0; i < adjSize; i++)
        {
            if(adj[i] < cEndInterior)
            {
                adjSizeInterior++;
                /*PetscSynchronizedPrintf(PETSC_COMM_WORLD, "adj[%d]= %d\n", i, adj[i]);*/
            }
        }
        /*
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "adjSizeInterior = %d\n", adjSizeInterior);
        PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
        */

        ierr = PetscMalloc3(adjSizeInterior*4, &arrayConVar, adjSizeInterior, &arrayDeltaX, adjSizeInterior, &arrayDeltaY); CHKERRQ(ierr);





        PetscFree3(arrayConVar, arrayDeltaX, arrayDeltaY); CHKERRQ(ierr);
        ierr = PetscFree(adj); CHKERRQ(ierr);



    }
    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArray(X, &x);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcLimiter"
PetscErrorCode CalcLimiter(DM dm)
{
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
    MPI_Comm comm;
    PetscErrorCode ierr;
    PetscInt i, j;

    char filename[PETSC_MAX_PATH_LEN] = "./cavity.exo";

    DM dm;
    PetscInt overlap = 2, dim = 2, numFields, numBC, numBCRegions, *bcField;
    PetscInt gradStart, limStart;
    IS  *bcPointIS, bcIS, innerIS;
    const PetscInt *bcId;
    Vec X, localX;

    PetscSection section;

    PetscInt timeStep = 0, maxTimeStep = 1;
    PetscViewer viewer;
    PetscBool hasLabel;

    ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
    comm = PETSC_COMM_WORLD;

    /* Reading the grid. */
    ierr = DMPlexCreateFromFile(comm, filename, PETSC_TRUE, &dm);CHKERRQ(ierr);

    /* Creating lables for boundaries. */
    ierr = DMCreateLabel(dm, "Face Sets");CHKERRQ(ierr);

    /* Distribute mesh over processes */
    {
        DM dmDist;
        ierr = DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);CHKERRQ(ierr);
        ierr = DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE);CHKERRQ(ierr);
        ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist);CHKERRQ(ierr);
        if (dmDist) 
        {
            ierr = DMDestroy(&dm);CHKERRQ(ierr);
            dm   = dmDist;
        }
    }

    {
        DM gdm;
        PetscInt numGhostCells;

        ierr = DMPlexConstructGhostCells(dm, "ghost", &numGhostCells, &gdm);CHKERRQ(ierr);
        ierr = DMDestroy(&dm);CHKERRQ(ierr);
        dm   = gdm;

        /*
        PetscSynchronizedPrintf(comm, "numGhostCells = %d\n", numGhostCells);
        PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
        */
    }

    /*
     * Create section.
     */
    {
        numFields  = 4;
        PetscInt numComp[4];
        PetscInt numDof[12];

        /* Create scalar fields rho and rhoE, a vector field rhoU */
        numComp[0] = 4;
        numComp[1] = 8;
        numComp[2] = 4;
        numComp[3] = 4;

        gradStart = 4; limStart = 12;
        for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0;

        /* Defining the fields on cell */
        numDof[0*(dim+1)+dim] = numComp[0]; 
        numDof[1*(dim+1)+dim] = numComp[1]; 
        numDof[2*(dim+1)+dim] = numComp[2];     
        numDof[3*(dim+1)+dim] = numComp[3];     

        /* Checking whether the boundaries are labeled. */
        ierr = DMHasLabel(dm, "Face Sets", &hasLabel);CHKERRQ(ierr);
        if(!hasLabel) return 0;

        /* Checking the number of boundaries.*/
        ierr = DMGetLabelSize(dm, "Face Sets", &numBC);CHKERRQ(ierr);
        /*
        PetscSynchronizedPrintf(comm, "numBC = %d\n", numBC);
        PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
        */

        ierr = DMGetLabelIdIS(dm, "Face Sets", &bcIS);CHKERRQ(ierr);
        ierr = ISGetIndices(bcIS, &bcId);CHKERRQ(ierr);

        for(i = 0; i < numBC; i++)
        {
            ierr = ISGetSize(bcIS, &numBCRegions); CHKERRQ(ierr);
            /*PetscSynchronizedPrintf(comm, "bcId[%d] = %d\n", i, bcId[i]);*/
        }

        ierr = ISRestoreIndices(bcIS, &bcId);CHKERRQ(ierr);
        ierr = ISDestroy(&bcIS);CHKERRQ(ierr);

        PetscMalloc(numBC * sizeof(PetscInt), &bcField);
        for(i = 0; i < numBC; i++) bcField[i] = 0;

        PetscMalloc(numBC * sizeof(IS), &bcPointIS);
        for(i = 0; i < numBC; i++) 
        {
            PetscInt bcPointSize;
            ierr = DMGetStratumIS(dm, "Face Sets", bcId[i], &bcPointIS[i]);CHKERRQ(ierr);
            ierr = ISGetSize(bcPointIS[i], &bcPointSize);CHKERRQ(ierr);
            /*PetscSynchronizedPrintf(comm, "bcPointSize[%d] = %d\n", i, bcPointSize);*/
        }
        /*PetscSynchronizedFlush(comm,PETSC_STDOUT);*/

        /* Create a PetscSection with this data layout */
        ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);CHKERRQ(ierr);
        for(i = 0; i < numBC; i++)
        {
            ierr = ISDestroy(&bcPointIS[i]);CHKERRQ(ierr);
        }
        PetscFree(bcPointIS); CHKERRQ(ierr);

        /* Name the Field variables */
        ierr = PetscSectionSetFieldName(section, 0, "conVar");CHKERRQ(ierr);
        ierr = PetscSectionSetFieldName(section, 1, "grad");CHKERRQ(ierr);
        ierr = PetscSectionSetFieldName(section, 2, "limiter");CHKERRQ(ierr);
        ierr = PetscSectionSetFieldName(section, 3, "flux");CHKERRQ(ierr);
        DMSetDefaultSection(dm, section);
        PetscSectionDestroy(&section);
    }


    /* Initial conditions */
    ierr = DMGetGlobalVector(dm, &X);CHKERRQ(ierr);
    ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);

    SetInitialCondition(dm, localX);

    ierr = DMLocalToGlobalBegin(dm,localX,INSERT_VALUES, X); CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(dm,localX,INSERT_VALUES, X); CHKERRQ(ierr);

    ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(dm, &X);CHKERRQ(ierr);

    ierr = DMGetGlobalVector(dm, &X);CHKERRQ(ierr);
    ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);

    while(timeStep < maxTimeStep)
    {
        CalcConVarGradient(dm, X, gradStart);
        CalcLimiter(dm);



        timeStep++;
    }



    ierr = PetscObjectSetName((PetscObject) localX, "sol.");CHKERRQ(ierr);
    ierr = OutputVTK(dm, "sol.vtu", &viewer);CHKERRQ(ierr);
    ierr = VecView(localX, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(dm, &X);CHKERRQ(ierr);

    PetscFree(bcField);

    VecDestroy(&localX);
    VecDestroy(&X);

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