<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Aptos;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
font-size:11.0pt;
font-family:"Aptos",sans-serif;
mso-ligatures:standardcontextual;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0cm;
margin-right:0cm;
margin-bottom:0cm;
margin-left:36.0pt;
font-size:11.0pt;
font-family:"Aptos",sans-serif;
mso-ligatures:standardcontextual;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Aptos",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:11.0pt;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 2.0cm 2.0cm 2.0cm;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:1210610862;
mso-list-type:hybrid;
mso-list-template-ids:-1057463690 -1110409478 67698691 67698693 67698689 67698691 67698693 67698689 67698691 67698693;}
@list l0:level1
{mso-level-start-at:0;
mso-level-number-format:bullet;
mso-level-text:-;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Aptos",sans-serif;
mso-fareast-font-family:Aptos;
mso-bidi-font-family:Arial;}
@list l0:level2
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level3
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
@list l0:level4
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Symbol;}
@list l0:level5
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level6
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
@list l0:level7
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Symbol;}
@list l0:level8
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level9
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
ol
{margin-bottom:0cm;}
ul
{margin-bottom:0cm;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="#467886" vlink="#96607D" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span lang="IT">Hi,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="IT"><o:p> </o:p></span></p>
<p class="MsoNormal">I’m using the toycode below to manage adaptivity using a DMForest, based on p4est.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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).
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I’m facing two main issues:<o:p></o:p></p>
<ul style="margin-top:0cm" type="disc">
<li class="MsoListParagraph" style="margin-left:0cm;mso-list:l0 level1 lfo1">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?<o:p></o:p></li><li class="MsoListParagraph" style="margin-left:0cm;mso-list:l0 level1 lfo1">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?<o:p></o:p></li></ul>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Below is the code, to be run with these options: <o:p></o:p></p>
<p class="MsoNormal">-dm_type p4est <o:p></o:p></p>
<p class="MsoNormal">-dm_forest_topology brick <o:p></o:p></p>
<p class="MsoNormal">-dm_p4est_brick_size 1,1 <o:p></o:p></p>
<p class="MsoNormal">-dm_view vtk:brick.vtu<o:p></o:p></p>
<p class="MsoNormal">-dm_forest_initial_refinement 2 <o:p></o:p></p>
<p class="MsoNormal">-dm_forest_minimum_refinement 1 <o:p></o:p></p>
<p class="MsoNormal">-dm_forest_maximum_refinement 5<o:p></o:p></p>
<div style="mso-element:para-border-div;border:none;border-bottom:double windowtext 2.25pt;padding:0cm 0cm 1.0pt 0cm">
<p class="MsoNormal" style="border:none;padding:0cm"><o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div style="mso-element:para-border-div;border:none;border-bottom:double windowtext 2.25pt;padding:0cm 0cm 1.0pt 0cm">
<p class="MsoNormal" style="border:none;padding:0cm">static char help[] = "Create and view a forest mesh\n\n";<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">#include <petscdmforest.h><o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">#include <petscdmplex.h><o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">#include <petscoptions.h><o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">static PetscErrorCode CreateAdaptLabel(DM dm, DM dmConv, PetscInt *nAdaptLoc)
<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">{<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> DMLabel adaptLabel;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscInt cStart, cEnd, c;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscFunctionBeginUser;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMGetCoordinatesLocalSetUp(dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMCreateLabel(dm, "adaptLabel"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMGetLabel(dm, "adaptLabel", &adaptLabel));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMForestGetCellChart(dm, &cStart, &cEnd));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> for (c = cStart; c < cEnd; ++c) {<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscReal centroid[3], volume, x, y;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMPlexComputeCellGeometryFVM(dmConv, c, &volume, centroid, NULL));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <span lang="FR">x = centroid[0];<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> y = centroid[1];<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> if (std::sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5))<0.25) {<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> ++nAdaptLoc[0];<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> } else {<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_KEEP));<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> ++nAdaptLoc[1];<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> }<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> }<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscFunctionReturn(PETSC_SUCCESS);<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR">}<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"><o:p> </o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR">static PetscErrorCode ForestToPlex(DM *dm, DM *dmConv)<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR">{<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscFunctionBeginUser;<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscCall(DMConvert(*dm, DMPLEX, dmConv));<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> PetscCall(DMLocalizeCoordinates(*dmConv));<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="FR"> </span>PetscCall(DMViewFromOptions(*dmConv, NULL, "-dm_conv_view"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMPlexCheckCellShape(*dmConv, PETSC_FALSE, PETSC_DETERMINE));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscFunctionReturn(PETSC_SUCCESS);<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">}<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">static PetscErrorCode AdaptMesh(DM *dm)<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">{<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> DM dmCur = *dm;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscBool hasLabel=PETSC_FALSE, adapt=PETSC_TRUE;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscInt adaptIter=0, maxAdaptIter=3;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscFunctionBeginUser;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> while (adapt) {<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> DM dmAdapt;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> DMLabel adaptLabel;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscInt nAdaptLoc[2]={0,0}, nAdapt[2]={0,0};<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> ++adaptIter;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPT ITER %d\n",adaptIter));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> DM dmConv;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(ForestToPlex(&dmCur,&dmConv));
<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(CreateAdaptLabel(dmCur,dmConv,nAdaptLoc));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCallMPI(MPIU_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dmCur)));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMGetLabel(dmCur, "adaptLabel", &adaptLabel));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cell to refine = %d\n",nAdapt[0]));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cell to keep = %d\n",nAdapt[1]));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> if (nAdapt[0]) {<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMAdaptLabel(dmCur, adaptLabel, &dmAdapt));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMHasLabel(dmAdapt, "adaptLabel", &hasLabel));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMDestroy(&dmCur));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> dmCur = dmAdapt;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> }<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> //PetscCall(DMLabelDestroy(&adaptLabel));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMDestroy(&dmConv));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> if (adaptIter==maxAdaptIter) adapt=PETSC_FALSE;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> }<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> *dm = dmCur;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscFunctionReturn(PETSC_SUCCESS);<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE">}<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE"><o:p> </o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE">int main(int argc, char **argv)<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE">{<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE"> DM dm;<o:p></o:p></span></p>
<p class="MsoNormal" style="border:none;padding:0cm"><span lang="DE"> </span>char typeString[256] = {'\0'};<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscViewer viewer = NULL;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscFunctionBeginUser;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscInitialize(&argc, &argv, NULL, help));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscStrncpy(typeString, DMFOREST, 256));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "DM Forest example options", NULL);<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscOptionsString("-dm_type", "The type of the dm", NULL, DMFOREST, typeString, sizeof(typeString), NULL));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscOptionsEnd();<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n ==== TOY CODE DMFOREST WITH AMR ====\n"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMSetType(dm, (DMType)typeString));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMSetFromOptions(dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMSetUp(dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> /* Adapt */<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPTIVITY PHASE STARTED\n"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(AdaptMesh(&dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscPrintf(PETSC_COMM_SELF,"\nADAPTIVITY PHASE ENDED\n\n"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> <o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscViewerDestroy(&viewer));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(DMDestroy(&dm));<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> PetscCall(PetscFinalize());<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"> return 0;<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm">}<o:p></o:p></p>
<p class="MsoNormal" style="border:none;padding:0cm"><o:p> </o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks a lot,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Silvia<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</body>
</html>