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

#include <petsc.h>
#include "exodusII.h"

#undef __FUNCT__
#define __FUNCT__ "EXOGetVarIndex"
PetscErrorCode EXOGetVarIndex(int exoid,ex_entity_type obj_type,const char name[],int *varIndex){
  PetscErrorCode    ierr;
  int               exoerr,num_vars,i,j;
  char              ext_name[MAX_STR_LENGTH+1],var_name[MAX_STR_LENGTH+1];
  const             int num_suffix = 5;
  char              *suffix[5];
  PetscBool         flg;

  PetscFunctionBeginUser;
  suffix[0] = (char *) "";
  suffix[1] = (char *) "_X";
  suffix[2] = (char *) "_XX";
  suffix[3] = (char *) "_1";
  suffix[4] = (char *) "_11";

  exoerr = ex_get_variable_param(exoid,obj_type,&num_vars);
  for (i = 0; i < num_vars; i++) {
	  exoerr = ex_get_variable_name(exoid,obj_type,i+1,var_name);
    for (j = 0; j < num_suffix; j++){
      ierr = PetscStrncpy(ext_name,name,MAX_STR_LENGTH);CHKERRQ(ierr);
      ierr = PetscStrncat(ext_name,suffix[j],MAX_STR_LENGTH);CHKERRQ(ierr);
      ierr = PetscStrcasecmp(ext_name,var_name,&flg);
      if (flg) {
        break;
      }       
    }
    if (flg) {
      *varIndex = i+1;
      break;
    }       
  }
  CHKMEMQ;
  if (varIndex != 0) {
    PetscFunctionReturn(1);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "VecViewNodal_PlexEXO"
PetscErrorCode VecViewNodal_PlexEXO(Vec v,DM dm,int exoid,int step){
  MPI_Comm          comm;
  PetscMPIInt       rank,numproc;
  Vec               vNatural,vComp;
  PetscReal         const *varray;
  PetscInt          xs,xm,bs,c;
  PetscBool         useNatural;
  const char        *vecname;
  int               offset;
  IS                compIS;
  PetscErrorCode    ierr;

  PetscFunctionBeginUser;
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&numproc);
  ierr = MPI_Comm_rank(comm,&rank);

  ierr = DMGetUseNatural(dm,&useNatural);CHKERRQ(ierr);
  if (numproc > 1 && useNatural) {
    ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);
    ierr = DMPlexGlobalToNaturalBegin(dm,v,vNatural);CHKERRQ(ierr);
    ierr = DMPlexGlobalToNaturalEnd(dm,v,vNatural);CHKERRQ(ierr);
  } else {
    vNatural = v;
  }

  /*
    Get the location of the variable in the exodus file
  */
  if (rank == 0) {
    ierr = PetscObjectGetName((PetscObject) v,&vecname);CHKERRQ(ierr);
    ierr = EXOGetVarIndex(exoid,EX_NODAL,vecname,&offset);
  }
  ierr = MPI_Bcast(&offset,1,MPI_INT,0,comm);
  if (!offset) {
    SETERRQ1(comm,PETSC_ERR_FILE_UNEXPECTED,"Cannot locate nodal variable %s in exodus file\n",vecname);
  }

  /*
    Write local chunk of the result in the exodus file
    exodus stores each component of a vector-valued field as a separate variable.
    We assume that they are stored sequentially
  */
  ierr = VecGetOwnershipRange(vNatural,&xs,&xm);CHKERRQ(ierr);
  ierr = VecGetBlockSize(vNatural,&bs);CHKERRQ(ierr);
  if (bs == 1) {
      ierr = VecGetArrayRead(vNatural,&varray);CHKERRQ(ierr);
      ierr = ex_put_partial_var(exoid,step,EX_NODAL,offset,1,xs+1,xm-xs,varray);
      ierr = VecRestoreArrayRead(vNatural,&varray);CHKERRQ(ierr);
  } else {
    ierr = ISCreateStride(comm,(xm-xs)/bs,xs,bs,&compIS);CHKERRQ(ierr);
    for (c = 0; c < bs; ++c) {
      ierr = ISStrideSetStride(compIS,(xm-xs)/bs,xs+c,bs);CHKERRQ(ierr);
      ierr = VecGetSubVector(vNatural,compIS,&vComp);CHKERRQ(ierr);
      ierr = VecGetArrayRead(vComp,&varray);CHKERRQ(ierr);
      ierr = ex_put_partial_var(exoid,step,EX_NODAL,offset+c,1,xs/bs+1,(xm-xs)/bs,varray);
      ierr = VecRestoreArrayRead(vComp,&varray);CHKERRQ(ierr);
      ierr = VecRestoreSubVector(vNatural,compIS,&vComp);CHKERRQ(ierr);
    }
    ierr = ISDestroy(&compIS);CHKERRQ(ierr);    
  }
  if (numproc > 1 && useNatural) {
    ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "VecLoadNodal_PlexEXO"
PetscErrorCode VecLoadNodal_PlexEXO(Vec v,DM dm,int exoid,int step){
  MPI_Comm          comm;
  PetscMPIInt       rank,numproc;
  Vec               vNatural,vComp;
  PetscReal         *varray;
  PetscInt          xs,xm,bs,c;
  PetscBool         useNatural;
  const char        *vecname;
  int               offset;
  IS                compIS;
  PetscErrorCode    ierr;

  PetscFunctionBeginUser;
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&numproc);
  ierr = MPI_Comm_rank(comm,&rank);

  ierr = DMGetUseNatural(dm,&useNatural);CHKERRQ(ierr);
  if (numproc > 1) {
    ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);
  } else {
    vNatural = v;
  }

  /*
    Get the location of the variable in the exodus file
  */
  if (rank == 0) {
    ierr = PetscObjectGetName((PetscObject) v,&vecname);CHKERRQ(ierr);
    ierr = EXOGetVarIndex(exoid,EX_NODAL,vecname,&offset);
  }
  ierr = MPI_Bcast(&offset,1,MPI_INT,0,comm);
  if (!offset) {
    SETERRQ1(comm,PETSC_ERR_FILE_UNEXPECTED,"Cannot locate nodal variable %s in exodus file\n",vecname);
  }

  /* 
    Read local chunk from the file
  */
  ierr = VecGetOwnershipRange(vNatural,&xs,&xm);CHKERRQ(ierr);
  ierr = VecGetBlockSize(vNatural,&bs);CHKERRQ(ierr);

  if (bs == 1) {
      ierr = VecGetArray(vNatural,&varray);CHKERRQ(ierr);
      ierr = ex_get_partial_var(exoid,step,EX_NODAL,offset,1,xs+1,xm-xs,varray);
      ierr = VecRestoreArray(vNatural,&varray);CHKERRQ(ierr);
  } else {
    ierr = ISCreateStride(comm,(xm-xs)/bs,xs,bs,&compIS);CHKERRQ(ierr);
    for (c = 0; c < bs; ++c) {
      ierr = ISStrideSetStride(compIS,(xm-xs)/bs,xs+c,bs);CHKERRQ(ierr);
      ierr = VecGetSubVector(vNatural,compIS,&vComp);CHKERRQ(ierr);
      ierr = VecGetArray(vComp,&varray);CHKERRQ(ierr);
      ierr = ex_get_partial_var(exoid,step,EX_NODAL,offset+c,1,xs/bs+1,(xm-xs)/bs,varray);
      ierr = VecRestoreArray(vComp,&varray);CHKERRQ(ierr);
      ierr = VecRestoreSubVector(vNatural,compIS,&vComp);CHKERRQ(ierr);
    }
    ierr = ISDestroy(&compIS);CHKERRQ(ierr);    
	}
  if (numproc > 1 && useNatural) {
  	ierr = VecSet(v,1000.0);
    ierr = DMPlexNaturalToGlobalBegin(dm,vNatural,v);CHKERRQ(ierr);
    ierr = DMPlexNaturalToGlobalEnd(dm,vNatural,v);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv) {  
  char              filename[2048];
  PetscBool         interpolate=PETSC_TRUE;
  PetscInt          p,pStart,pEnd;
  PetscInt          d,sdim,numStrata;


  DM                dmUAS,dmU,dmAlpha,dmSigma,dmDist;
  PetscSection      sectionUAS;
  Vec               UAS,U,Alpha,Sigma;
  /* 
    dof layout: cell,face,edge,vertex 
    ordered by increasing depth in the DAG
  */
  PetscInt          dofU2D[]     = {2,0,0};
  PetscInt          dofAlpha2D[] = {1,0,0};
  PetscInt          dofSigma2D[] = {0,0,3};

  PetscInt          dofU3D[]     = {3,0,0,0};
  PetscInt          dofAlpha3D[] = {1,0,0,0};
  PetscInt          dofSigma3D[] = {0,0,0,6};

  PetscInt          *dofU,*dofAlpha,*dofSigma;

  PetscMPIInt       rank,numproc;
  PetscErrorCode    ierr;
  PetscInt          i;


  /* 
    exodus I/O related variables
  */
  int               CPU_word_size;
  int               IO_word_size;
  int               exoid;
  float             EXO_version;
  int               EXO_mode,nstep=1,step=1;
  MPI_Info          mpi_info = MPI_INFO_NULL;


  PetscInitialize(&argc,&argv,NULL,help);
  ex_opts(EX_VERBOSE+EX_DEBUG);
  
  PetscOptionsBegin(PETSC_COMM_WORLD,"","options","none");
  ierr = PetscOptionsString("-i","filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
  PetscOptionsEnd();
  
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  MPI_Comm_size(PETSC_COMM_WORLD,&numproc);

  ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,filename,interpolate,&dmUAS);CHKERRQ(ierr);
  
  /* Making and set default section */
  ierr = PetscSectionCreate(PETSC_COMM_WORLD,&sectionUAS);CHKERRQ(ierr);
  ierr = PetscSectionSetNumFields(sectionUAS,3);CHKERRQ(ierr);
  ierr = DMPlexGetChart(dmUAS,&pStart,&pEnd);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(sectionUAS,pStart,pEnd);CHKERRQ(ierr);

  ierr = DMGetDimension(dmUAS,&sdim);CHKERRQ(ierr);
  ierr = PetscSectionSetFieldComponents(sectionUAS,0,sdim);CHKERRQ(ierr);
  ierr = PetscSectionSetFieldComponents(sectionUAS,1,1);CHKERRQ(ierr);
  ierr = PetscSectionSetFieldComponents(sectionUAS,2,sdim*(sdim+1)/2);CHKERRQ(ierr);

  if (sdim == 2) {
    dofU     = dofU2D;
    dofAlpha = dofAlpha2D;
    dofSigma = dofSigma2D;
  } else {
    dofU     = dofU3D;
    dofAlpha = dofAlpha3D;
    dofSigma = dofSigma3D;
  }

  for (d = 0; d < sdim+1; d++) {
    ierr = DMPlexGetDepthStratum(dmUAS,d,&pStart,&pEnd);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"depth %d: range %d-%d layout %d %d %d \n",d,pStart,pEnd,dofU[d],dofAlpha[d],dofSigma[d]);
    for (p = pStart; p < pEnd; p++) {
      ierr = PetscSectionSetDof(sectionUAS,p,dofU[d]+dofAlpha[d]+dofSigma[d]);CHKERRQ(ierr);
      ierr = PetscSectionSetFieldDof(sectionUAS,p,0,dofU[d]);CHKERRQ(ierr);
      ierr = PetscSectionSetFieldDof(sectionUAS,p,1,dofAlpha[d]);CHKERRQ(ierr);
      ierr = PetscSectionSetFieldDof(sectionUAS,p,2,dofSigma[d]);CHKERRQ(ierr);
    }
  }
  ierr = PetscSectionSetUp(sectionUAS);CHKERRQ(ierr);
  ierr = DMSetDefaultSection(dmUAS,sectionUAS);CHKERRQ(ierr);
  ierr = DMSetFromOptions(dmUAS);CHKERRQ(ierr);

  ierr = PetscSectionView(sectionUAS,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  
  ierr = DMDestroy(&dmUAS);CHKERRQ(ierr);
  PetscFinalize();
  return 0;
}