[petsc-dev] DMPlexComputeGeometryFEM

[Yann Jobic]

Hi all,

I made a very simple example which uses the DMPlexComputeGeometryFEM 
function.

However, i've got an error :
./testComputeGeomFEM -heat_petscspace_degree 0
   Quadrature of order 1 on 1 points (dim 2)
[0]PETSC ERROR: --------------------- Error Message 
--------------------------------------------------------------
[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 7
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.1, unknown
[0]PETSC ERROR: 
/home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/testComputeGeomFEM 
on a  named luke by jobic Fri Nov 25 14:08:34 2022
[0]PETSC ERROR: Configure options 
--prefix=/local/lib/petsc/3.18/p0/gcc/nompi --with-mpi=0 
--with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est 
--download-triangle=1 --with-single-library=0 --with-large-file-io=1 
--with-shared-libraries=0 -CFLAGS=" -g -O0" -CXXFLAGS=" -g -O0" 
-FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc
[0]PETSC ERROR: #1 DMPlexComputeCellGeometryFEM() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:2158
[0]PETSC ERROR: #2 DMPlexComputeGeometryFEM() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:2524
[0]PETSC ERROR: #3 main() at 
/home/jobic/projet/fe-utils/petsc/fetools/src/testComputeGeomFEM.c:31
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -heat_petscspace_degree 0
[0]PETSC ERROR: ----------------End of Error Message -------send entire 
error message to petsc-maint at mcs.anl.gov----------

Do you know where it comes from ?

Thanks,

Yann
-------------- next part --------------
static char help[] = "Test the function \n\n";

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

extern double get_volume_DM_v2(DM dm, Vec cellgeom);

int main(int argc, char **argv) {
  DM              dm;
  Vec             cellgeom;      /* solution, residual vectors */
  PetscFE         fe;
  PetscReal       vol;
  PetscInt        dim=2,NumComp=1;
  PetscInt        cells[3]={2,2,2};
  PetscBool       isSimplice=PETSC_TRUE;
  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(DMPlexComputeGeometryFEM(dm, &cellgeom));
  vol = get_volume_DM_v2(dm,cellgeom);
  PetscPrintf(PETSC_COMM_WORLD, "Volume %f\n",vol);

  PetscCall(VecDestroy(&cellgeom));
  PetscCall(DMDestroy(&dm));

  PetscCall(PetscFinalize());
  return 0;
}

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

  DMGetDS(dm,&prob);
  PetscDSGetDiscretization(prob,0,(PetscObject*) &fe);
  PetscFEGetQuadrature(fe,&quad);
  PetscQuadratureView(quad,PETSC_VIEWER_STDOUT_WORLD);
  PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);

  DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
  VecGetArray(cellgeom, &cgeom);
  for (c = cStart; c < cEnd; ++c) {
    PetscFEGeom *cg;
    DMPlexPointLocalRef(dm, c, cgeom, &cg);
    for (q = 0; q < Nq; ++q) {
      if (cg->detJ < 0) {
        ierr = PetscPrintf(PETSC_COMM_WORLD, "ATTENTION, volume non correct !!\n\n");
        ierr = PetscFinalize();
        return -1;
      }
      vol += *cg->detJ*quadWeights[q*qNc];
    }
  }

  ierr = MPI_Allreduce(&vol, &volG, 1,MPIU_REAL, MPI_SUM,PETSC_COMM_WORLD );CHKERRQ(ierr);

  return volG;
}