[petsc-users] grid adaptivity with dmforest

[Matthew Knepley]

On Wed, Feb 19, 2025 at 11:18 AM Preda Silvia via petsc-users <
petsc-users at mcs.anl.gov> wrote:

> 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?
>
>
Toby understands this better than me. Toby, two questions:

  1) Can we catch that error and just return?

  2) Will DMProjectField() work in the case of multiple refinements like
this?

  Thanks,

      Matt


> 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
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ZSCq5o18nPkr_fesXZsTE1dDngSGlwHv1_gJXKWrkgdaUqTlCeNxjF-VSyRK8bjgaHXIu9LJw-fH5fP0nQ8Z$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ZSCq5o18nPkr_fesXZsTE1dDngSGlwHv1_gJXKWrkgdaUqTlCeNxjF-VSyRK8bjgaHXIu9LJw-fH5Wo8ZisR$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250219/7b933fcd/attachment.html>