static char help[] = "\n";

#include <petscdmplex.h>
#include <petscsf.h>
#include <exodusII.h>
#include <petsc-private/dmimpl.h>     /*I      "petscdm.h"     I*/
#include <petscsf.h>
#include <petscviewerhdf5.h>
#include <petscdmplex.h>   

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
	PetscErrorCode      ierr;
	DM                  dm,distDM;
	char                ifilename[PETSC_MAX_PATH_LEN],ofilename[PETSC_MAX_PATH_LEN],prefix[PETSC_MAX_PATH_LEN];
  char                *ext;
	PetscBool           flg;
	int                 CPU_word_size = 0,IO_word_size = 0,exoid;
	float               version;  
	PetscMPIInt         numproc,rank;
	PetscViewer         hdf5Viewer;
	PetscInt            dim;
	PetscInt            numFields = 2;
	PetscInt            numComp[2] = {1,1};
	PetscInt            numDof[6] = {1,0,0,
	                                 0,0,1}; /*{Vertex,Edge,Cell} */
	PetscInt            bcFields[1] = {0},numBC=0;
	IS                  bcPoints[1] = {NULL},*ISlist;
	PetscSection        seqSection,distSection;
	Vec                 V;
	PetscSF             pointSF;
	


	ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
	ierr = PetscOptionsGetString(NULL,"-i",ifilename,sizeof(ifilename),&flg);CHKERRQ(ierr);
  ierr = PetscStrrchr(ifilename,'.',&ext);
  ierr = PetscStrncpy(prefix,ifilename,ext-ifilename);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"prefix is %s\n",prefix);CHKERRQ(ierr);
	
	ierr = MPI_Comm_size(PETSC_COMM_WORLD,&numproc);
	ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

	/* 
	  create DM from exo
	*/

	ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,ifilename,PETSC_TRUE,&dm);CHKERRQ(ierr);
	ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr);
	ierr = DMPlexCreateSection(dm,dim,numFields,numComp,numDof,numBC,bcFields,bcPoints,NULL,&seqSection);CHKERRQ(ierr);
	ierr = DMSetDefaultSection(dm,seqSection);CHKERRQ(ierr);
	
	if (numproc == 1) {
	  ierr = PetscSNPrintf(ofilename,FILENAME_MAX,"%s_seq.h5",prefix);CHKERRQ(ierr);
	  ierr = PetscPrintf(PETSC_COMM_WORLD,"Writing to %s\n",ofilename);CHKERRQ(ierr);
    ierr = PetscViewerHDF5Open(PetscObjectComm((PetscObject) dm),ofilename,FILE_MODE_WRITE,&hdf5Viewer); CHKERRQ(ierr);
    ierr = DMView_Plex(dm,hdf5Viewer); CHKERRQ(ierr);
    ierr = PetscViewerHDF5SetTimestep(hdf5Viewer,0);

  	ierr = DMGetGlobalVector(dm,&V);CHKERRQ(ierr);
    ierr = VecSet(V,0.0);CHKERRQ(ierr);
    ierr = VecView(V,hdf5Viewer);CHKERRQ(ierr);
    ierr = PetscViewerHDF5SetTimestep(hdf5Viewer,1);
    ierr = VecSet(V,1.0);CHKERRQ(ierr);
    ierr = VecView(V,hdf5Viewer);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&hdf5Viewer);CHKERRQ(ierr);
	}


	if (numproc > 1) {
		ierr = DMPlexDistribute(dm,PETSC_NULL,0,&pointSF,&distDM);CHKERRQ(ierr);
		ierr = DMPlexCreateSection(distDM,dim,numFields,numComp,numDof,numBC,bcFields,bcPoints,NULL,&distSection);CHKERRQ(ierr);
		ierr = DMSetDefaultSection(distDM,distSection);CHKERRQ(ierr);
	  ierr = PetscSNPrintf(ofilename,FILENAME_MAX,"%s_dist.h5",prefix);CHKERRQ(ierr);
		ierr = PetscViewerHDF5Open(PetscObjectComm((PetscObject) distDM),ofilename,FILE_MODE_WRITE,&hdf5Viewer); CHKERRQ(ierr);
		ierr = DMView(distDM,hdf5Viewer); CHKERRQ(ierr);

    /* 
      I suspect that I also need to save time values, but I can't figure out how to do that.
    */
  	ierr = DMGetGlobalVector(distDM,&V);CHKERRQ(ierr);
    ierr = VecSet(V,0.0);CHKERRQ(ierr);
    ierr = VecView(V,hdf5Viewer);CHKERRQ(ierr);
    ierr = PetscViewerHDF5SetTimestep(hdf5Viewer,1);
    ierr = VecSet(V,1.0);CHKERRQ(ierr);
    ierr = VecView(V,hdf5Viewer);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(distDM,&V);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&hdf5Viewer);CHKERRQ(ierr);
		ierr = DMDestroy(&distDM);CHKERRQ(ierr);
	}


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