static char help[] = "Test the function \n\n";

#include <petscdmplex.h>
#include <petscds.h>
#include <petscfe.h>

static PetscErrorCode get_volume_DM_v2(DM dm, PetscFEGeom *cellgeom, PetscReal *volG)
{
  PetscDS         prob;
  PetscFE         fe;
  PetscQuadrature quad;
  const PetscReal *quadWeights;
  PetscInt        c,q,qNc,Nq,cellHeight=0,cStart,cEnd;
  PetscReal       vol=0.;

  PetscFunctionBeginUser;
  DMGetDS(dm,&prob);
  PetscDSGetDiscretization(prob,0,(PetscObject*) &fe);
  PetscFEGetQuadrature(fe,&quad);
  PetscQuadratureView(quad,PETSC_VIEWER_STDOUT_WORLD);
  PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
  PetscCheck(cellgeom->numPoints == Nq, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Cell geometry does not match FE quadrature");

  DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
  for (c = cStart; c < cEnd; ++c) {
    for (q = 0; q < Nq; ++q) {
      const PetscReal detJ = cellgeom->detJ[(c - cStart) * Nq + q];

      PetscCheck(detJ >= 0., PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "ATTENTION, volume non correct !!\n\n");
      vol += detJ*quadWeights[q*qNc];
    }
  }

  PetscCallMPI(MPI_Allreduce(&vol, volG, 1,MPIU_REAL, MPI_SUM,PETSC_COMM_WORLD));
  PetscFunctionReturn(0);
}

int main(int argc, char **argv) {
  DM              dm;
  PetscFE         fe;
  PetscReal       vol;
  PetscInt        dim=2,NumComp=1;
  PetscInt        cells[3]={2,2,2};
  PetscBool       isSimplice=PETSC_TRUE;
  DMField         coordField;
  IS              cellIS;
  PetscFEGeom    *fegeom;
  PetscQuadrature quad;

  PetscCall(PetscInitialize(&argc,&argv,(char *)0,help));

  PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, isSimplice, cells, NULL, NULL, NULL, PETSC_TRUE, &dm));
  PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim,NumComp,PETSC_FALSE,"heat_",PETSC_DECIDE,&fe));
  PetscCall(PetscFEViewFromOptions(fe,NULL,"-fe_field_view"));
  PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
  PetscFEGetQuadrature(fe,&quad);
  PetscQuadratureView(quad,PETSC_VIEWER_STDOUT_WORLD);
  PetscCall(PetscFEDestroy(&fe));
  PetscCall(DMCreateDS(dm));

  PetscCall(DMGetCoordinateField(dm, &coordField));
  PetscCall(DMGetStratumIS(dm, "depth", dim, &cellIS));
  PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, PETSC_FALSE, &fegeom));
  PetscCall(ISDestroy(&cellIS));

  PetscCall(get_volume_DM_v2(dm, fegeom, &vol));
  PetscPrintf(PETSC_COMM_WORLD, "Volume %f\n",vol);

  PetscCall(PetscFEGeomDestroy(&fegeom));
  PetscCall(DMDestroy(&dm));

  PetscCall(PetscFinalize());
  return 0;
}