static char help[] = "EXOTest5 \n\n";

#include <petsc.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv) {
  DM                dm,dmDist,cdm;
  PetscSF           migrationSF;
  char              filename[2048];
  PetscBool         interpolate=PETSC_FALSE,flag;
  
  PetscInt          i,c,dim,cdim,numStrata,closureSize;
  PetscInt          pStart,pEnd,cStart,cEnd,vStart,vEnd,eStart,eEnd,step,nstep=1;
  PetscSection      section,coordSection;
  PetscInt          dof=1,totalDofs[] = {1,0,0}; // Vertex,edge,face
  
  Vec               global,coord;
  PetscMPIInt       rank,numproc;
  PetscErrorCode    ierr;
  PetscScalar       *cval,*val;

  PetscInitialize(&argc,&argv,NULL,help);
  
  PetscOptionsBegin(PETSC_COMM_WORLD,"","options","none");
  ierr = PetscOptionsString("-i","filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-dof","number of dof per vertex","",dof,&dof,&flag);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-nstep","number of time steps","",nstep,&nstep,&flag);CHKERRQ(ierr);
  PetscOptionsEnd();
  
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  MPI_Comm_size(PETSC_COMM_WORLD,&numproc);

  ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,filename,interpolate,&dm);CHKERRQ(ierr);
  
  /* Making and set default section */
  ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
  ierr = DMPlexGetDepth(dm,&numStrata);CHKERRQ(ierr);


  ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
  ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr);
  ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"\ndim         : %d\n",dim);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"depth       : %d\n",numStrata);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"cell stratum: %d-%d\n",cStart,cEnd);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"edge stratum: %d-%d\n",eStart,eEnd);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"vert stratum: %d-%d\n\n",vStart,vEnd);CHKERRQ(ierr);

  
  ierr = PetscSectionCreate(PETSC_COMM_WORLD,&section);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(section,pStart,pEnd);CHKERRQ(ierr);
  
  for (i = vStart; i < vEnd; ++i) {
    ierr = PetscSectionSetDof(section,i,totalDofs[0]*dof);CHKERRQ(ierr);
  }
  if (cStart != eStart){
    for (i = eStart; i < eEnd; ++i) {
      ierr = PetscSectionSetDof(section,i,totalDofs[1]*dof);CHKERRQ(ierr);
    }
  }
  for (i = cStart; i < cEnd; ++i) {
    ierr = PetscSectionSetDof(section,i,totalDofs[2]*dof);CHKERRQ(ierr);
  }
  ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
  ierr = DMSetDefaultSection(dm,section);CHKERRQ(ierr);
  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);

  /* Distribute */
  ierr = DMSetUseNatural(dm,PETSC_TRUE);CHKERRQ(ierr);
  ierr = DMPlexDistribute(dm,0,&migrationSF,&dmDist);CHKERRQ(ierr);
  if (dmDist) {
    ierr = DMDestroy(&dm);CHKERRQ(ierr);
    dm   = dmDist;
  }

  ierr = DMGetGlobalVector(dm,&global);CHKERRQ(ierr);
  //ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
  //ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr);
  //ierr = DMGetCoordinates(dm,&coord);CHKERRQ(ierr);

  ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
  ierr = DMPlexGetDepthStratum(dm,1,&eStart,&eEnd);CHKERRQ(ierr);
  ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);
  for (i = pStart; i < pEnd; i++) {
    cdim = 0;
    //cval = NULL;
    //ierr = DMPlexVecGetClosure(cdm,NULL,coord,i,&cdim,&cval);CHKERRQ(ierr);
    //ierr = DMPlexVecGetClosure(dm,coordSection,coord,i,&cdim,&cval);CHKERRQ(ierr);
    ierr = DMPlexVecGetClosure(dm,NULL,global,i,&cdim,&cval);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"point %d cdim %d\n",i,cdim);CHKERRQ(ierr);
    for (step = 0; step < cdim; step++) {
      cval[step] = i+step*.0001;
    }
    ierr = DMPlexVecRestoreClosure(dm,NULL,global,i,&cdim,&cval);CHKERRQ(ierr);
    //ierr = DMPlexVecRestoreClosure(dm,coordSection,coord,i,&cdim,&cval);CHKERRQ(ierr);
    //ierr = DMPlexVecRestoreClosure(cdm,NULL,coord,i,&cdim,&cval);CHKERRQ(ierr);
  }
  ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);
  ierr = DMRestoreGlobalVector(dm,&global);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  return 0;
}