static char help[] = "test1.c: test of the new generation ExodusII readers: loads and displays a mesh\n\n";

#include <petscsys.h>
#include <petscdmmesh.hh>
#include "meshexodus.h"
#include "exodusII.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main (int argc,char ** argv) {
  DM              dmSeq,dm;
  PetscBool       inflag;
  char            infilename[PETSC_MAX_PATH_LEN+1];
  PetscErrorCode  ierr;
  PetscInt        my_num_cells,my_num_vertices;
  int             rank,numproc;
  int             CPU_word_size = 0;
  int             IO_word_size  = 0;
  int             exoid;
  float           version;


  ierr = PetscInitialize(&argc,&argv,(char *)0,help);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&numproc);
  ierr = PetscOptionsGetString(PETSC_NULL,"-i",infilename,PETSC_MAX_PATH_LEN,&inflag);CHKERRQ(ierr);
  if (!inflag) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"No file name given\n");CHKERRQ(ierr);
    PetscFinalize();
    return 0;
  }

  /*
    Reads a mesh
  */  
  if (rank == 0) {
    exoid = ex_open(infilename,EX_READ,&CPU_word_size,&IO_word_size,&version);CHKERRQ(!exoid);
  }
  ierr = myDMMeshCreateExodusNG(PETSC_COMM_WORLD,exoid,&dmSeq);CHKERRQ(ierr);
  if (rank == 0) {
    ierr = ex_close(exoid);CHKERRQ(ierr);
  }
  
  /*
    Partition it if necessary
  */
  if (numproc >1) {
    ierr = DMMeshCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
    ierr = DMMeshDistribute(dmSeq,PETSC_NULL,&dm);CHKERRQ(ierr);
    ierr = DMDestroy(&dmSeq);CHKERRQ(ierr);
  } else {
    dm = dmSeq;
  }

  /*
    Displays the body mesh
  */
  {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"\ndm:\n");CHKERRQ(ierr);
    ierr = DMView(dm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ALE::Obj<PETSC_MESH_TYPE> mesh;
    ierr = DMMeshGetMesh(dm,mesh);CHKERRQ(ierr);
    my_num_cells    = mesh->heightStratum(0)->size();
    my_num_vertices = mesh->depthStratum(0)->size();
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%i]: Number of vertices in mesh: \t%i\n",rank,my_num_vertices);CHKERRQ(ierr);
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%i]: Number of cells in mesh: \t%i\n",rank,my_num_cells);CHKERRQ(ierr);
    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD);CHKERRQ(ierr);
    
    mesh->view("mesh");
    const ALE::Obj<PETSC_MESH_TYPE::label_type>& cellSets = mesh->getLabel("Cell Sets");
    const ALE::Obj<PETSC_MESH_TYPE::label_type>& vertexSets = mesh->getLabel("Vertex Sets");
    cellSets->view("Cell sets");
    vertexSets->view("Vertex sets");
  }
  ierr = DMDestroy(&dm);CHKERRQ(ierr);  
  ierr = PetscFinalize();
  return 0;
}