static char help[] = "Help message \n\n";

#define PETSC_APPLE_FRAMEWORK

#include <petsc.h>
#include <stdio.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv) {
  DM                dm,dmDist;
  PetscSF           migrationSF;
  PetscInt          pStart,pEnd,cStart,cEnd,vStart,vEnd,eStart,eEnd;
  PetscInt          numDofs[] = {2,0,0}; // Vertex, edge, face
  PetscSection      section;
  PetscInt          vecSize,v,e,c;
  Vec               globalVec;
  char              filename[2048];
  PetscBool         interpolate=PETSC_FALSE;
  PetscErrorCode    ierr;
  
  ierr = PetscInitialize(&argc,&argv,NULL,help);
  
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","options","none");CHKERRQ(ierr);
  ierr = PetscOptionsString("-i","filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-interpolate","Generate intermediate mesh elements","",interpolate,&interpolate,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  
  ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,filename,interpolate,&dm);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;
  }

  /* Create Section */
  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 = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(section,pStart,pEnd);CHKERRQ(ierr);
  
  for (v = vStart; v < vEnd; ++v) {
    ierr = PetscSectionSetDof(section, v, numDofs[0]);CHKERRQ(ierr);
  }
  if (cStart != eStart){
    for (e = eStart; e < eEnd; ++e) {
      ierr = PetscSectionSetDof(section, e, numDofs[1]);CHKERRQ(ierr);
    }
  }
  for (c = cStart; c < cEnd; ++c) {
    ierr = PetscSectionSetDof(section, c, numDofs[2]);CHKERRQ(ierr);
  }
  ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
  ierr = DMSetDefaultSection(dm,section);CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&section);
  
  /* Get global Vector */
  ierr = DMGetGlobalVector(dm, &globalVec);CHKERRQ(ierr);
//  ierr = DMGetGlobalVector(dm, &globalVec);CHKERRQ(ierr); // Uncomment this
  ierr = VecGetSize(globalVec, &vecSize);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "\nGlobal VecSize: %d\n",vecSize);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(dm, &globalVec);CHKERRQ(ierr);
  
  ierr = DMDestroy(&dm);
  ierr = PetscFinalize();
  return ierr;

}