static char help[] = "testing periodic bc\n";
#include <petsc.h>

#define nvar 2

int main(int argc, char **argv)
{
   PetscErrorCode ierr;
   DM             dm, dmDist = NULL;
   PetscSection   sec;
   PetscInt       dim = 2, numFields = 1, overlap = 1, numBC, i;
   PetscInt       numComp[1]; // number of components per field
   PetscInt       numDof[3];
   PetscBool      interpolate = PETSC_TRUE, useCone = PETSC_TRUE, useClosure = PETSC_TRUE;
   PetscMPIInt    rank, size;

   ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
   MPI_Comm_size(PETSC_COMM_WORLD, &size);
   
   // read mesh from gmsh file
   ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, "ug_periodic.msh", interpolate, &dm); CHKERRQ(ierr);
   // set periodicity information (periodic in x and y)
   //const DMBoundaryType bd[] = {DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC};
   //ierr = DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, bd); CHKERRQ(ierr); 

   // set partitioner type to Parmetis
   PetscPartitioner part;
   ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr);
   ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE); CHKERRQ(ierr);
   ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 

   // distribute mesh over processes;
   ierr = DMSetBasicAdjacency(dm, useCone, useClosure); CHKERRQ(ierr);
   ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);
   if(dmDist) 
   {
      ierr = DMDestroy(&dm); CHKERRQ(ierr);
      dm = dmDist;
   }

   ierr = DMSetNumFields(dm, numFields); CHKERRQ(ierr);
   ierr = DMSetAdjacency(dm, 0, useCone, useClosure); CHKERRQ(ierr);
   // setup FV
   PetscFV fv;
   ierr = PetscFVCreate(PETSC_COMM_WORLD, &fv); CHKERRQ(ierr);
   ierr = PetscFVSetType(fv, PETSCFVUPWIND); CHKERRQ(ierr);
   ierr = PetscFVSetNumComponents(fv, nvar); CHKERRQ(ierr);
   ierr = PetscFVSetSpatialDimension(fv, dim); CHKERRQ(ierr);
   ierr = PetscFVSetFromOptions(fv); CHKERRQ(ierr);
   ierr = PetscFVSetUp(fv); CHKERRQ(ierr); 

   // create field for solution u 
   numComp[0] = nvar; // six components
   for(i=0; i<numFields*(dim+1); ++i) numDof[i] = 0;
   // define u at cell centers
   numDof[0*(dim+1)+dim] = nvar;
   
   numBC = 0;

   // create a PetscSection with this layout
   ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC,
                              NULL, NULL, NULL, NULL, &sec); CHKERRQ(ierr);
   // name the field and its components
   ierr = PetscSectionSetFieldName(sec, 0, "U"); CHKERRQ(ierr);

   // set the section for dm
   ierr = DMSetLocalSection(dm, sec); CHKERRQ(ierr); 

   ierr = DMSetUp(dm); CHKERRQ(ierr);
   ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);   

   // get geometric data
   DM dmFace, dmCell, dummy;
   Vec cellGeom, faceGeom;
   const PetscScalar *fgeom, *cgeom; // arrays corresponding to vectors
   PetscFVCellGeom *cg; // struct with cell geometry information
   PetscFVFaceGeom *fg; // struct with face geometry information
   ierr = DMPlexGetDataFVM(dm, fv, &cellGeom, &faceGeom, &dummy); CHKERRQ(ierr);
   ierr = VecGetDM(faceGeom, &dmFace); CHKERRQ(ierr);
   ierr = VecGetDM(cellGeom, &dmCell); CHKERRQ(ierr);
   ierr = VecGetArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
   ierr = VecGetArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
   // print geometric information
   PetscInt f, fStart, fEnd, ss;
   const PetscInt *supp; // support
   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);
   for(f=fStart; f<fEnd; ++f)
   {
      ierr = DMPlexGetSupportSize(dm, f, &ss); CHKERRQ(ierr);
      ierr = DMPlexGetSupport(dm, f, &supp); CHKERRQ(ierr);
      // face geometry
      ierr = DMPlexPointLocalRead(dmFace, f, fgeom, &fg); CHKERRQ(ierr);
      PetscPrintf(PETSC_COMM_WORLD, "face: %d, centroid = %f, %f, normal = %f, %f\n", f, fg->centroid[0], fg->centroid[1], fg->normal[0], fg->normal[1]);
      PetscReal len = PetscSqrtReal(PetscPowReal(fg->normal[0],2) + PetscPowReal(fg->normal[1],2));
      if(PetscAbs(len - 2.5) > 1.0e-10)
      {
         PetscPrintf(PETSC_COMM_SELF,"===> Face length incorrect = %f, should be 2.5\n", len);
      }
      // cell geometry
      for(PetscInt i=0; i<ss; ++i)
      {
         ierr = DMPlexPointLocalRead(dmCell, supp[i], cgeom, &cg); CHKERRQ(ierr); 
         PetscPrintf(PETSC_COMM_WORLD, "support[%d] = %d, cent = %f, %f, area = %f\n", i, supp[i], cg->centroid[0], cg->centroid[1], cg->volume);
      }
      PetscPrintf(PETSC_COMM_WORLD, "\n");
   }

   // cleanup
   ierr = PetscFVDestroy(&fv); CHKERRQ(ierr);
   ierr = VecRestoreArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
   ierr = VecRestoreArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
   ierr = DMDestroy(&dm); CHKERRQ(ierr);
   ierr = PetscFinalize(); CHKERRQ(ierr);
   return ierr;
}