On Wed, Feb 19, 2025 at 11:18 AM Preda Silvia via petsc-users <
petsc-users@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$
>
