static char help[] = "Test load second order gmsh file\n\n";

#include <petscdmplex.h>

int main(int argc, char **argv)
{
  char           filename[256];
  DM             dm,cdm;
  PetscFE        fe;
  PetscFE        fe_old;
  PetscBool      interpolate = PETSC_TRUE;
  PetscBool      is_set = PETSC_FALSE;
  PetscReal      lo[3], hi[3];
  PetscInt       i, dim;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Load gmsh file",NULL);CHKERRQ(ierr);
  ierr = PetscOptionsString("-filename","Gmsh file name","",filename,filename,sizeof(filename),&is_set);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-interpolate","Interpolate the loaded mesh (default: TRUE)","",interpolate,&interpolate,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  
  if (is_set) {
    ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,filename,interpolate,&dm);CHKERRQ(ierr);
    ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
    ierr = DMGetBoundingBox(dm,lo,hi);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,"Old Bounding Box:\n");CHKERRQ(ierr);
    for (i=0;i<dim;i++) {
      ierr = PetscPrintf(PETSC_COMM_SELF,"  %d: lo = %g hi = %g\n",i,(double)lo[i],(double)hi[i]);CHKERRQ(ierr);
    }

    ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
    ierr = DMGetField(cdm,0,NULL,(PetscObject*)&fe_old);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) fe_old, "OldCoordinatesFE");CHKERRQ(ierr);
    ierr = PetscFEViewFromOptions(fe_old,NULL,"-old_fe_view");CHKERRQ(ierr);

    ierr = PetscFECreateLagrange(PETSC_COMM_SELF,dim,dim,PETSC_TRUE,2,PETSC_DETERMINE,&fe);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) fe, "NewCoordinatesFE");CHKERRQ(ierr);
    ierr = PetscFEViewFromOptions(fe,NULL,"-new_fe_view");CHKERRQ(ierr);

    ierr = DMProjectCoordinates(dm, fe);CHKERRQ(ierr);
    ierr = DMGetBoundingBox(dm,lo,hi);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,"New Bounding Box:\n");CHKERRQ(ierr);
    for (i=0;i<dim;i++) {
      ierr = PetscPrintf(PETSC_COMM_SELF,"  %d: lo = %g hi = %g\n",i,(double)lo[i],(double)hi[i]);CHKERRQ(ierr);
    }

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