[petsc-users] grid adaptivity with dmforest

[Preda Silvia via petsc-users]

Hi,

I'm using the toycode below to manage adaptivity using a DMForest, based on 
p4est.

The code simply performs three adaptivity steps on a square grid [0,1]x[0,1], 
uniformly refined at level 2 at the beginning. The minimum and the maximum 
level of refinement are set to 1 and 5, respectively. During the adaptivity 
procedure, a quadrant is refined if its centroid is inside the circle of radius 
0.25, centred in (0.5,0.5).

I'm facing two main issues:

  *   As far as I understood, when the label asks to refine a quadrant which is 
already at the maximum level, the adaption procedure errors out instead of 
ignoring the request. Thus I need to check the quadrant level before setting 
the adaptivity label. How can I get access to the quadrant-level of a cell?
  *   Once the mesh is adapted, data need to be projected on the new mesh and 
in my method I'll need to write an elaborate ad-hoc procedure for that. How can 
I access the map that provides the correspondence between the indexes of 
outcoming quadrants and the indexes of incoming ones?

Below is the code, to be run with these options:
-dm_type p4est
-dm_forest_topology brick
-dm_p4est_brick_size 1,1
-dm_view vtk:brick.vtu
-dm_forest_initial_refinement 2
-dm_forest_minimum_refinement 1
-dm_forest_maximum_refinement 5

static char help[] = "Create and view a forest mesh\n\n";

#include <petscdmforest.h>
#include <petscdmplex.h>
#include <petscoptions.h>

static PetscErrorCode CreateAdaptLabel(DM dm, DM dmConv, PetscInt *nAdaptLoc)
{
  DMLabel  adaptLabel;
  PetscInt cStart, cEnd, c;

  PetscFunctionBeginUser;
  PetscCall(DMGetCoordinatesLocalSetUp(dm));
  PetscCall(DMCreateLabel(dm, "adaptLabel"));
  PetscCall(DMGetLabel(dm, "adaptLabel", &adaptLabel));
  PetscCall(DMForestGetCellChart(dm, &cStart, &cEnd));
  for (c = cStart; c < cEnd; ++c) {
    PetscReal centroid[3], volume, x, y;
    PetscCall(DMPlexComputeCellGeometryFVM(dmConv, c, &volume, centroid, NULL));
    x = centroid[0];
    y = centroid[1];
    if (std::sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5))<0.25) {
      PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));
      ++nAdaptLoc[0];
    } else {
      PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_KEEP));
      ++nAdaptLoc[1];
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ForestToPlex(DM *dm, DM *dmConv)
{
  PetscFunctionBeginUser;
  PetscCall(DMConvert(*dm, DMPLEX, dmConv));
  PetscCall(DMLocalizeCoordinates(*dmConv));
  PetscCall(DMViewFromOptions(*dmConv, NULL, "-dm_conv_view"));
  PetscCall(DMPlexCheckCellShape(*dmConv, PETSC_FALSE, PETSC_DETERMINE));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode AdaptMesh(DM *dm)
{
  DM              dmCur = *dm;
  PetscBool       hasLabel=PETSC_FALSE, adapt=PETSC_TRUE;
  PetscInt        adaptIter=0, maxAdaptIter=3;

  PetscFunctionBeginUser;
  while (adapt) {
    DM       dmAdapt;
    DMLabel  adaptLabel;
    PetscInt nAdaptLoc[2]={0,0}, nAdapt[2]={0,0};

    ++adaptIter;
    PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPT ITER %d\n",adaptIter));

    DM dmConv;
    PetscCall(ForestToPlex(&dmCur,&dmConv));
    PetscCall(CreateAdaptLabel(dmCur,dmConv,nAdaptLoc));
    PetscCallMPI(MPIU_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, 
PetscObjectComm((PetscObject)dmCur)));
    PetscCall(DMGetLabel(dmCur, "adaptLabel", &adaptLabel));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cell to refine = %d\n",nAdapt[0]));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cell to keep = %d\n",nAdapt[1]));

    if (nAdapt[0]) {
      PetscCall(DMAdaptLabel(dmCur, adaptLabel, &dmAdapt));
      PetscCall(DMHasLabel(dmAdapt, "adaptLabel", &hasLabel));
      PetscCall(DMDestroy(&dmCur));
      PetscCall(DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view"));
      dmCur = dmAdapt;
    }
    //PetscCall(DMLabelDestroy(&adaptLabel));
    PetscCall(DMDestroy(&dmConv));
    if (adaptIter==maxAdaptIter) adapt=PETSC_FALSE;
  }
  *dm = dmCur;
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  DM          dm;
  char        typeString[256] = {'\0'};
  PetscViewer viewer          = NULL;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
  PetscCall(PetscStrncpy(typeString, DMFOREST, 256));
  PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "DM Forest example options", NULL);
  PetscCall(PetscOptionsString("-dm_type", "The type of the dm", NULL, 
DMFOREST, typeString, sizeof(typeString), NULL));
  PetscOptionsEnd();

  PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n ==== TOY CODE DMFOREST WITH AMR 
====\n"));

  PetscCall(DMSetType(dm, (DMType)typeString));
  PetscCall(DMSetFromOptions(dm));
  PetscCall(DMSetUp(dm));

  /* Adapt */
  PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPTIVITY PHASE STARTED\n"));
  PetscCall(AdaptMesh(&dm));
  PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPTIVITY PHASE ENDED\n\n"));

  PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
  PetscCall(PetscViewerDestroy(&viewer));

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


Thanks a lot,

Silvia