[petsc-users] Refine DMPlex with a Refinement Function

[Daniel Mckinnell via petsc-users]

Hi,

I have been having some trouble trying to refine a DMPlex object using a 
Refinement Function. I have been working with reference to the code discussed 
here: 
https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2019-April/038341.html 
including the code from the GitHub directory mentioned. The error I have been 
getting is as follows:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: No grid refiner of dimension 2 registered
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.11.2, unknown
[0]PETSC ERROR: ./fpf on a arch-linux2-c-debug named cromars by daniel Thu Aug  
1 15:40:18 2019
[0]PETSC ERROR: Configure options −−download−f2cblaslapack=yes 
−−with−debugging=1 −−download−metis=yes −−download−parmetis=yes 
−−with−fortran−bindings=0 −−with−python=0 -download-openmpi=1 −−with−c−support 
−−with−clanguage=cxx
[0]PETSC ERROR: #1 DMPlexRefine_Internal() line 215 in 
/home/daniel/petsc/src/dm/impls/plex/plexadapt.c
[0]PETSC ERROR: #2 DMRefine_Plex() line 10381 in 
/home/daniel/petsc/src/dm/impls/plex/plexrefine.c
[0]PETSC ERROR: #3 DMRefine() line 1881 in 
/home/daniel/petsc/src/dm/interface/dm.c

Any help would be greatly appreciated.

Thanks,

Daniel
#include <petscdmplex.h>
#include <iostream>

PetscErrorCode testrefinement(const PetscReal xloc[], PetscReal *volumelimit)
{
  PetscErrorCode ierr;

  const PetscReal testcentroid[3] = {0.25, 0.25, 0.25}; 

  double coord[3]= {xloc[0],xloc[1],xloc[2]};
  //transform the point and return the intensity value
  *volumelimit = 0.125;

  bool refineelement = (PetscSqrtReal((coord[0] - testcentroid[0] )*(coord[0] - testcentroid[0] ) + (coord[1] - testcentroid[1] )*(coord[1] - testcentroid[1] )+(coord[2] - testcentroid[2] )*(coord[2] - testcentroid[2] ))  < 0.05)  ;

  // refine at the edges
  if (refineelement)
  {
    *volumelimit = 0.01;
  }
   
  PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
DM             dm, dmDist = NULL, refinedm;
PetscSF        sf;
Vec            u;
PetscInt       dim = 3, xcells = 2, ycells = 2, zcells = 2;
PetscBool      interpolate = PETSC_TRUE;
PetscMPIInt    rank;

try {
     PetscInitialize(&argc, &argv, NULL, NULL);
}
catch (const std::exception& e) {return(-1);};
PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);

DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, (&xcells,&ycells,&zcells), NULL, NULL, NULL, interpolate, &dm);
DMView(dm, PETSC_VIEWER_STDOUT_WORLD);


DMPlexDistribute(dm, 0, &sf, &dmDist);
if (dmDist) {DMDestroy(&dm); dm = dmDist;}
MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

DMPlexSetRefinementUniform(dm, PETSC_FALSE);
DMPlexSetRefinementFunction(dm, testrefinement);

DMRefine(dm, PETSC_COMM_WORLD, &refinedm);
  if (refinedm) {
    DMDestroy(&dm);
    dm  = refinedm;
}

DMGetCoordinates(dm, &u);
VecView(u, PETSC_VIEWER_STDOUT_WORLD);
DMSetCoordinates(dm, u);

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