static char help[] = "Trying to understand DMPlexComputeCellGeometryAffineFEM";
#include <petsc.h>

int main(int argc,char **argv)
{
  DM             dm;
  char           filename[2048];
  Vec            coordLoc;
  PetscReal      v0[3];
  PetscReal      J[9];
  PetscReal      invJ[9];
  PetscReal      detJ;
  PetscInt       dim;

  PetscCall(PetscInitialize(&argc,&argv,NULL,help));
  PetscOptionsBegin(PETSC_COMM_WORLD,"","TestDMPlexComputeCellGeometryAffineFEM options","none");
  PetscCall(PetscOptionsString("-i","filename to read","",filename,filename,sizeof(filename),NULL));
  PetscOptionsEnd();

  PetscCall(PetscPrintf(PETSC_COMM_WORLD," filename %s\n",filename));
  PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, NULL, PETSC_FALSE, &dm));

  PetscCall(DMGetCoordinatesLocal(dm,&coordLoc));
  PetscCall(VecView(coordLoc,PETSC_VIEWER_STDOUT_WORLD));

  PetscCall(DMGetDimension(dm,&dim));
  PetscCall(DMPlexComputeCellGeometryAffineFEM(dm,0,v0,J,invJ,&detJ));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,"v0\n"));
  PetscCall(PetscRealView(dim,v0,PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,"J\n"));
  for (int i = 0; i < dim; i++){
    PetscCall(PetscRealView(dim,&J[i*dim],PETSC_VIEWER_STDOUT_WORLD));
  }
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,"invJ\n"));
  for (int i = 0; i < dim; i++){
    PetscCall(PetscRealView(dim,&invJ[i*dim],PETSC_VIEWER_STDOUT_WORLD));
  }
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,"detJ : %g\n",detJ));


  PetscCall(DMDestroy(&dm));
  PetscCall(PetscFinalize());
  return 0;
}