Actual source code: pforest.h

  1: #include <petscds.h>
  2: #include <petsc/private/dmimpl.h>
  3: #include <petsc/private/dmforestimpl.h>
  4: #include <petsc/private/dmpleximpl.h>
  5: #include <petsc/private/dmlabelimpl.h>
  6: #include <petsc/private/viewerimpl.h>
  7: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
  8: #include "petsc_p4est_package.h"

 10: #if defined(PETSC_HAVE_P4EST)

 12:   #if !defined(P4_TO_P8)
 13:     #include <p4est.h>
 14:     #include <p4est_extended.h>
 15:     #include <p4est_geometry.h>
 16:     #include <p4est_ghost.h>
 17:     #include <p4est_lnodes.h>
 18:     #include <p4est_vtk.h>
 19:     #include <p4est_plex.h>
 20:     #include <p4est_bits.h>
 21:     #include <p4est_algorithms.h>
 22:   #else
 23:     #include <p8est.h>
 24:     #include <p8est_extended.h>
 25:     #include <p8est_geometry.h>
 26:     #include <p8est_ghost.h>
 27:     #include <p8est_lnodes.h>
 28:     #include <p8est_vtk.h>
 29:     #include <p8est_plex.h>
 30:     #include <p8est_bits.h>
 31:     #include <p8est_algorithms.h>
 32:   #endif

 34: typedef enum {
 35:   PATTERN_HASH,
 36:   PATTERN_FRACTAL,
 37:   PATTERN_CORNER,
 38:   PATTERN_CENTER,
 39:   PATTERN_COUNT
 40: } DMRefinePattern;
 41: static const char *DMRefinePatternName[PATTERN_COUNT] = {"hash", "fractal", "corner", "center"};

 43: typedef struct _DMRefinePatternCtx {
 44:   PetscInt       corner;
 45:   PetscBool      fractal[P4EST_CHILDREN];
 46:   PetscReal      hashLikelihood;
 47:   PetscInt       maxLevel;
 48:   p4est_refine_t refine_fn;
 49: } DMRefinePatternCtx;

 51: static int DMRefinePattern_Corner(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
 52: {
 53:   p4est_quadrant_t    root, rootcorner;
 54:   DMRefinePatternCtx *ctx;

 56:   ctx = (DMRefinePatternCtx *)p4est->user_pointer;
 57:   if (quadrant->level >= ctx->maxLevel) return 0;

 59:   root.x = root.y = 0;
 60:   #if defined(P4_TO_P8)
 61:   root.z = 0;
 62:   #endif
 63:   root.level = 0;
 64:   p4est_quadrant_corner_descendant(&root, &rootcorner, ctx->corner, quadrant->level);
 65:   if (p4est_quadrant_is_equal(quadrant, &rootcorner)) return 1;
 66:   return 0;
 67: }

 69: static int DMRefinePattern_Center(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
 70: {
 71:   int                 cid;
 72:   p4est_quadrant_t    ancestor, ancestorcorner;
 73:   DMRefinePatternCtx *ctx;

 75:   ctx = (DMRefinePatternCtx *)p4est->user_pointer;
 76:   if (quadrant->level >= ctx->maxLevel) return 0;
 77:   if (quadrant->level <= 1) return 1;

 79:   p4est_quadrant_ancestor(quadrant, 1, &ancestor);
 80:   cid = p4est_quadrant_child_id(&ancestor);
 81:   p4est_quadrant_corner_descendant(&ancestor, &ancestorcorner, P4EST_CHILDREN - 1 - cid, quadrant->level);
 82:   if (p4est_quadrant_is_equal(quadrant, &ancestorcorner)) return 1;
 83:   return 0;
 84: }

 86: static int DMRefinePattern_Fractal(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
 87: {
 88:   int                 cid;
 89:   DMRefinePatternCtx *ctx;

 91:   ctx = (DMRefinePatternCtx *)p4est->user_pointer;
 92:   if (quadrant->level >= ctx->maxLevel) return 0;
 93:   if (!quadrant->level) return 1;
 94:   cid = p4est_quadrant_child_id(quadrant);
 95:   if (ctx->fractal[cid ^ ((int)(quadrant->level % P4EST_CHILDREN))]) return 1;
 96:   return 0;
 97: }

 99:   /* simplified from MurmurHash3 by Austin Appleby */
100:   #define DMPROT32(x, y) ((x << y) | (x >> (32 - y)))
101: static uint32_t DMPforestHash(const uint32_t *blocks, uint32_t nblocks)
102: {
103:   uint32_t c1   = 0xcc9e2d51;
104:   uint32_t c2   = 0x1b873593;
105:   uint32_t r1   = 15;
106:   uint32_t r2   = 13;
107:   uint32_t m    = 5;
108:   uint32_t n    = 0xe6546b64;
109:   uint32_t hash = 0;
110:   int      len  = nblocks * 4;
111:   uint32_t i;

113:   for (i = 0; i < nblocks; i++) {
114:     uint32_t k;

116:     k = blocks[i];
117:     k *= c1;
118:     k = DMPROT32(k, r1);
119:     k *= c2;

121:     hash ^= k;
122:     hash = DMPROT32(hash, r2) * m + n;
123:   }

125:   hash ^= len;
126:   hash ^= (hash >> 16);
127:   hash *= 0x85ebca6b;
128:   hash ^= (hash >> 13);
129:   hash *= 0xc2b2ae35;
130:   hash ^= (hash >> 16);

132:   return hash;
133: }

135:   #if defined(UINT32_MAX)
136:     #define DMP4EST_HASH_MAX UINT32_MAX
137:   #else
138:     #define DMP4EST_HASH_MAX ((uint32_t)0xffffffff)
139:   #endif

141: static int DMRefinePattern_Hash(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
142: {
143:   uint32_t            data[5];
144:   uint32_t            result;
145:   DMRefinePatternCtx *ctx;

147:   ctx = (DMRefinePatternCtx *)p4est->user_pointer;
148:   if (quadrant->level >= ctx->maxLevel) return 0;
149:   data[0] = ((uint32_t)quadrant->level) << 24;
150:   data[1] = (uint32_t)which_tree;
151:   data[2] = (uint32_t)quadrant->x;
152:   data[3] = (uint32_t)quadrant->y;
153:   #if defined(P4_TO_P8)
154:   data[4] = (uint32_t)quadrant->z;
155:   #endif

157:   result = DMPforestHash(data, 2 + P4EST_DIM);
158:   if (((double)result / (double)DMP4EST_HASH_MAX) < ctx->hashLikelihood) return 1;
159:   return 0;
160: }

162:   #define DMConvert_pforest_plex _infix_pforest(DMConvert, _plex)
163: static PetscErrorCode DMConvert_pforest_plex(DM, DMType, DM *);

165:   #define DMFTopology_pforest _append_pforest(DMFTopology)
166: typedef struct {
167:   PetscInt              refct;
168:   p4est_connectivity_t *conn;
169:   p4est_geometry_t     *geom;
170:   PetscInt             *tree_face_to_uniq; /* p4est does not explicitly enumerate facets, but we must to keep track of labels */
171: } DMFTopology_pforest;

173:   #define DM_Forest_pforest _append_pforest(DM_Forest)
174: typedef struct {
175:   DMFTopology_pforest *topo;
176:   p4est_t             *forest;
177:   p4est_ghost_t       *ghost;
178:   p4est_lnodes_t      *lnodes;
179:   PetscBool            partition_for_coarsening;
180:   PetscBool            coarsen_hierarchy;
181:   PetscBool            labelsFinalized;
182:   PetscBool            adaptivitySuccess;
183:   PetscInt             cLocalStart;
184:   PetscInt             cLocalEnd;
185:   DM                   plex;
186:   char                *ghostName;
187:   PetscSF              pointAdaptToSelfSF;
188:   PetscSF              pointSelfToAdaptSF;
189:   PetscInt            *pointAdaptToSelfCids;
190:   PetscInt            *pointSelfToAdaptCids;
191: } DM_Forest_pforest;

193:   #define DM_Forest_geometry_pforest _append_pforest(DM_Forest_geometry)
194: typedef struct {
195:   DM base;
196:   PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
197:   void             *mapCtx;
198:   PetscInt          coordDim;
199:   p4est_geometry_t *inner;
200: } DM_Forest_geometry_pforest;

202:   #define GeometryMapping_pforest _append_pforest(GeometryMapping)
203: static void GeometryMapping_pforest(p4est_geometry_t *geom, p4est_topidx_t which_tree, const double abc[3], double xyz[3])
204: {
205:   DM_Forest_geometry_pforest *geom_pforest = (DM_Forest_geometry_pforest *)geom->user;
206:   PetscReal                   PetscABC[3]  = {0.};
207:   PetscReal                   PetscXYZ[3]  = {0.};
208:   PetscInt                    i, d = PetscMin(3, geom_pforest->coordDim);
209:   double                      ABC[3];
210:   PetscErrorCode              ierr;

212:   (geom_pforest->inner->X)(geom_pforest->inner, which_tree, abc, ABC);

214:   for (i = 0; i < d; i++) PetscABC[i] = ABC[i];
215:   (geom_pforest->map)(geom_pforest->base, (PetscInt)which_tree, geom_pforest->coordDim, PetscABC, PetscXYZ, geom_pforest->mapCtx);
216:   PETSC_P4EST_ASSERT(!ierr);
217:   for (i = 0; i < d; i++) xyz[i] = PetscXYZ[i];
218: }

220:   #define GeometryDestroy_pforest _append_pforest(GeometryDestroy)
221: static void GeometryDestroy_pforest(p4est_geometry_t *geom)
222: {
223:   DM_Forest_geometry_pforest *geom_pforest = (DM_Forest_geometry_pforest *)geom->user;
224:   PetscErrorCode              ierr;

226:   p4est_geometry_destroy(geom_pforest->inner);
227:   PetscFree(geom->user);
228:   PETSC_P4EST_ASSERT(!ierr);
229:   PetscFree(geom);
230:   PETSC_P4EST_ASSERT(!ierr);
231: }

233:   #define DMFTopologyDestroy_pforest _append_pforest(DMFTopologyDestroy)
234: static PetscErrorCode DMFTopologyDestroy_pforest(DMFTopology_pforest **topo)
235: {
236:   if (!(*topo)) return 0;
237:   if (--((*topo)->refct) > 0) {
238:     *topo = NULL;
239:     return 0;
240:   }
241:   if ((*topo)->geom) PetscCallP4est(p4est_geometry_destroy, ((*topo)->geom));
242:   PetscCallP4est(p4est_connectivity_destroy, ((*topo)->conn));
243:   PetscFree((*topo)->tree_face_to_uniq);
244:   PetscFree(*topo);
245:   *topo = NULL;
246:   return 0;
247: }

249: static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *, PetscInt **);

251:   #define DMFTopologyCreateBrick_pforest _append_pforest(DMFTopologyCreateBrick)
252: static PetscErrorCode DMFTopologyCreateBrick_pforest(DM dm, PetscInt N[], PetscInt P[], PetscReal B[], DMFTopology_pforest **topo, PetscBool useMorton)
253: {
254:   double  *vertices;
255:   PetscInt i, numVerts;

258:   PetscNew(topo);

260:   (*topo)->refct = 1;
261:   #if !defined(P4_TO_P8)
262:   PetscCallP4estReturn((*topo)->conn, p4est_connectivity_new_brick, ((int)N[0], (int)N[1], (P[0] == DM_BOUNDARY_NONE) ? 0 : 1, (P[1] == DM_BOUNDARY_NONE) ? 0 : 1));
263:   #else
264:   PetscCallP4estReturn((*topo)->conn, p8est_connectivity_new_brick, ((int)N[0], (int)N[1], (int)N[2], (P[0] == DM_BOUNDARY_NONE) ? 0 : 1, (P[1] == DM_BOUNDARY_NONE) ? 0 : 1, (P[2] == DM_BOUNDARY_NONE) ? 0 : 1));
265:   #endif
266:   numVerts = (*topo)->conn->num_vertices;
267:   vertices = (*topo)->conn->vertices;
268:   for (i = 0; i < 3 * numVerts; i++) {
269:     PetscInt j = i % 3;

271:     vertices[i] = B[2 * j] + (vertices[i] / N[j]) * (B[2 * j + 1] - B[2 * j]);
272:   }
273:   (*topo)->geom = NULL;
274:   PforestConnectivityEnumerateFacets((*topo)->conn, &(*topo)->tree_face_to_uniq);
275:   return 0;
276: }

278:   #define DMFTopologyCreate_pforest _append_pforest(DMFTopologyCreate)
279: static PetscErrorCode DMFTopologyCreate_pforest(DM dm, DMForestTopology topologyName, DMFTopology_pforest **topo)
280: {
281:   const char *name = (const char *)topologyName;
282:   const char *prefix;
283:   PetscBool   isBrick, isShell, isSphere, isMoebius;

288:   PetscStrcmp(name, "brick", &isBrick);
289:   PetscStrcmp(name, "shell", &isShell);
290:   PetscStrcmp(name, "sphere", &isSphere);
291:   PetscStrcmp(name, "moebius", &isMoebius);
292:   PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
293:   if (isBrick) {
294:     PetscBool flgN, flgP, flgM, flgB, useMorton = PETSC_TRUE, periodic = PETSC_FALSE;
295:     PetscInt  N[3] = {2, 2, 2}, P[3] = {0, 0, 0}, nretN = P4EST_DIM, nretP = P4EST_DIM, nretB = 2 * P4EST_DIM, i;
296:     PetscReal B[6] = {0.0, 1.0, 0.0, 1.0, 0.0, 1.0}, Lstart[3] = {0., 0., 0.}, L[3] = {-1.0, -1.0, -1.0}, maxCell[3] = {-1.0, -1.0, -1.0};

298:     if (dm->setfromoptionscalled) {
299:       PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_size", N, &nretN, &flgN);
300:       PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_periodicity", P, &nretP, &flgP);
301:       PetscOptionsGetRealArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_bounds", B, &nretB, &flgB);
302:       PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_use_morton_curve", &useMorton, &flgM);
306:     }
307:     for (i = 0; i < P4EST_DIM; i++) {
308:       P[i]     = (P[i] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE);
309:       periodic = (PetscBool)(P[i] || periodic);
310:       if (!flgB) B[2 * i + 1] = N[i];
311:       if (P[i]) {
312:         Lstart[i]  = B[2 * i + 0];
313:         L[i]       = B[2 * i + 1] - B[2 * i + 0];
314:         maxCell[i] = 1.1 * (L[i] / N[i]);
315:       }
316:     }
317:     DMFTopologyCreateBrick_pforest(dm, N, P, B, topo, useMorton);
318:     if (periodic) DMSetPeriodicity(dm, maxCell, Lstart, L);
319:   } else {
320:     PetscNew(topo);

322:     (*topo)->refct = 1;
323:     PetscCallP4estReturn((*topo)->conn, p4est_connectivity_new_byname, (name));
324:     (*topo)->geom = NULL;
325:     if (isMoebius) DMSetCoordinateDim(dm, 3);
326:   #if defined(P4_TO_P8)
327:     if (isShell) {
328:       PetscReal R2 = 1., R1 = .55;

330:       if (dm->setfromoptionscalled) {
331:         PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_shell_outer_radius", &R2, NULL);
332:         PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_shell_inner_radius", &R1, NULL);
333:       }
334:       PetscCallP4estReturn((*topo)->geom, p8est_geometry_new_shell, ((*topo)->conn, R2, R1));
335:     } else if (isSphere) {
336:       PetscReal R2 = 1., R1 = 0.191728, R0 = 0.039856;

338:       if (dm->setfromoptionscalled) {
339:         PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_outer_radius", &R2, NULL);
340:         PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_inner_radius", &R1, NULL);
341:         PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_core_radius", &R0, NULL);
342:       }
343:       PetscCallP4estReturn((*topo)->geom, p8est_geometry_new_sphere, ((*topo)->conn, R2, R1, R0));
344:     }
345:   #endif
346:     PforestConnectivityEnumerateFacets((*topo)->conn, &(*topo)->tree_face_to_uniq);
347:   }
348:   return 0;
349: }

351:   #define DMConvert_plex_pforest _append_pforest(DMConvert_plex)
352: static PetscErrorCode DMConvert_plex_pforest(DM dm, DMType newtype, DM *pforest)
353: {
354:   MPI_Comm  comm;
355:   PetscBool isPlex;
356:   PetscInt  dim;
357:   void     *ctx;


361:   comm = PetscObjectComm((PetscObject)dm);
362:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
364:   DMGetDimension(dm, &dim);
366:   DMCreate(comm, pforest);
367:   DMSetType(*pforest, DMPFOREST);
368:   DMForestSetBaseDM(*pforest, dm);
369:   DMGetApplicationContext(dm, &ctx);
370:   DMSetApplicationContext(*pforest, ctx);
371:   DMCopyDisc(dm, *pforest);
372:   return 0;
373: }

375:   #define DMForestDestroy_pforest _append_pforest(DMForestDestroy)
376: static PetscErrorCode DMForestDestroy_pforest(DM dm)
377: {
378:   DM_Forest         *forest  = (DM_Forest *)dm->data;
379:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;

382:   if (pforest->lnodes) PetscCallP4est(p4est_lnodes_destroy, (pforest->lnodes));
383:   pforest->lnodes = NULL;
384:   if (pforest->ghost) PetscCallP4est(p4est_ghost_destroy, (pforest->ghost));
385:   pforest->ghost = NULL;
386:   if (pforest->forest) PetscCallP4est(p4est_destroy, (pforest->forest));
387:   pforest->forest = NULL;
388:   DMFTopologyDestroy_pforest(&pforest->topo);
389:   PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", NULL);
390:   PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", NULL);
391:   PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL);
392:   PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", NULL);
393:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL);
394:   PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", NULL);
395:   PetscFree(pforest->ghostName);
396:   DMDestroy(&pforest->plex);
397:   PetscSFDestroy(&pforest->pointAdaptToSelfSF);
398:   PetscSFDestroy(&pforest->pointSelfToAdaptSF);
399:   PetscFree(pforest->pointAdaptToSelfCids);
400:   PetscFree(pforest->pointSelfToAdaptCids);
401:   PetscFree(forest->data);
402:   return 0;
403: }

405:   #define DMForestTemplate_pforest _append_pforest(DMForestTemplate)
406: static PetscErrorCode DMForestTemplate_pforest(DM dm, DM tdm)
407: {
408:   DM_Forest_pforest *pforest  = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
409:   DM_Forest_pforest *tpforest = (DM_Forest_pforest *)((DM_Forest *)tdm->data)->data;

411:   if (pforest->topo) pforest->topo->refct++;
412:   DMFTopologyDestroy_pforest(&(tpforest->topo));
413:   tpforest->topo = pforest->topo;
414:   return 0;
415: }

417:   #define DMPlexCreateConnectivity_pforest _append_pforest(DMPlexCreateConnectivity)
418: static PetscErrorCode DMPlexCreateConnectivity_pforest(DM, p4est_connectivity_t **, PetscInt **);

420: typedef struct _PforestAdaptCtx {
421:   PetscInt  maxLevel;
422:   PetscInt  minLevel;
423:   PetscInt  currLevel;
424:   PetscBool anyChange;
425: } PforestAdaptCtx;

427: static int pforest_coarsen_currlevel(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
428: {
429:   PforestAdaptCtx *ctx       = (PforestAdaptCtx *)p4est->user_pointer;
430:   PetscInt         minLevel  = ctx->minLevel;
431:   PetscInt         currLevel = ctx->currLevel;

433:   if (quadrants[0]->level <= minLevel) return 0;
434:   return (int)((PetscInt)quadrants[0]->level == currLevel);
435: }

437: static int pforest_coarsen_uniform(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
438: {
439:   PforestAdaptCtx *ctx      = (PforestAdaptCtx *)p4est->user_pointer;
440:   PetscInt         minLevel = ctx->minLevel;

442:   return (int)((PetscInt)quadrants[0]->level > minLevel);
443: }

445: static int pforest_coarsen_flag_any(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
446: {
447:   PetscInt         i;
448:   PetscBool        any      = PETSC_FALSE;
449:   PforestAdaptCtx *ctx      = (PforestAdaptCtx *)p4est->user_pointer;
450:   PetscInt         minLevel = ctx->minLevel;

452:   if (quadrants[0]->level <= minLevel) return 0;
453:   for (i = 0; i < P4EST_CHILDREN; i++) {
454:     if (quadrants[i]->p.user_int == DM_ADAPT_KEEP) {
455:       any = PETSC_FALSE;
456:       break;
457:     }
458:     if (quadrants[i]->p.user_int == DM_ADAPT_COARSEN) {
459:       any = PETSC_TRUE;
460:       break;
461:     }
462:   }
463:   return any ? 1 : 0;
464: }

466: static int pforest_coarsen_flag_all(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
467: {
468:   PetscInt         i;
469:   PetscBool        all      = PETSC_TRUE;
470:   PforestAdaptCtx *ctx      = (PforestAdaptCtx *)p4est->user_pointer;
471:   PetscInt         minLevel = ctx->minLevel;

473:   if (quadrants[0]->level <= minLevel) return 0;
474:   for (i = 0; i < P4EST_CHILDREN; i++) {
475:     if (quadrants[i]->p.user_int != DM_ADAPT_COARSEN) {
476:       all = PETSC_FALSE;
477:       break;
478:     }
479:   }
480:   return all ? 1 : 0;
481: }

483: static void pforest_init_determine(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
484: {
485:   quadrant->p.user_int = DM_ADAPT_DETERMINE;
486: }

488: static int pforest_refine_uniform(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
489: {
490:   PforestAdaptCtx *ctx      = (PforestAdaptCtx *)p4est->user_pointer;
491:   PetscInt         maxLevel = ctx->maxLevel;

493:   return ((PetscInt)quadrant->level < maxLevel);
494: }

496: static int pforest_refine_flag(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
497: {
498:   PforestAdaptCtx *ctx      = (PforestAdaptCtx *)p4est->user_pointer;
499:   PetscInt         maxLevel = ctx->maxLevel;

501:   if ((PetscInt)quadrant->level >= maxLevel) return 0;

503:   return (quadrant->p.user_int == DM_ADAPT_REFINE);
504: }

506: static PetscErrorCode DMPforestComputeLocalCellTransferSF_loop(p4est_t *p4estFrom, PetscInt FromOffset, p4est_t *p4estTo, PetscInt ToOffset, p4est_topidx_t flt, p4est_topidx_t llt, PetscInt *toFineLeavesCount, PetscInt *toLeaves, PetscSFNode *fromRoots, PetscInt *fromFineLeavesCount, PetscInt *fromLeaves, PetscSFNode *toRoots)
507: {
508:   PetscMPIInt    rank = p4estFrom->mpirank;
509:   p4est_topidx_t t;
510:   PetscInt       toFineLeaves = 0, fromFineLeaves = 0;

512:   for (t = flt; t <= llt; t++) { /* count roots and leaves */
513:     p4est_tree_t     *treeFrom  = &(((p4est_tree_t *)p4estFrom->trees->array)[t]);
514:     p4est_tree_t     *treeTo    = &(((p4est_tree_t *)p4estTo->trees->array)[t]);
515:     p4est_quadrant_t *firstFrom = &treeFrom->first_desc;
516:     p4est_quadrant_t *firstTo   = &treeTo->first_desc;
517:     PetscInt          numFrom   = (PetscInt)treeFrom->quadrants.elem_count;
518:     PetscInt          numTo     = (PetscInt)treeTo->quadrants.elem_count;
519:     p4est_quadrant_t *quadsFrom = (p4est_quadrant_t *)treeFrom->quadrants.array;
520:     p4est_quadrant_t *quadsTo   = (p4est_quadrant_t *)treeTo->quadrants.array;
521:     PetscInt          currentFrom, currentTo;
522:     PetscInt          treeOffsetFrom = (PetscInt)treeFrom->quadrants_offset;
523:     PetscInt          treeOffsetTo   = (PetscInt)treeTo->quadrants_offset;
524:     int               comp;

526:     PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (firstFrom, firstTo));

529:     for (currentFrom = 0, currentTo = 0; currentFrom < numFrom && currentTo < numTo;) {
530:       p4est_quadrant_t *quadFrom = &quadsFrom[currentFrom];
531:       p4est_quadrant_t *quadTo   = &quadsTo[currentTo];

533:       if (quadFrom->level == quadTo->level) {
534:         if (toLeaves) {
535:           toLeaves[toFineLeaves]        = currentTo + treeOffsetTo + ToOffset;
536:           fromRoots[toFineLeaves].rank  = rank;
537:           fromRoots[toFineLeaves].index = currentFrom + treeOffsetFrom + FromOffset;
538:         }
539:         toFineLeaves++;
540:         currentFrom++;
541:         currentTo++;
542:       } else {
543:         int fromIsAncestor;

545:         PetscCallP4estReturn(fromIsAncestor, p4est_quadrant_is_ancestor, (quadFrom, quadTo));
546:         if (fromIsAncestor) {
547:           p4est_quadrant_t lastDesc;

549:           if (toLeaves) {
550:             toLeaves[toFineLeaves]        = currentTo + treeOffsetTo + ToOffset;
551:             fromRoots[toFineLeaves].rank  = rank;
552:             fromRoots[toFineLeaves].index = currentFrom + treeOffsetFrom + FromOffset;
553:           }
554:           toFineLeaves++;
555:           currentTo++;
556:           PetscCallP4est(p4est_quadrant_last_descendant, (quadFrom, &lastDesc, quadTo->level));
557:           PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (quadTo, &lastDesc));
558:           if (comp) currentFrom++;
559:         } else {
560:           p4est_quadrant_t lastDesc;

562:           if (fromLeaves) {
563:             fromLeaves[fromFineLeaves]    = currentFrom + treeOffsetFrom + FromOffset;
564:             toRoots[fromFineLeaves].rank  = rank;
565:             toRoots[fromFineLeaves].index = currentTo + treeOffsetTo + ToOffset;
566:           }
567:           fromFineLeaves++;
568:           currentFrom++;
569:           PetscCallP4est(p4est_quadrant_last_descendant, (quadTo, &lastDesc, quadFrom->level));
570:           PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (quadFrom, &lastDesc));
571:           if (comp) currentTo++;
572:         }
573:       }
574:     }
575:   }
576:   *toFineLeavesCount   = toFineLeaves;
577:   *fromFineLeavesCount = fromFineLeaves;
578:   return 0;
579: }

581: /* Compute the maximum level across all the trees */
582: static PetscErrorCode DMPforestGetRefinementLevel(DM dm, PetscInt *lev)
583: {
584:   p4est_topidx_t     t, flt, llt;
585:   DM_Forest         *forest      = (DM_Forest *)dm->data;
586:   DM_Forest_pforest *pforest     = (DM_Forest_pforest *)forest->data;
587:   PetscInt           maxlevelloc = 0;
588:   p4est_t           *p4est;

592:   p4est = pforest->forest;
593:   flt   = p4est->first_local_tree;
594:   llt   = p4est->last_local_tree;
595:   for (t = flt; t <= llt; t++) {
596:     p4est_tree_t *tree = &(((p4est_tree_t *)p4est->trees->array)[t]);
597:     maxlevelloc        = PetscMax((PetscInt)tree->maxlevel, maxlevelloc);
598:   }
599:   MPIU_Allreduce(&maxlevelloc, lev, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
600:   return 0;
601: }

603: /* Puts identity in coarseToFine */
604: /* assumes a matching partition */
605: static PetscErrorCode DMPforestComputeLocalCellTransferSF(MPI_Comm comm, p4est_t *p4estFrom, PetscInt FromOffset, p4est_t *p4estTo, PetscInt ToOffset, PetscSF *fromCoarseToFine, PetscSF *toCoarseFromFine)
606: {
607:   p4est_topidx_t flt, llt;
608:   PetscSF        fromCoarse, toCoarse;
609:   PetscInt       numRootsFrom, numRootsTo, numLeavesFrom, numLeavesTo;
610:   PetscInt      *fromLeaves = NULL, *toLeaves = NULL;
611:   PetscSFNode   *fromRoots = NULL, *toRoots = NULL;

613:   flt = p4estFrom->first_local_tree;
614:   llt = p4estFrom->last_local_tree;
615:   PetscSFCreate(comm, &fromCoarse);
616:   if (toCoarseFromFine) PetscSFCreate(comm, &toCoarse);
617:   numRootsFrom = p4estFrom->local_num_quadrants + FromOffset;
618:   numRootsTo   = p4estTo->local_num_quadrants + ToOffset;
619:   DMPforestComputeLocalCellTransferSF_loop(p4estFrom, FromOffset, p4estTo, ToOffset, flt, llt, &numLeavesTo, NULL, NULL, &numLeavesFrom, NULL, NULL);
620:   PetscMalloc1(numLeavesTo, &toLeaves);
621:   PetscMalloc1(numLeavesTo, &fromRoots);
622:   if (toCoarseFromFine) {
623:     PetscMalloc1(numLeavesFrom, &fromLeaves);
624:     PetscMalloc1(numLeavesFrom, &fromRoots);
625:   }
626:   DMPforestComputeLocalCellTransferSF_loop(p4estFrom, FromOffset, p4estTo, ToOffset, flt, llt, &numLeavesTo, toLeaves, fromRoots, &numLeavesFrom, fromLeaves, toRoots);
627:   if (!ToOffset && (numLeavesTo == numRootsTo)) { /* compress */
628:     PetscFree(toLeaves);
629:     PetscSFSetGraph(fromCoarse, numRootsFrom, numLeavesTo, NULL, PETSC_OWN_POINTER, fromRoots, PETSC_OWN_POINTER);
630:   } else PetscSFSetGraph(fromCoarse, numRootsFrom, numLeavesTo, toLeaves, PETSC_OWN_POINTER, fromRoots, PETSC_OWN_POINTER);
631:   *fromCoarseToFine = fromCoarse;
632:   if (toCoarseFromFine) {
633:     PetscSFSetGraph(toCoarse, numRootsTo, numLeavesFrom, fromLeaves, PETSC_OWN_POINTER, toRoots, PETSC_OWN_POINTER);
634:     *toCoarseFromFine = toCoarse;
635:   }
636:   return 0;
637: }

639: /* range of processes whose B sections overlap this ranks A section */
640: static PetscErrorCode DMPforestComputeOverlappingRanks(PetscMPIInt size, PetscMPIInt rank, p4est_t *p4estA, p4est_t *p4estB, PetscInt *startB, PetscInt *endB)
641: {
642:   p4est_quadrant_t *myCoarseStart = &(p4estA->global_first_position[rank]);
643:   p4est_quadrant_t *myCoarseEnd   = &(p4estA->global_first_position[rank + 1]);
644:   p4est_quadrant_t *globalFirstB  = p4estB->global_first_position;

646:   *startB = -1;
647:   *endB   = -1;
648:   if (p4estA->local_num_quadrants) {
649:     PetscInt lo, hi, guess;
650:     /* binary search to find interval containing myCoarseStart */
651:     lo    = 0;
652:     hi    = size;
653:     guess = rank;
654:     while (1) {
655:       int startCompMy, myCompEnd;

657:       PetscCallP4estReturn(startCompMy, p4est_quadrant_compare_piggy, (&globalFirstB[guess], myCoarseStart));
658:       PetscCallP4estReturn(myCompEnd, p4est_quadrant_compare_piggy, (myCoarseStart, &globalFirstB[guess + 1]));
659:       if (startCompMy <= 0 && myCompEnd < 0) {
660:         *startB = guess;
661:         break;
662:       } else if (startCompMy > 0) { /* guess is to high */
663:         hi = guess;
664:       } else { /* guess is to low */
665:         lo = guess + 1;
666:       }
667:       guess = lo + (hi - lo) / 2;
668:     }
669:     /* reset bounds, but not guess */
670:     lo = 0;
671:     hi = size;
672:     while (1) {
673:       int startCompMy, myCompEnd;

675:       PetscCallP4estReturn(startCompMy, p4est_quadrant_compare_piggy, (&globalFirstB[guess], myCoarseEnd));
676:       PetscCallP4estReturn(myCompEnd, p4est_quadrant_compare_piggy, (myCoarseEnd, &globalFirstB[guess + 1]));
677:       if (startCompMy < 0 && myCompEnd <= 0) { /* notice that the comparison operators are different from above */
678:         *endB = guess + 1;
679:         break;
680:       } else if (startCompMy >= 0) { /* guess is to high */
681:         hi = guess;
682:       } else { /* guess is to low */
683:         lo = guess + 1;
684:       }
685:       guess = lo + (hi - lo) / 2;
686:     }
687:   }
688:   return 0;
689: }

691: static PetscErrorCode DMPforestGetPlex(DM, DM *);

693:   #define DMSetUp_pforest _append_pforest(DMSetUp)
694: static PetscErrorCode DMSetUp_pforest(DM dm)
695: {
696:   DM_Forest         *forest  = (DM_Forest *)dm->data;
697:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
698:   DM                 base, adaptFrom;
699:   DMForestTopology   topoName;
700:   PetscSF            preCoarseToFine = NULL, coarseToPreFine = NULL;
701:   PforestAdaptCtx    ctx;

703:   ctx.minLevel  = PETSC_MAX_INT;
704:   ctx.maxLevel  = 0;
705:   ctx.currLevel = 0;
706:   ctx.anyChange = PETSC_FALSE;
707:   /* sanity check */
708:   DMForestGetAdaptivityForest(dm, &adaptFrom);
709:   DMForestGetBaseDM(dm, &base);
710:   DMForestGetTopology(dm, &topoName);

713:   /* === Step 1: DMFTopology === */
714:   if (adaptFrom) { /* reference already created topology */
715:     PetscBool          ispforest;
716:     DM_Forest         *aforest  = (DM_Forest *)adaptFrom->data;
717:     DM_Forest_pforest *apforest = (DM_Forest_pforest *)aforest->data;

719:     PetscObjectTypeCompare((PetscObject)adaptFrom, DMPFOREST, &ispforest);
722:     DMSetUp(adaptFrom);
723:     DMForestGetBaseDM(dm, &base);
724:     DMForestGetTopology(dm, &topoName);
725:   } else if (base) { /* construct a connectivity from base */
726:     PetscBool isPlex, isDA;

728:     PetscObjectGetName((PetscObject)base, &topoName);
729:     DMForestSetTopology(dm, topoName);
730:     PetscObjectTypeCompare((PetscObject)base, DMPLEX, &isPlex);
731:     PetscObjectTypeCompare((PetscObject)base, DMDA, &isDA);
732:     if (isPlex) {
733:       MPI_Comm              comm = PetscObjectComm((PetscObject)dm);
734:       PetscInt              depth;
735:       PetscMPIInt           size;
736:       p4est_connectivity_t *conn = NULL;
737:       DMFTopology_pforest  *topo;
738:       PetscInt             *tree_face_to_uniq = NULL;

740:       DMPlexGetDepth(base, &depth);
741:       if (depth == 1) {
742:         DM connDM;

744:         DMPlexInterpolate(base, &connDM);
745:         base = connDM;
746:         DMForestSetBaseDM(dm, base);
747:         DMDestroy(&connDM);
749:       MPI_Comm_size(comm, &size);
750:       if (size > 1) {
751:         DM      dmRedundant;
752:         PetscSF sf;

754:         DMPlexGetRedundantDM(base, &sf, &dmRedundant);
756:         PetscObjectCompose((PetscObject)dmRedundant, "_base_migration_sf", (PetscObject)sf);
757:         PetscSFDestroy(&sf);
758:         base = dmRedundant;
759:         DMForestSetBaseDM(dm, base);
760:         DMDestroy(&dmRedundant);
761:       }
762:       DMViewFromOptions(base, NULL, "-dm_p4est_base_view");
763:       DMPlexCreateConnectivity_pforest(base, &conn, &tree_face_to_uniq);
764:       PetscNew(&topo);
765:       topo->refct = 1;
766:       topo->conn  = conn;
767:       topo->geom  = NULL;
768:       {
769:         PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
770:         void *mapCtx;

772:         DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx);
773:         if (map) {
774:           DM_Forest_geometry_pforest *geom_pforest;
775:           p4est_geometry_t           *geom;

777:           PetscNew(&geom_pforest);
778:           DMGetCoordinateDim(dm, &geom_pforest->coordDim);
779:           geom_pforest->map    = map;
780:           geom_pforest->mapCtx = mapCtx;
781:           PetscCallP4estReturn(geom_pforest->inner, p4est_geometry_new_connectivity, (conn));
782:           PetscNew(&geom);
783:           geom->name    = topoName;
784:           geom->user    = geom_pforest;
785:           geom->X       = GeometryMapping_pforest;
786:           geom->destroy = GeometryDestroy_pforest;
787:           topo->geom    = geom;
788:         }
789:       }
790:       topo->tree_face_to_uniq = tree_face_to_uniq;
791:       pforest->topo           = topo;
793:   #if 0
794:       PetscInt N[3], P[3];

796:       /* get the sizes, periodicities */
797:       /* ... */
798:                                                                   /* don't use Morton order */
799:       DMFTopologyCreateBrick_pforest(dm,N,P,&pforest->topo,PETSC_FALSE);
800:   #endif
801:     {
802:       PetscInt numLabels, l;

804:       DMGetNumLabels(base, &numLabels);
805:       for (l = 0; l < numLabels; l++) {
806:         PetscBool   isDepth, isGhost, isVTK, isDim, isCellType;
807:         DMLabel     label, labelNew;
808:         PetscInt    defVal;
809:         const char *name;

811:         DMGetLabelName(base, l, &name);
812:         DMGetLabelByNum(base, l, &label);
813:         PetscStrcmp(name, "depth", &isDepth);
814:         if (isDepth) continue;
815:         PetscStrcmp(name, "dim", &isDim);
816:         if (isDim) continue;
817:         PetscStrcmp(name, "celltype", &isCellType);
818:         if (isCellType) continue;
819:         PetscStrcmp(name, "ghost", &isGhost);
820:         if (isGhost) continue;
821:         PetscStrcmp(name, "vtk", &isVTK);
822:         if (isVTK) continue;
823:         DMCreateLabel(dm, name);
824:         DMGetLabel(dm, name, &labelNew);
825:         DMLabelGetDefaultValue(label, &defVal);
826:         DMLabelSetDefaultValue(labelNew, defVal);
827:       }
828:       /* map dm points (internal plex) to base
829:          we currently create the subpoint_map for the entire hierarchy, starting from the finest forest
830:          and propagating back to the coarsest
831:          This is not an optimal approach, since we need the map only on the coarsest level
832:          during DMForestTransferVecFromBase */
833:       DMForestGetMinimumRefinement(dm, &l);
834:       if (!l) DMCreateLabel(dm, "_forest_base_subpoint_map");
835:     }
836:   } else { /* construct from topology name */
837:     DMFTopology_pforest *topo;

839:     DMFTopologyCreate_pforest(dm, topoName, &topo);
840:     pforest->topo = topo;
841:     /* TODO: construct base? */
842:   }

844:   /* === Step 2: get the leaves of the forest === */
845:   if (adaptFrom) { /* start with the old forest */
846:     DMLabel            adaptLabel;
847:     PetscInt           defaultValue;
848:     PetscInt           numValues, numValuesGlobal, cLocalStart, count;
849:     DM_Forest         *aforest  = (DM_Forest *)adaptFrom->data;
850:     DM_Forest_pforest *apforest = (DM_Forest_pforest *)aforest->data;
851:     PetscBool          computeAdaptSF;
852:     p4est_topidx_t     flt, llt, t;

854:     flt         = apforest->forest->first_local_tree;
855:     llt         = apforest->forest->last_local_tree;
856:     cLocalStart = apforest->cLocalStart;
857:     DMForestGetComputeAdaptivitySF(dm, &computeAdaptSF);
858:     PetscCallP4estReturn(pforest->forest, p4est_copy, (apforest->forest, 0)); /* 0 indicates no data copying */
859:     DMForestGetAdaptivityLabel(dm, &adaptLabel);
860:     if (adaptLabel) {
861:       /* apply the refinement/coarsening by flags, plus minimum/maximum refinement */
862:       DMLabelGetNumValues(adaptLabel, &numValues);
863:       MPI_Allreduce(&numValues, &numValuesGlobal, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)adaptFrom));
864:       DMLabelGetDefaultValue(adaptLabel, &defaultValue);
865:       if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN_LAST) { /* uniform coarsen of the last level only (equivalent to DM_ADAPT_COARSEN for conforming grids)  */
866:         DMForestGetMinimumRefinement(dm, &ctx.minLevel);
867:         DMPforestGetRefinementLevel(dm, &ctx.currLevel);
868:         pforest->forest->user_pointer = (void *)&ctx;
869:         PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_currlevel, NULL));
870:         pforest->forest->user_pointer = (void *)dm;
871:         PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
872:         /* we will have to change the offset after we compute the overlap */
873:         if (computeAdaptSF) DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL);
874:       } else if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN) { /* uniform coarsen */
875:         DMForestGetMinimumRefinement(dm, &ctx.minLevel);
876:         pforest->forest->user_pointer = (void *)&ctx;
877:         PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_uniform, NULL));
878:         pforest->forest->user_pointer = (void *)dm;
879:         PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
880:         /* we will have to change the offset after we compute the overlap */
881:         if (computeAdaptSF) DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL);
882:       } else if (!numValuesGlobal && defaultValue == DM_ADAPT_REFINE) { /* uniform refine */
883:         DMForestGetMaximumRefinement(dm, &ctx.maxLevel);
884:         pforest->forest->user_pointer = (void *)&ctx;
885:         PetscCallP4est(p4est_refine, (pforest->forest, 0, pforest_refine_uniform, NULL));
886:         pforest->forest->user_pointer = (void *)dm;
887:         PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
888:         /* we will have to change the offset after we compute the overlap */
889:         if (computeAdaptSF) DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), apforest->forest, apforest->cLocalStart, pforest->forest, 0, &preCoarseToFine, NULL);
890:       } else if (numValuesGlobal) {
891:         p4est_t                   *p4est = pforest->forest;
892:         PetscInt                  *cellFlags;
893:         DMForestAdaptivityStrategy strategy;
894:         PetscSF                    cellSF;
895:         PetscInt                   c, cStart, cEnd;
896:         PetscBool                  adaptAny;

898:         DMForestGetMaximumRefinement(dm, &ctx.maxLevel);
899:         DMForestGetMinimumRefinement(dm, &ctx.minLevel);
900:         DMForestGetAdaptivityStrategy(dm, &strategy);
901:         PetscStrncmp(strategy, "any", 3, &adaptAny);
902:         DMForestGetCellChart(adaptFrom, &cStart, &cEnd);
903:         DMForestGetCellSF(adaptFrom, &cellSF);
904:         PetscMalloc1(cEnd - cStart, &cellFlags);
905:         for (c = cStart; c < cEnd; c++) DMLabelGetValue(adaptLabel, c, &cellFlags[c - cStart]);
906:         if (cellSF) {
907:           if (adaptAny) {
908:             PetscSFReduceBegin(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MAX);
909:             PetscSFReduceEnd(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MAX);
910:           } else {
911:             PetscSFReduceBegin(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MIN);
912:             PetscSFReduceEnd(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MIN);
913:           }
914:         }
915:         for (t = flt, count = cLocalStart; t <= llt; t++) {
916:           p4est_tree_t     *tree     = &(((p4est_tree_t *)p4est->trees->array)[t]);
917:           PetscInt          numQuads = (PetscInt)tree->quadrants.elem_count, i;
918:           p4est_quadrant_t *quads    = (p4est_quadrant_t *)tree->quadrants.array;

920:           for (i = 0; i < numQuads; i++) {
921:             p4est_quadrant_t *q = &quads[i];
922:             q->p.user_int       = cellFlags[count++];
923:           }
924:         }
925:         PetscFree(cellFlags);

927:         pforest->forest->user_pointer = (void *)&ctx;
928:         if (adaptAny) PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_flag_any, pforest_init_determine));
929:         else PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_flag_all, pforest_init_determine));
930:         PetscCallP4est(p4est_refine, (pforest->forest, 0, pforest_refine_flag, NULL));
931:         pforest->forest->user_pointer = (void *)dm;
932:         PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
933:         if (computeAdaptSF) DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), apforest->forest, apforest->cLocalStart, pforest->forest, 0, &preCoarseToFine, &coarseToPreFine);
934:       }
935:       for (t = flt, count = cLocalStart; t <= llt; t++) {
936:         p4est_tree_t     *atree     = &(((p4est_tree_t *)apforest->forest->trees->array)[t]);
937:         p4est_tree_t     *tree      = &(((p4est_tree_t *)pforest->forest->trees->array)[t]);
938:         PetscInt          anumQuads = (PetscInt)atree->quadrants.elem_count, i;
939:         PetscInt          numQuads  = (PetscInt)tree->quadrants.elem_count;
940:         p4est_quadrant_t *aquads    = (p4est_quadrant_t *)atree->quadrants.array;
941:         p4est_quadrant_t *quads     = (p4est_quadrant_t *)tree->quadrants.array;

943:         if (anumQuads != numQuads) {
944:           ctx.anyChange = PETSC_TRUE;
945:         } else {
946:           for (i = 0; i < numQuads; i++) {
947:             p4est_quadrant_t *aq = &aquads[i];
948:             p4est_quadrant_t *q  = &quads[i];

950:             if (aq->level != q->level) {
951:               ctx.anyChange = PETSC_TRUE;
952:               break;
953:             }
954:           }
955:         }
956:         if (ctx.anyChange) break;
957:       }
958:     }
959:     {
960:       PetscInt numLabels, l;

962:       DMGetNumLabels(adaptFrom, &numLabels);
963:       for (l = 0; l < numLabels; l++) {
964:         PetscBool   isDepth, isCellType, isGhost, isVTK;
965:         DMLabel     label, labelNew;
966:         PetscInt    defVal;
967:         const char *name;

969:         DMGetLabelName(adaptFrom, l, &name);
970:         DMGetLabelByNum(adaptFrom, l, &label);
971:         PetscStrcmp(name, "depth", &isDepth);
972:         if (isDepth) continue;
973:         PetscStrcmp(name, "celltype", &isCellType);
974:         if (isCellType) continue;
975:         PetscStrcmp(name, "ghost", &isGhost);
976:         if (isGhost) continue;
977:         PetscStrcmp(name, "vtk", &isVTK);
978:         if (isVTK) continue;
979:         DMCreateLabel(dm, name);
980:         DMGetLabel(dm, name, &labelNew);
981:         DMLabelGetDefaultValue(label, &defVal);
982:         DMLabelSetDefaultValue(labelNew, defVal);
983:       }
984:     }
985:   } else { /* initial */
986:     PetscInt initLevel, minLevel;
987:   #if defined(PETSC_HAVE_MPIUNI)
988:     sc_MPI_Comm comm = sc_MPI_COMM_WORLD;
989:   #else
990:     MPI_Comm comm = PetscObjectComm((PetscObject)dm);
991:   #endif

993:     DMForestGetInitialRefinement(dm, &initLevel);
994:     DMForestGetMinimumRefinement(dm, &minLevel);
995:     PetscCallP4estReturn(pforest->forest, p4est_new_ext,
996:                          (comm, pforest->topo->conn, 0, /* minimum number of quadrants per processor */
997:                           initLevel,                    /* level of refinement */
998:                           1,                            /* uniform refinement */
999:                           0,                            /* we don't allocate any per quadrant data */
1000:                           NULL,                         /* there is no special quadrant initialization */
1001:                           (void *)dm));                 /* this dm is the user context */

1003:     if (initLevel > minLevel) pforest->coarsen_hierarchy = PETSC_TRUE;
1004:     if (dm->setfromoptionscalled) {
1005:       PetscBool   flgPattern, flgFractal;
1006:       PetscInt    corner = 0;
1007:       PetscInt    corners[P4EST_CHILDREN], ncorner = P4EST_CHILDREN;
1008:       PetscReal   likelihood = 1. / P4EST_DIM;
1009:       PetscInt    pattern;
1010:       const char *prefix;

1012:       PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1013:       PetscOptionsGetEList(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_pattern", DMRefinePatternName, PATTERN_COUNT, &pattern, &flgPattern);
1014:       PetscOptionsGetInt(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_corner", &corner, NULL);
1015:       PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_fractal_corners", corners, &ncorner, &flgFractal);
1016:       PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_hash_likelihood", &likelihood, NULL);

1018:       if (flgPattern) {
1019:         DMRefinePatternCtx *ctx;
1020:         PetscInt            maxLevel;

1022:         DMForestGetMaximumRefinement(dm, &maxLevel);
1023:         PetscNew(&ctx);
1024:         ctx->maxLevel = PetscMin(maxLevel, P4EST_QMAXLEVEL);
1025:         if (initLevel + ctx->maxLevel > minLevel) pforest->coarsen_hierarchy = PETSC_TRUE;
1026:         switch (pattern) {
1027:         case PATTERN_HASH:
1028:           ctx->refine_fn      = DMRefinePattern_Hash;
1029:           ctx->hashLikelihood = likelihood;
1030:           break;
1031:         case PATTERN_CORNER:
1032:           ctx->corner    = corner;
1033:           ctx->refine_fn = DMRefinePattern_Corner;
1034:           break;
1035:         case PATTERN_CENTER:
1036:           ctx->refine_fn = DMRefinePattern_Center;
1037:           break;
1038:         case PATTERN_FRACTAL:
1039:           if (flgFractal) {
1040:             PetscInt i;

1042:             for (i = 0; i < ncorner; i++) ctx->fractal[corners[i]] = PETSC_TRUE;
1043:           } else {
1044:   #if !defined(P4_TO_P8)
1045:             ctx->fractal[0] = ctx->fractal[1] = ctx->fractal[2] = PETSC_TRUE;
1046:   #else
1047:             ctx->fractal[0] = ctx->fractal[3] = ctx->fractal[5] = ctx->fractal[6] = PETSC_TRUE;
1048:   #endif
1049:           }
1050:           ctx->refine_fn = DMRefinePattern_Fractal;
1051:           break;
1052:         default:
1053:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Not a valid refinement pattern");
1054:         }

1056:         pforest->forest->user_pointer = (void *)ctx;
1057:         PetscCallP4est(p4est_refine, (pforest->forest, 1, ctx->refine_fn, NULL));
1058:         PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
1059:         PetscFree(ctx);
1060:         pforest->forest->user_pointer = (void *)dm;
1061:       }
1062:     }
1063:   }
1064:   if (pforest->coarsen_hierarchy) {
1065:     PetscInt initLevel, currLevel, minLevel;

1067:     DMPforestGetRefinementLevel(dm, &currLevel);
1068:     DMForestGetInitialRefinement(dm, &initLevel);
1069:     DMForestGetMinimumRefinement(dm, &minLevel);
1070:     if (currLevel > minLevel) {
1071:       DM_Forest_pforest *coarse_pforest;
1072:       DMLabel            coarsen;
1073:       DM                 coarseDM;

1075:       DMForestTemplate(dm, MPI_COMM_NULL, &coarseDM);
1076:       DMForestSetAdaptivityPurpose(coarseDM, DM_ADAPT_COARSEN);
1077:       DMLabelCreate(PETSC_COMM_SELF, "coarsen", &coarsen);
1078:       DMLabelSetDefaultValue(coarsen, DM_ADAPT_COARSEN);
1079:       DMForestSetAdaptivityLabel(coarseDM, coarsen);
1080:       DMLabelDestroy(&coarsen);
1081:       DMSetCoarseDM(dm, coarseDM);
1082:       PetscObjectDereference((PetscObject)coarseDM);
1083:       initLevel = currLevel == initLevel ? initLevel - 1 : initLevel;
1084:       DMForestSetInitialRefinement(coarseDM, initLevel);
1085:       DMForestSetMinimumRefinement(coarseDM, minLevel);
1086:       coarse_pforest                    = (DM_Forest_pforest *)((DM_Forest *)coarseDM->data)->data;
1087:       coarse_pforest->coarsen_hierarchy = PETSC_TRUE;
1088:     }
1089:   }

1091:   { /* repartitioning and overlap */
1092:     PetscMPIInt size, rank;

1094:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
1095:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1096:     if ((size > 1) && (pforest->partition_for_coarsening || forest->cellWeights || forest->weightCapacity != 1. || forest->weightsFactor != 1.)) {
1097:       PetscBool      copyForest  = PETSC_FALSE;
1098:       p4est_t       *forest_copy = NULL;
1099:       p4est_gloidx_t shipped     = 0;

1101:       if (preCoarseToFine || coarseToPreFine) copyForest = PETSC_TRUE;
1102:       if (copyForest) PetscCallP4estReturn(forest_copy, p4est_copy, (pforest->forest, 0));

1104:       if (!forest->cellWeights && forest->weightCapacity == 1. && forest->weightsFactor == 1.) {
1105:         PetscCallP4estReturn(shipped, p4est_partition_ext, (pforest->forest, (int)pforest->partition_for_coarsening, NULL));
1106:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Non-uniform partition cases not implemented yet");
1107:       if (shipped) ctx.anyChange = PETSC_TRUE;
1108:       if (forest_copy) {
1109:         if (preCoarseToFine || coarseToPreFine) {
1110:           PetscSF        repartSF; /* repartSF has roots in the old partition */
1111:           PetscInt       pStart = -1, pEnd = -1, p;
1112:           PetscInt       numRoots, numLeaves;
1113:           PetscSFNode   *repartRoots;
1114:           p4est_gloidx_t postStart  = pforest->forest->global_first_quadrant[rank];
1115:           p4est_gloidx_t postEnd    = pforest->forest->global_first_quadrant[rank + 1];
1116:           p4est_gloidx_t partOffset = postStart;

1118:           numRoots  = (PetscInt)(forest_copy->global_first_quadrant[rank + 1] - forest_copy->global_first_quadrant[rank]);
1119:           numLeaves = (PetscInt)(postEnd - postStart);
1120:           DMPforestComputeOverlappingRanks(size, rank, pforest->forest, forest_copy, &pStart, &pEnd);
1121:           PetscMalloc1((PetscInt)pforest->forest->local_num_quadrants, &repartRoots);
1122:           for (p = pStart; p < pEnd; p++) {
1123:             p4est_gloidx_t preStart = forest_copy->global_first_quadrant[p];
1124:             p4est_gloidx_t preEnd   = forest_copy->global_first_quadrant[p + 1];
1125:             PetscInt       q;

1127:             if (preEnd == preStart) continue;
1129:             preEnd = preEnd > postEnd ? postEnd : preEnd;
1130:             for (q = partOffset; q < preEnd; q++) {
1131:               repartRoots[q - postStart].rank  = p;
1132:               repartRoots[q - postStart].index = partOffset - preStart;
1133:             }
1134:             partOffset = preEnd;
1135:           }
1136:           PetscSFCreate(PetscObjectComm((PetscObject)dm), &repartSF);
1137:           PetscSFSetGraph(repartSF, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, repartRoots, PETSC_OWN_POINTER);
1138:           PetscSFSetUp(repartSF);
1139:           if (preCoarseToFine) {
1140:             PetscSF         repartSFembed, preCoarseToFineNew;
1141:             PetscInt        nleaves;
1142:             const PetscInt *leaves;

1144:             PetscSFSetUp(preCoarseToFine);
1145:             PetscSFGetGraph(preCoarseToFine, NULL, &nleaves, &leaves, NULL);
1146:             if (leaves) {
1147:               PetscSFCreateEmbeddedRootSF(repartSF, nleaves, leaves, &repartSFembed);
1148:             } else {
1149:               repartSFembed = repartSF;
1150:               PetscObjectReference((PetscObject)repartSFembed);
1151:             }
1152:             PetscSFCompose(preCoarseToFine, repartSFembed, &preCoarseToFineNew);
1153:             PetscSFDestroy(&preCoarseToFine);
1154:             PetscSFDestroy(&repartSFembed);
1155:             preCoarseToFine = preCoarseToFineNew;
1156:           }
1157:           if (coarseToPreFine) {
1158:             PetscSF repartSFinv, coarseToPreFineNew;

1160:             PetscSFCreateInverseSF(repartSF, &repartSFinv);
1161:             PetscSFCompose(repartSFinv, coarseToPreFine, &coarseToPreFineNew);
1162:             PetscSFDestroy(&coarseToPreFine);
1163:             PetscSFDestroy(&repartSFinv);
1164:             coarseToPreFine = coarseToPreFineNew;
1165:           }
1166:           PetscSFDestroy(&repartSF);
1167:         }
1168:         PetscCallP4est(p4est_destroy, (forest_copy));
1169:       }
1170:     }
1171:     if (size > 1) {
1172:       PetscInt overlap;

1174:       DMForestGetPartitionOverlap(dm, &overlap);

1176:       if (adaptFrom) {
1177:         PetscInt aoverlap;

1179:         DMForestGetPartitionOverlap(adaptFrom, &aoverlap);
1180:         if (aoverlap != overlap) ctx.anyChange = PETSC_TRUE;
1181:       }

1183:       if (overlap > 0) {
1184:         PetscInt i, cLocalStart;
1185:         PetscInt cEnd;
1186:         PetscSF  preCellSF = NULL, cellSF = NULL;

1188:         PetscCallP4estReturn(pforest->ghost, p4est_ghost_new, (pforest->forest, P4EST_CONNECT_FULL));
1189:         PetscCallP4estReturn(pforest->lnodes, p4est_lnodes_new, (pforest->forest, pforest->ghost, -P4EST_DIM));
1190:         PetscCallP4est(p4est_ghost_support_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));
1191:         for (i = 1; i < overlap; i++) PetscCallP4est(p4est_ghost_expand_by_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));

1193:         cLocalStart = pforest->cLocalStart = pforest->ghost->proc_offsets[rank];
1194:         cEnd                               = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[size];

1196:         /* shift sfs by cLocalStart, expand by cell SFs */
1197:         if (preCoarseToFine || coarseToPreFine) {
1198:           if (adaptFrom) DMForestGetCellSF(adaptFrom, &preCellSF);
1199:           dm->setupcalled = PETSC_TRUE;
1200:           DMForestGetCellSF(dm, &cellSF);
1201:         }
1202:         if (preCoarseToFine) {
1203:           PetscSF            preCoarseToFineNew;
1204:           PetscInt           nleaves, nroots, *leavesNew, i, nleavesNew;
1205:           const PetscInt    *leaves;
1206:           const PetscSFNode *remotes;
1207:           PetscSFNode       *remotesAll;

1209:           PetscSFSetUp(preCoarseToFine);
1210:           PetscSFGetGraph(preCoarseToFine, &nroots, &nleaves, &leaves, &remotes);
1211:           PetscMalloc1(cEnd, &remotesAll);
1212:           for (i = 0; i < cEnd; i++) {
1213:             remotesAll[i].rank  = -1;
1214:             remotesAll[i].index = -1;
1215:           }
1216:           for (i = 0; i < nleaves; i++) remotesAll[(leaves ? leaves[i] : i) + cLocalStart] = remotes[i];
1217:           PetscSFSetUp(cellSF);
1218:           PetscSFBcastBegin(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE);
1219:           PetscSFBcastEnd(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE);
1220:           nleavesNew = 0;
1221:           for (i = 0; i < nleaves; i++) {
1222:             if (remotesAll[i].rank >= 0) nleavesNew++;
1223:           }
1224:           PetscMalloc1(nleavesNew, &leavesNew);
1225:           nleavesNew = 0;
1226:           for (i = 0; i < nleaves; i++) {
1227:             if (remotesAll[i].rank >= 0) {
1228:               leavesNew[nleavesNew] = i;
1229:               if (i > nleavesNew) remotesAll[nleavesNew] = remotesAll[i];
1230:               nleavesNew++;
1231:             }
1232:           }
1233:           PetscSFCreate(PetscObjectComm((PetscObject)dm), &preCoarseToFineNew);
1234:           if (nleavesNew < cEnd) {
1235:             PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, leavesNew, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES);
1236:           } else { /* all cells are leaves */
1237:             PetscFree(leavesNew);
1238:             PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, NULL, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES);
1239:           }
1240:           PetscFree(remotesAll);
1241:           PetscSFDestroy(&preCoarseToFine);
1242:           preCoarseToFine = preCoarseToFineNew;
1243:           preCoarseToFine = preCoarseToFineNew;
1244:         }
1245:         if (coarseToPreFine) {
1246:           PetscSF            coarseToPreFineNew;
1247:           PetscInt           nleaves, nroots, i, nleavesCellSF, nleavesExpanded, *leavesNew;
1248:           const PetscInt    *leaves;
1249:           const PetscSFNode *remotes;
1250:           PetscSFNode       *remotesNew, *remotesNewRoot, *remotesExpanded;

1252:           PetscSFSetUp(coarseToPreFine);
1253:           PetscSFGetGraph(coarseToPreFine, &nroots, &nleaves, &leaves, &remotes);
1254:           PetscSFGetGraph(preCellSF, NULL, &nleavesCellSF, NULL, NULL);
1255:           PetscMalloc1(nroots, &remotesNewRoot);
1256:           PetscMalloc1(nleaves, &remotesNew);
1257:           for (i = 0; i < nroots; i++) {
1258:             remotesNewRoot[i].rank  = rank;
1259:             remotesNewRoot[i].index = i + cLocalStart;
1260:           }
1261:           PetscSFBcastBegin(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE);
1262:           PetscSFBcastEnd(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE);
1263:           PetscFree(remotesNewRoot);
1264:           PetscMalloc1(nleavesCellSF, &remotesExpanded);
1265:           for (i = 0; i < nleavesCellSF; i++) {
1266:             remotesExpanded[i].rank  = -1;
1267:             remotesExpanded[i].index = -1;
1268:           }
1269:           for (i = 0; i < nleaves; i++) remotesExpanded[leaves ? leaves[i] : i] = remotesNew[i];
1270:           PetscFree(remotesNew);
1271:           PetscSFBcastBegin(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE);
1272:           PetscSFBcastEnd(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE);

1274:           nleavesExpanded = 0;
1275:           for (i = 0; i < nleavesCellSF; i++) {
1276:             if (remotesExpanded[i].rank >= 0) nleavesExpanded++;
1277:           }
1278:           PetscMalloc1(nleavesExpanded, &leavesNew);
1279:           nleavesExpanded = 0;
1280:           for (i = 0; i < nleavesCellSF; i++) {
1281:             if (remotesExpanded[i].rank >= 0) {
1282:               leavesNew[nleavesExpanded] = i;
1283:               if (i > nleavesExpanded) remotesExpanded[nleavesExpanded] = remotes[i];
1284:               nleavesExpanded++;
1285:             }
1286:           }
1287:           PetscSFCreate(PetscObjectComm((PetscObject)dm), &coarseToPreFineNew);
1288:           if (nleavesExpanded < nleavesCellSF) {
1289:             PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, leavesNew, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES);
1290:           } else {
1291:             PetscFree(leavesNew);
1292:             PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, NULL, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES);
1293:           }
1294:           PetscFree(remotesExpanded);
1295:           PetscSFDestroy(&coarseToPreFine);
1296:           coarseToPreFine = coarseToPreFineNew;
1297:         }
1298:       }
1299:     }
1300:   }
1301:   forest->preCoarseToFine = preCoarseToFine;
1302:   forest->coarseToPreFine = coarseToPreFine;
1303:   dm->setupcalled         = PETSC_TRUE;
1304:   MPI_Allreduce(&ctx.anyChange, &(pforest->adaptivitySuccess), 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
1305:   DMPforestGetPlex(dm, NULL);
1306:   return 0;
1307: }

1309:   #define DMForestGetAdaptivitySuccess_pforest _append_pforest(DMForestGetAdaptivitySuccess)
1310: static PetscErrorCode DMForestGetAdaptivitySuccess_pforest(DM dm, PetscBool *success)
1311: {
1312:   DM_Forest         *forest;
1313:   DM_Forest_pforest *pforest;

1315:   forest   = (DM_Forest *)dm->data;
1316:   pforest  = (DM_Forest_pforest *)forest->data;
1317:   *success = pforest->adaptivitySuccess;
1318:   return 0;
1319: }

1321:   #define DMView_ASCII_pforest _append_pforest(DMView_ASCII)
1322: static PetscErrorCode DMView_ASCII_pforest(PetscObject odm, PetscViewer viewer)
1323: {
1324:   DM dm = (DM)odm;

1328:   DMSetUp(dm);
1329:   switch (viewer->format) {
1330:   case PETSC_VIEWER_DEFAULT:
1331:   case PETSC_VIEWER_ASCII_INFO: {
1332:     PetscInt    dim;
1333:     const char *name;

1335:     PetscObjectGetName((PetscObject)dm, &name);
1336:     DMGetDimension(dm, &dim);
1337:     if (name) PetscViewerASCIIPrintf(viewer, "Forest %s in %" PetscInt_FMT " dimensions:\n", name, dim);
1338:     else PetscViewerASCIIPrintf(viewer, "Forest in %" PetscInt_FMT " dimensions:\n", dim);
1339:   } /* fall through */
1340:   case PETSC_VIEWER_ASCII_INFO_DETAIL:
1341:   case PETSC_VIEWER_LOAD_BALANCE: {
1342:     DM plex;

1344:     DMPforestGetPlex(dm, &plex);
1345:     DMView(plex, viewer);
1346:   } break;
1347:   default:
1348:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1349:   }
1350:   return 0;
1351: }

1353:   #define DMView_VTK_pforest _append_pforest(DMView_VTK)
1354: static PetscErrorCode DMView_VTK_pforest(PetscObject odm, PetscViewer viewer)
1355: {
1356:   DM                 dm      = (DM)odm;
1357:   DM_Forest         *forest  = (DM_Forest *)dm->data;
1358:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
1359:   PetscBool          isvtk;
1360:   PetscReal          vtkScale = 1. - PETSC_MACHINE_EPSILON;
1361:   PetscViewer_VTK   *vtk      = (PetscViewer_VTK *)viewer->data;
1362:   const char        *name;
1363:   char              *filenameStrip = NULL;
1364:   PetscBool          hasExt;
1365:   size_t             len;
1366:   p4est_geometry_t  *geom;

1370:   DMSetUp(dm);
1371:   geom = pforest->topo->geom;
1372:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1374:   switch (viewer->format) {
1375:   case PETSC_VIEWER_VTK_VTU:
1377:     name = vtk->filename;
1378:     PetscStrlen(name, &len);
1379:     PetscStrcasecmp(name + len - 4, ".vtu", &hasExt);
1380:     if (hasExt) {
1381:       PetscStrallocpy(name, &filenameStrip);
1382:       filenameStrip[len - 4] = '\0';
1383:       name                   = filenameStrip;
1384:     }
1385:     if (!pforest->topo->geom) PetscCallP4estReturn(geom, p4est_geometry_new_connectivity, (pforest->topo->conn));
1386:     {
1387:       p4est_vtk_context_t *pvtk;
1388:       int                  footerr;

1390:       PetscCallP4estReturn(pvtk, p4est_vtk_context_new, (pforest->forest, name));
1391:       PetscCallP4est(p4est_vtk_context_set_geom, (pvtk, geom));
1392:       PetscCallP4est(p4est_vtk_context_set_scale, (pvtk, (double)vtkScale));
1393:       PetscCallP4estReturn(pvtk, p4est_vtk_write_header, (pvtk));
1395:       PetscCallP4estReturn(pvtk, p4est_vtk_write_cell_dataf,
1396:                            (pvtk, 1, /* write tree */
1397:                             1,       /* write level */
1398:                             1,       /* write rank */
1399:                             0,       /* do not wrap rank */
1400:                             0,       /* no scalar fields */
1401:                             0,       /* no vector fields */
1402:                             pvtk));
1404:       PetscCallP4estReturn(footerr, p4est_vtk_write_footer, (pvtk));
1406:     }
1407:     if (!pforest->topo->geom) PetscCallP4est(p4est_geometry_destroy, (geom));
1408:     PetscFree(filenameStrip);
1409:     break;
1410:   default:
1411:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1412:   }
1413:   return 0;
1414: }

1416:   #define DMView_HDF5_pforest _append_pforest(DMView_HDF5)
1417: static PetscErrorCode DMView_HDF5_pforest(DM dm, PetscViewer viewer)
1418: {
1419:   DM plex;

1421:   DMSetUp(dm);
1422:   DMPforestGetPlex(dm, &plex);
1423:   DMView(plex, viewer);
1424:   return 0;
1425: }

1427:   #define DMView_GLVis_pforest _append_pforest(DMView_GLVis)
1428: static PetscErrorCode DMView_GLVis_pforest(DM dm, PetscViewer viewer)
1429: {
1430:   DM plex;

1432:   DMSetUp(dm);
1433:   DMPforestGetPlex(dm, &plex);
1434:   DMView(plex, viewer);
1435:   return 0;
1436: }

1438:   #define DMView_pforest _append_pforest(DMView)
1439: static PetscErrorCode DMView_pforest(DM dm, PetscViewer viewer)
1440: {
1441:   PetscBool isascii, isvtk, ishdf5, isglvis;

1445:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
1446:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1447:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1448:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis);
1449:   if (isascii) {
1450:     DMView_ASCII_pforest((PetscObject)dm, viewer);
1451:   } else if (isvtk) {
1452:     DMView_VTK_pforest((PetscObject)dm, viewer);
1453:   } else if (ishdf5) {
1454:     DMView_HDF5_pforest(dm, viewer);
1455:   } else if (isglvis) {
1456:     DMView_GLVis_pforest(dm, viewer);
1457:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Viewer not supported (not VTK, HDF5, or GLVis)");
1458:   return 0;
1459: }

1461: static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *conn, PetscInt **tree_face_to_uniq)
1462: {
1463:   PetscInt *ttf, f, t, g, count;
1464:   PetscInt  numFacets;

1466:   numFacets = conn->num_trees * P4EST_FACES;
1467:   PetscMalloc1(numFacets, &ttf);
1468:   for (f = 0; f < numFacets; f++) ttf[f] = -1;
1469:   for (g = 0, count = 0, t = 0; t < conn->num_trees; t++) {
1470:     for (f = 0; f < P4EST_FACES; f++, g++) {
1471:       if (ttf[g] == -1) {
1472:         PetscInt ng;

1474:         ttf[g]  = count++;
1475:         ng      = conn->tree_to_tree[g] * P4EST_FACES + (conn->tree_to_face[g] % P4EST_FACES);
1476:         ttf[ng] = ttf[g];
1477:       }
1478:     }
1479:   }
1480:   *tree_face_to_uniq = ttf;
1481:   return 0;
1482: }

1484: static PetscErrorCode DMPlexCreateConnectivity_pforest(DM dm, p4est_connectivity_t **connOut, PetscInt **tree_face_to_uniq)
1485: {
1486:   p4est_topidx_t numTrees, numVerts, numCorns, numCtt;
1487:   PetscSection   ctt;
1488:   #if defined(P4_TO_P8)
1489:   p4est_topidx_t numEdges, numEtt;
1490:   PetscSection   ett;
1491:   PetscInt       eStart, eEnd, e, ettSize;
1492:   PetscInt       vertOff = 1 + P4EST_FACES + P8EST_EDGES;
1493:   PetscInt       edgeOff = 1 + P4EST_FACES;
1494:   #else
1495:   PetscInt vertOff = 1 + P4EST_FACES;
1496:   #endif
1497:   p4est_connectivity_t *conn;
1498:   PetscInt              cStart, cEnd, c, vStart, vEnd, v, fStart, fEnd, f;
1499:   PetscInt             *star = NULL, *closure = NULL, closureSize, starSize, cttSize;
1500:   PetscInt             *ttf;

1502:   /* 1: count objects, allocate */
1503:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1504:   P4estTopidxCast(cEnd - cStart, &numTrees);
1505:   numVerts = P4EST_CHILDREN * numTrees;
1506:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1507:   P4estTopidxCast(vEnd - vStart, &numCorns);
1508:   PetscSectionCreate(PETSC_COMM_SELF, &ctt);
1509:   PetscSectionSetChart(ctt, vStart, vEnd);
1510:   for (v = vStart; v < vEnd; v++) {
1511:     PetscInt s;

1513:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1514:     for (s = 0; s < starSize; s++) {
1515:       PetscInt p = star[2 * s];

1517:       if (p >= cStart && p < cEnd) {
1518:         /* we want to count every time cell p references v, so we see how many times it comes up in the closure.  This
1519:          * only protects against periodicity problems */
1520:         DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1522:         for (c = 0; c < P4EST_CHILDREN; c++) {
1523:           PetscInt cellVert = closure[2 * (c + vertOff)];

1526:           if (cellVert == v) PetscSectionAddDof(ctt, v, 1);
1527:         }
1528:         DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1529:       }
1530:     }
1531:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1532:   }
1533:   PetscSectionSetUp(ctt);
1534:   PetscSectionGetStorageSize(ctt, &cttSize);
1535:   P4estTopidxCast(cttSize, &numCtt);
1536:   #if defined(P4_TO_P8)
1537:   DMPlexGetSimplexOrBoxCells(dm, P4EST_DIM - 1, &eStart, &eEnd);
1538:   P4estTopidxCast(eEnd - eStart, &numEdges);
1539:   PetscSectionCreate(PETSC_COMM_SELF, &ett);
1540:   PetscSectionSetChart(ett, eStart, eEnd);
1541:   for (e = eStart; e < eEnd; e++) {
1542:     PetscInt s;

1544:     DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);
1545:     for (s = 0; s < starSize; s++) {
1546:       PetscInt p = star[2 * s];

1548:       if (p >= cStart && p < cEnd) {
1549:         /* we want to count every time cell p references e, so we see how many times it comes up in the closure.  This
1550:          * only protects against periodicity problems */
1551:         DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1553:         for (c = 0; c < P8EST_EDGES; c++) {
1554:           PetscInt cellEdge = closure[2 * (c + edgeOff)];

1557:           if (cellEdge == e) PetscSectionAddDof(ett, e, 1);
1558:         }
1559:         DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1560:       }
1561:     }
1562:     DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);
1563:   }
1564:   PetscSectionSetUp(ett);
1565:   PetscSectionGetStorageSize(ett, &ettSize);
1566:   P4estTopidxCast(ettSize, &numEtt);

1568:   /* This routine allocates space for the arrays, which we fill below */
1569:   PetscCallP4estReturn(conn, p8est_connectivity_new, (numVerts, numTrees, numEdges, numEtt, numCorns, numCtt));
1570:   #else
1571:   PetscCallP4estReturn(conn, p4est_connectivity_new, (numVerts, numTrees, numCorns, numCtt));
1572:   #endif

1574:   /* 2: visit every face, determine neighboring cells(trees) */
1575:   DMPlexGetSimplexOrBoxCells(dm, 1, &fStart, &fEnd);
1576:   PetscMalloc1((cEnd - cStart) * P4EST_FACES, &ttf);
1577:   for (f = fStart; f < fEnd; f++) {
1578:     PetscInt        numSupp, s;
1579:     PetscInt        myFace[2] = {-1, -1};
1580:     PetscInt        myOrnt[2] = {PETSC_MIN_INT, PETSC_MIN_INT};
1581:     const PetscInt *supp;

1583:     DMPlexGetSupportSize(dm, f, &numSupp);
1585:     DMPlexGetSupport(dm, f, &supp);

1587:     for (s = 0; s < numSupp; s++) {
1588:       PetscInt p = supp[s];

1590:       if (p >= cEnd) {
1591:         numSupp--;
1592:         if (s) supp = &supp[1 - s];
1593:         break;
1594:       }
1595:     }
1596:     for (s = 0; s < numSupp; s++) {
1597:       PetscInt        p = supp[s], i;
1598:       PetscInt        numCone;
1599:       DMPolytopeType  ct;
1600:       const PetscInt *cone;
1601:       const PetscInt *ornt;
1602:       PetscInt        orient = PETSC_MIN_INT;

1604:       DMPlexGetConeSize(dm, p, &numCone);
1606:       DMPlexGetCone(dm, p, &cone);
1607:       DMPlexGetCellType(dm, cone[0], &ct);
1608:       DMPlexGetConeOrientation(dm, p, &ornt);
1609:       for (i = 0; i < P4EST_FACES; i++) {
1610:         if (cone[i] == f) {
1611:           orient = DMPolytopeConvertNewOrientation_Internal(ct, ornt[i]);
1612:           break;
1613:         }
1614:       }
1616:       if (p < cStart || p >= cEnd) {
1617:         DMPolytopeType ct;
1618:         DMPlexGetCellType(dm, p, &ct);
1619:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " (%s) should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], cStart, cEnd);
1620:       }
1621:       ttf[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = f - fStart;
1622:       if (numSupp == 1) {
1623:         /* boundary faces indicated by self reference */
1624:         conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = p - cStart;
1625:         conn->tree_to_face[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = (int8_t)PetscFaceToP4estFace[i];
1626:       } else {
1627:         const PetscInt N = P4EST_CHILDREN / 2;

1629:         conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = supp[1 - s] - cStart;
1630:         myFace[s]                                                                = PetscFaceToP4estFace[i];
1631:         /* get the orientation of cell p in p4est-type closure to facet f, by composing the p4est-closure to
1632:          * petsc-closure permutation and the petsc-closure to facet orientation */
1633:         myOrnt[s] = DihedralCompose(N, orient, DMPolytopeConvertNewOrientation_Internal(ct, P4estFaceToPetscOrnt[myFace[s]]));
1634:       }
1635:     }
1636:     if (numSupp == 2) {
1637:       for (s = 0; s < numSupp; s++) {
1638:         PetscInt       p = supp[s];
1639:         PetscInt       orntAtoB;
1640:         PetscInt       p4estOrient;
1641:         const PetscInt N = P4EST_CHILDREN / 2;

1643:         /* composing the forward permutation with the other cell's inverse permutation gives the self-to-neighbor
1644:          * permutation of this cell-facet's cone */
1645:         orntAtoB = DihedralCompose(N, DihedralInvert(N, myOrnt[1 - s]), myOrnt[s]);

1647:         /* convert cone-description permutation (i.e., edges around facet) to cap-description permutation (i.e.,
1648:          * vertices around facet) */
1649:   #if !defined(P4_TO_P8)
1650:         p4estOrient = orntAtoB < 0 ? -(orntAtoB + 1) : orntAtoB;
1651:   #else
1652:         {
1653:           PetscInt firstVert      = orntAtoB < 0 ? ((-orntAtoB) % N) : orntAtoB;
1654:           PetscInt p4estFirstVert = firstVert < 2 ? firstVert : (firstVert ^ 1);

1656:           /* swap bits */
1657:           p4estOrient = ((myFace[s] <= myFace[1 - s]) || (orntAtoB < 0)) ? p4estFirstVert : ((p4estFirstVert >> 1) | ((p4estFirstVert & 1) << 1));
1658:         }
1659:   #endif
1660:         /* encode neighbor face and orientation in tree_to_face per p4est_connectivity standard (see
1661:          * p4est_connectivity.h, p8est_connectivity.h) */
1662:         conn->tree_to_face[P4EST_FACES * (p - cStart) + myFace[s]] = (int8_t)myFace[1 - s] + p4estOrient * P4EST_FACES;
1663:       }
1664:     }
1665:   }

1667:   #if defined(P4_TO_P8)
1668:   /* 3: visit every edge */
1669:   conn->ett_offset[0] = 0;
1670:   for (e = eStart; e < eEnd; e++) {
1671:     PetscInt off, s;

1673:     PetscSectionGetOffset(ett, e, &off);
1674:     conn->ett_offset[e - eStart] = (p4est_topidx_t)off;
1675:     DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);
1676:     for (s = 0; s < starSize; s++) {
1677:       PetscInt p = star[2 * s];

1679:       if (p >= cStart && p < cEnd) {
1680:         DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1682:         for (c = 0; c < P8EST_EDGES; c++) {
1683:           PetscInt       cellEdge = closure[2 * (c + edgeOff)];
1684:           PetscInt       cellOrnt = closure[2 * (c + edgeOff) + 1];
1685:           DMPolytopeType ct;

1687:           DMPlexGetCellType(dm, cellEdge, &ct);
1688:           cellOrnt = DMPolytopeConvertNewOrientation_Internal(ct, cellOrnt);
1689:           if (cellEdge == e) {
1690:             PetscInt p4estEdge = PetscEdgeToP4estEdge[c];
1691:             PetscInt totalOrient;

1693:             /* compose p4est-closure to petsc-closure permutation and petsc-closure to edge orientation */
1694:             totalOrient = DihedralCompose(2, cellOrnt, DMPolytopeConvertNewOrientation_Internal(DM_POLYTOPE_SEGMENT, P4estEdgeToPetscOrnt[p4estEdge]));
1695:             /* p4est orientations are positive: -2 => 1, -1 => 0 */
1696:             totalOrient             = (totalOrient < 0) ? -(totalOrient + 1) : totalOrient;
1697:             conn->edge_to_tree[off] = (p4est_locidx_t)(p - cStart);
1698:             /* encode cell-edge and orientation in edge_to_edge per p8est_connectivity standart (see
1699:              * p8est_connectivity.h) */
1700:             conn->edge_to_edge[off++]                                  = (int8_t)p4estEdge + P8EST_EDGES * totalOrient;
1701:             conn->tree_to_edge[P8EST_EDGES * (p - cStart) + p4estEdge] = e - eStart;
1702:           }
1703:         }
1704:         DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1705:       }
1706:     }
1707:     DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star);
1708:   }
1709:   PetscSectionDestroy(&ett);
1710:   #endif

1712:   /* 4: visit every vertex */
1713:   conn->ctt_offset[0] = 0;
1714:   for (v = vStart; v < vEnd; v++) {
1715:     PetscInt off, s;

1717:     PetscSectionGetOffset(ctt, v, &off);
1718:     conn->ctt_offset[v - vStart] = (p4est_topidx_t)off;
1719:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1720:     for (s = 0; s < starSize; s++) {
1721:       PetscInt p = star[2 * s];

1723:       if (p >= cStart && p < cEnd) {
1724:         DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1726:         for (c = 0; c < P4EST_CHILDREN; c++) {
1727:           PetscInt cellVert = closure[2 * (c + vertOff)];

1729:           if (cellVert == v) {
1730:             PetscInt p4estVert = PetscVertToP4estVert[c];

1732:             conn->corner_to_tree[off]                                       = (p4est_locidx_t)(p - cStart);
1733:             conn->corner_to_corner[off++]                                   = (int8_t)p4estVert;
1734:             conn->tree_to_corner[P4EST_CHILDREN * (p - cStart) + p4estVert] = v - vStart;
1735:           }
1736:         }
1737:         DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
1738:       }
1739:     }
1740:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1741:   }
1742:   PetscSectionDestroy(&ctt);

1744:   /* 5: Compute the coordinates */
1745:   {
1746:     PetscInt coordDim;

1748:     DMGetCoordinateDim(dm, &coordDim);
1749:     DMGetCoordinatesLocalSetUp(dm);
1750:     for (c = cStart; c < cEnd; c++) {
1751:       PetscInt           dof;
1752:       PetscBool          isDG;
1753:       PetscScalar       *cellCoords = NULL;
1754:       const PetscScalar *array;

1756:       DMPlexGetCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords);
1758:       for (v = 0; v < P4EST_CHILDREN; v++) {
1759:         PetscInt i, lim = PetscMin(3, coordDim);
1760:         PetscInt p4estVert = PetscVertToP4estVert[v];

1762:         conn->tree_to_vertex[P4EST_CHILDREN * (c - cStart) + v] = P4EST_CHILDREN * (c - cStart) + v;
1763:         /* p4est vertices are always embedded in R^3 */
1764:         for (i = 0; i < 3; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = 0.;
1765:         for (i = 0; i < lim; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = PetscRealPart(cellCoords[v * coordDim + i]);
1766:       }
1767:       DMPlexRestoreCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords);
1768:     }
1769:   }

1771:   #if defined(P4EST_ENABLE_DEBUG)
1773:   #endif

1775:   *connOut = conn;

1777:   *tree_face_to_uniq = ttf;

1779:   return 0;
1780: }

1782: static PetscErrorCode locidx_to_PetscInt(sc_array_t *array)
1783: {
1784:   sc_array_t *newarray;
1785:   size_t      zz, count = array->elem_count;


1789:   if (sizeof(p4est_locidx_t) == sizeof(PetscInt)) return 0;

1791:   newarray = sc_array_new_size(sizeof(PetscInt), array->elem_count);
1792:   for (zz = 0; zz < count; zz++) {
1793:     p4est_locidx_t il = *((p4est_locidx_t *)sc_array_index(array, zz));
1794:     PetscInt      *ip = (PetscInt *)sc_array_index(newarray, zz);

1796:     *ip = (PetscInt)il;
1797:   }

1799:   sc_array_reset(array);
1800:   sc_array_init_size(array, sizeof(PetscInt), count);
1801:   sc_array_copy(array, newarray);
1802:   sc_array_destroy(newarray);
1803:   return 0;
1804: }

1806: static PetscErrorCode coords_double_to_PetscScalar(sc_array_t *array, PetscInt dim)
1807: {
1808:   sc_array_t *newarray;
1809:   size_t      zz, count = array->elem_count;

1812:   #if !defined(PETSC_USE_COMPLEX)
1813:   if (sizeof(double) == sizeof(PetscScalar) && dim == 3) return 0;
1814:   #endif

1816:   newarray = sc_array_new_size(dim * sizeof(PetscScalar), array->elem_count);
1817:   for (zz = 0; zz < count; zz++) {
1818:     int          i;
1819:     double      *id = (double *)sc_array_index(array, zz);
1820:     PetscScalar *ip = (PetscScalar *)sc_array_index(newarray, zz);

1822:     for (i = 0; i < dim; i++) ip[i] = 0.;
1823:     for (i = 0; i < PetscMin(dim, 3); i++) ip[i] = (PetscScalar)id[i];
1824:   }

1826:   sc_array_reset(array);
1827:   sc_array_init_size(array, dim * sizeof(PetscScalar), count);
1828:   sc_array_copy(array, newarray);
1829:   sc_array_destroy(newarray);
1830:   return 0;
1831: }

1833: static PetscErrorCode locidx_pair_to_PetscSFNode(sc_array_t *array)
1834: {
1835:   sc_array_t *newarray;
1836:   size_t      zz, count = array->elem_count;


1840:   newarray = sc_array_new_size(sizeof(PetscSFNode), array->elem_count);
1841:   for (zz = 0; zz < count; zz++) {
1842:     p4est_locidx_t *il = (p4est_locidx_t *)sc_array_index(array, zz);
1843:     PetscSFNode    *ip = (PetscSFNode *)sc_array_index(newarray, zz);

1845:     ip->rank  = (PetscInt)il[0];
1846:     ip->index = (PetscInt)il[1];
1847:   }

1849:   sc_array_reset(array);
1850:   sc_array_init_size(array, sizeof(PetscSFNode), count);
1851:   sc_array_copy(array, newarray);
1852:   sc_array_destroy(newarray);
1853:   return 0;
1854: }

1856: static PetscErrorCode P4estToPlex_Local(p4est_t *p4est, DM *plex)
1857: {
1858:   {
1859:     sc_array_t    *points_per_dim    = sc_array_new(sizeof(p4est_locidx_t));
1860:     sc_array_t    *cone_sizes        = sc_array_new(sizeof(p4est_locidx_t));
1861:     sc_array_t    *cones             = sc_array_new(sizeof(p4est_locidx_t));
1862:     sc_array_t    *cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
1863:     sc_array_t    *coords            = sc_array_new(3 * sizeof(double));
1864:     sc_array_t    *children          = sc_array_new(sizeof(p4est_locidx_t));
1865:     sc_array_t    *parents           = sc_array_new(sizeof(p4est_locidx_t));
1866:     sc_array_t    *childids          = sc_array_new(sizeof(p4est_locidx_t));
1867:     sc_array_t    *leaves            = sc_array_new(sizeof(p4est_locidx_t));
1868:     sc_array_t    *remotes           = sc_array_new(2 * sizeof(p4est_locidx_t));
1869:     p4est_locidx_t first_local_quad;

1871:     PetscCallP4est(p4est_get_plex_data, (p4est, P4EST_CONNECT_FULL, 0, &first_local_quad, points_per_dim, cone_sizes, cones, cone_orientations, coords, children, parents, childids, leaves, remotes));

1873:     locidx_to_PetscInt(points_per_dim);
1874:     locidx_to_PetscInt(cone_sizes);
1875:     locidx_to_PetscInt(cones);
1876:     locidx_to_PetscInt(cone_orientations);
1877:     coords_double_to_PetscScalar(coords, P4EST_DIM);

1879:     DMPlexCreate(PETSC_COMM_SELF, plex);
1880:     DMSetDimension(*plex, P4EST_DIM);
1881:     DMPlexCreateFromDAG(*plex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array);
1882:     DMPlexConvertOldOrientations_Internal(*plex);
1883:     sc_array_destroy(points_per_dim);
1884:     sc_array_destroy(cone_sizes);
1885:     sc_array_destroy(cones);
1886:     sc_array_destroy(cone_orientations);
1887:     sc_array_destroy(coords);
1888:     sc_array_destroy(children);
1889:     sc_array_destroy(parents);
1890:     sc_array_destroy(childids);
1891:     sc_array_destroy(leaves);
1892:     sc_array_destroy(remotes);
1893:   }
1894:   return 0;
1895: }

1897:   #define DMReferenceTreeGetChildSymmetry_pforest _append_pforest(DMReferenceTreeGetChildSymmetry)
1898: static PetscErrorCode DMReferenceTreeGetChildSymmetry_pforest(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
1899: {
1900:   PetscInt coneSize, dStart, dEnd, vStart, vEnd, dim, ABswap, oAvert, oBvert, ABswapVert;

1902:   if (parentOrientA == parentOrientB) {
1903:     if (childOrientB) *childOrientB = childOrientA;
1904:     if (childB) *childB = childA;
1905:     return 0;
1906:   }
1907:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1908:   if (childA >= vStart && childA < vEnd) { /* vertices (always in the middle) are invariant under rotation */
1909:     if (childOrientB) *childOrientB = 0;
1910:     if (childB) *childB = childA;
1911:     return 0;
1912:   }
1913:   for (dim = 0; dim < 3; dim++) {
1914:     DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd);
1915:     if (parent >= dStart && parent <= dEnd) break;
1916:   }
1919:   if (childA < dStart || childA >= dEnd) { /* a 1-cell in a 2-cell */
1920:     /* this is a lower-dimensional child: bootstrap */
1921:     PetscInt        size, i, sA = -1, sB, sOrientB, sConeSize;
1922:     const PetscInt *supp, *coneA, *coneB, *oA, *oB;

1924:     DMPlexGetSupportSize(dm, childA, &size);
1925:     DMPlexGetSupport(dm, childA, &supp);

1927:     /* find a point sA in supp(childA) that has the same parent */
1928:     for (i = 0; i < size; i++) {
1929:       PetscInt sParent;

1931:       sA = supp[i];
1932:       if (sA == parent) continue;
1933:       DMPlexGetTreeParent(dm, sA, &sParent, NULL);
1934:       if (sParent == parent) break;
1935:     }
1937:     /* find out which point sB is in an equivalent position to sA under
1938:      * parentOrientB */
1939:     DMReferenceTreeGetChildSymmetry_pforest(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB);
1940:     DMPlexGetConeSize(dm, sA, &sConeSize);
1941:     DMPlexGetCone(dm, sA, &coneA);
1942:     DMPlexGetCone(dm, sB, &coneB);
1943:     DMPlexGetConeOrientation(dm, sA, &oA);
1944:     DMPlexGetConeOrientation(dm, sB, &oB);
1945:     /* step through the cone of sA in natural order */
1946:     for (i = 0; i < sConeSize; i++) {
1947:       if (coneA[i] == childA) {
1948:         /* if childA is at position i in coneA,
1949:          * then we want the point that is at sOrientB*i in coneB */
1950:         PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
1951:         if (childB) *childB = coneB[j];
1952:         if (childOrientB) {
1953:           DMPolytopeType ct;
1954:           PetscInt       oBtrue;

1956:           DMPlexGetConeSize(dm, childA, &coneSize);
1957:           /* compose sOrientB and oB[j] */
1959:           ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1960:           /* we may have to flip an edge */
1961:           oBtrue        = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientationInv(ct, -1, oB[j]);
1962:           oBtrue        = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
1963:           ABswap        = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
1964:           *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
1965:         }
1966:         break;
1967:       }
1968:     }
1970:     return 0;
1971:   }
1972:   /* get the cone size and symmetry swap */
1973:   DMPlexGetConeSize(dm, parent, &coneSize);
1974:   ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
1975:   if (dim == 2) {
1976:     /* orientations refer to cones: we want them to refer to vertices:
1977:      * if it's a rotation, they are the same, but if the order is reversed, a
1978:      * permutation that puts side i first does *not* put vertex i first */
1979:     oAvert     = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
1980:     oBvert     = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
1981:     ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
1982:   } else {
1983:     oAvert     = parentOrientA;
1984:     oBvert     = parentOrientB;
1985:     ABswapVert = ABswap;
1986:   }
1987:   if (childB) {
1988:     /* assume that each child corresponds to a vertex, in the same order */
1989:     PetscInt        p, posA = -1, numChildren, i;
1990:     const PetscInt *children;

1992:     /* count which position the child is in */
1993:     DMPlexGetTreeChildren(dm, parent, &numChildren, &children);
1994:     for (i = 0; i < numChildren; i++) {
1995:       p = children[i];
1996:       if (p == childA) {
1997:         if (dim == 1) {
1998:           posA = i;
1999:         } else { /* 2D Morton to rotation */
2000:           posA = (i & 2) ? (i ^ 1) : i;
2001:         }
2002:         break;
2003:       }
2004:     }
2005:     if (posA >= coneSize) {
2006:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find childA in children of parent");
2007:     } else {
2008:       /* figure out position B by applying ABswapVert */
2009:       PetscInt posB, childIdB;

2011:       posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
2012:       if (dim == 1) {
2013:         childIdB = posB;
2014:       } else { /* 2D rotation to Morton */
2015:         childIdB = (posB & 2) ? (posB ^ 1) : posB;
2016:       }
2017:       if (childB) *childB = children[childIdB];
2018:     }
2019:   }
2020:   if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
2021:   return 0;
2022: }

2024:   #define DMCreateReferenceTree_pforest _append_pforest(DMCreateReferenceTree)
2025: static PetscErrorCode DMCreateReferenceTree_pforest(MPI_Comm comm, DM *dm)
2026: {
2027:   p4est_connectivity_t *refcube;
2028:   p4est_t              *root, *refined;
2029:   DM                    dmRoot, dmRefined;
2030:   DM_Plex              *mesh;
2031:   PetscMPIInt           rank;
2032:   #if defined(PETSC_HAVE_MPIUNI)
2033:   sc_MPI_Comm comm_self = sc_MPI_COMM_SELF;
2034:   #else
2035:   MPI_Comm comm_self = PETSC_COMM_SELF;
2036:   #endif

2038:   PetscCallP4estReturn(refcube, p4est_connectivity_new_byname, ("unit"));
2039:   { /* [-1,1]^d geometry */
2040:     PetscInt i, j;

2042:     for (i = 0; i < P4EST_CHILDREN; i++) {
2043:       for (j = 0; j < 3; j++) {
2044:         refcube->vertices[3 * i + j] *= 2.;
2045:         refcube->vertices[3 * i + j] -= 1.;
2046:       }
2047:     }
2048:   }
2049:   PetscCallP4estReturn(root, p4est_new, (comm_self, refcube, 0, NULL, NULL));
2050:   PetscCallP4estReturn(refined, p4est_new_ext, (comm_self, refcube, 0, 1, 1, 0, NULL, NULL));
2051:   P4estToPlex_Local(root, &dmRoot);
2052:   P4estToPlex_Local(refined, &dmRefined);
2053:   {
2054:   #if !defined(P4_TO_P8)
2055:     PetscInt nPoints   = 25;
2056:     PetscInt perm[25]  = {0, 1, 2, 3, 4, 12, 8, 14, 6, 9, 15, 5, 13, 10, 7, 11, 16, 22, 20, 24, 17, 21, 18, 23, 19};
2057:     PetscInt ident[25] = {0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0, 0, 0, 5, 6, 7, 8, 1, 2, 3, 4, 0};
2058:   #else
2059:     PetscInt nPoints    = 125;
2060:     PetscInt perm[125]  = {0,  1,  2,  3,  4,  5,  6,  7,  8,  32, 16, 36, 24, 40, 12, 17,  37,  25,  41,  9,   33,  20,  26, 42,  13,  21,  27,  43,  10,  34,  18,  38,  28,  14,  19,  39,  29,  11,  35,  22,  30, 15,
2061:                            23, 31, 44, 84, 76, 92, 52, 86, 68, 94, 60, 78, 70, 96, 45, 85,  77,  93,  54,  72,  62,  74,  46, 80,  53,  87,  69,  95,  64,  82,  47,  81,  55,  73,  66,  48,  88,  56,  90,  61,  79, 71,
2062:                            97, 49, 89, 58, 63, 75, 50, 57, 91, 65, 83, 51, 59, 67, 98, 106, 110, 122, 114, 120, 118, 124, 99, 111, 115, 119, 100, 107, 116, 121, 101, 117, 102, 108, 112, 123, 103, 113, 104, 109, 105};
2063:     PetscInt ident[125] = {0,  0,  0,  0,  0,  0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 0, 0, 0, 0,  0,  0,  0,  0,  0,  0,  0,  0, 7, 7, 8,  8,  9,  9,  10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16,
2064:                            16, 17, 17, 18, 18, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 0, 0, 0, 0, 0, 0, 19, 20, 21, 22, 23, 24, 25, 26, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1,  2,  3,  4,  5,  6,  0};

2066:   #endif
2067:     IS permIS;
2068:     DM dmPerm;

2070:     ISCreateGeneral(PETSC_COMM_SELF, nPoints, perm, PETSC_USE_POINTER, &permIS);
2071:     DMPlexPermute(dmRefined, permIS, &dmPerm);
2072:     if (dmPerm) {
2073:       DMDestroy(&dmRefined);
2074:       dmRefined = dmPerm;
2075:     }
2076:     ISDestroy(&permIS);
2077:     {
2078:       PetscInt p;
2079:       DMCreateLabel(dmRoot, "identity");
2080:       DMCreateLabel(dmRefined, "identity");
2081:       for (p = 0; p < P4EST_INSUL; p++) DMSetLabelValue(dmRoot, "identity", p, p);
2082:       for (p = 0; p < nPoints; p++) DMSetLabelValue(dmRefined, "identity", p, ident[p]);
2083:     }
2084:   }
2085:   DMPlexCreateReferenceTree_Union(dmRoot, dmRefined, "identity", dm);
2086:   mesh                   = (DM_Plex *)(*dm)->data;
2087:   mesh->getchildsymmetry = DMReferenceTreeGetChildSymmetry_pforest;
2088:   MPI_Comm_rank(comm, &rank);
2089:   if (rank == 0) {
2090:     DMViewFromOptions(dmRoot, NULL, "-dm_p4est_ref_root_view");
2091:     DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_refined_view");
2092:     DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_tree_view");
2093:   }
2094:   DMDestroy(&dmRefined);
2095:   DMDestroy(&dmRoot);
2096:   PetscCallP4est(p4est_destroy, (refined));
2097:   PetscCallP4est(p4est_destroy, (root));
2098:   PetscCallP4est(p4est_connectivity_destroy, (refcube));
2099:   return 0;
2100: }

2102: static PetscErrorCode DMShareDiscretization(DM dmA, DM dmB)
2103: {
2104:   void     *ctx;
2105:   PetscInt  num;
2106:   PetscReal val;

2108:   DMGetApplicationContext(dmA, &ctx);
2109:   DMSetApplicationContext(dmB, ctx);
2110:   DMCopyDisc(dmA, dmB);
2111:   DMGetOutputSequenceNumber(dmA, &num, &val);
2112:   DMSetOutputSequenceNumber(dmB, num, val);
2113:   if (dmB->localSection != dmA->localSection || dmB->globalSection != dmA->globalSection) {
2114:     DMClearLocalVectors(dmB);
2115:     PetscObjectReference((PetscObject)dmA->localSection);
2116:     PetscSectionDestroy(&(dmB->localSection));
2117:     dmB->localSection = dmA->localSection;
2118:     DMClearGlobalVectors(dmB);
2119:     PetscObjectReference((PetscObject)dmA->globalSection);
2120:     PetscSectionDestroy(&(dmB->globalSection));
2121:     dmB->globalSection = dmA->globalSection;
2122:     PetscObjectReference((PetscObject)dmA->defaultConstraint.section);
2123:     PetscSectionDestroy(&(dmB->defaultConstraint.section));
2124:     dmB->defaultConstraint.section = dmA->defaultConstraint.section;
2125:     PetscObjectReference((PetscObject)dmA->defaultConstraint.mat);
2126:     MatDestroy(&(dmB->defaultConstraint.mat));
2127:     dmB->defaultConstraint.mat = dmA->defaultConstraint.mat;
2128:     if (dmA->map) PetscLayoutReference(dmA->map, &dmB->map);
2129:   }
2130:   if (dmB->sectionSF != dmA->sectionSF) {
2131:     PetscObjectReference((PetscObject)dmA->sectionSF);
2132:     PetscSFDestroy(&dmB->sectionSF);
2133:     dmB->sectionSF = dmA->sectionSF;
2134:   }
2135:   return 0;
2136: }

2138: /* Get an SF that broadcasts a coarse-cell covering of the local fine cells */
2139: static PetscErrorCode DMPforestGetCellCoveringSF(MPI_Comm comm, p4est_t *p4estC, p4est_t *p4estF, PetscInt cStart, PetscInt cEnd, PetscSF *coveringSF)
2140: {
2141:   PetscInt     startF, endF, startC, endC, p, nLeaves;
2142:   PetscSFNode *leaves;
2143:   PetscSF      sf;
2144:   PetscInt    *recv, *send;
2145:   PetscMPIInt  tag;
2146:   MPI_Request *recvReqs, *sendReqs;
2147:   PetscSection section;

2149:   DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estF, p4estC, &startC, &endC);
2150:   PetscMalloc2(2 * (endC - startC), &recv, endC - startC, &recvReqs);
2151:   PetscCommGetNewTag(comm, &tag);
2152:   for (p = startC; p < endC; p++) {
2153:     recvReqs[p - startC] = MPI_REQUEST_NULL;                                        /* just in case we don't initiate a receive */
2154:     if (p4estC->global_first_quadrant[p] == p4estC->global_first_quadrant[p + 1]) { /* empty coarse partition */
2155:       recv[2 * (p - startC)]     = 0;
2156:       recv[2 * (p - startC) + 1] = 0;
2157:       continue;
2158:     }

2160:     MPI_Irecv(&recv[2 * (p - startC)], 2, MPIU_INT, p, tag, comm, &recvReqs[p - startC]);
2161:   }
2162:   DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estC, p4estF, &startF, &endF);
2163:   PetscMalloc2(2 * (endF - startF), &send, endF - startF, &sendReqs);
2164:   /* count the quadrants rank will send to each of [startF,endF) */
2165:   for (p = startF; p < endF; p++) {
2166:     p4est_quadrant_t *myFineStart = &p4estF->global_first_position[p];
2167:     p4est_quadrant_t *myFineEnd   = &p4estF->global_first_position[p + 1];
2168:     PetscInt          tStart      = (PetscInt)myFineStart->p.which_tree;
2169:     PetscInt          tEnd        = (PetscInt)myFineEnd->p.which_tree;
2170:     PetscInt          firstCell = -1, lastCell = -1;
2171:     p4est_tree_t     *treeStart = &(((p4est_tree_t *)p4estC->trees->array)[tStart]);
2172:     p4est_tree_t     *treeEnd   = (size_t)tEnd < p4estC->trees->elem_count ? &(((p4est_tree_t *)p4estC->trees->array)[tEnd]) : NULL;
2173:     ssize_t           overlapIndex;

2175:     sendReqs[p - startF] = MPI_REQUEST_NULL; /* just in case we don't initiate a send */
2176:     if (p4estF->global_first_quadrant[p] == p4estF->global_first_quadrant[p + 1]) continue;

2178:     /* locate myFineStart in (or before) a cell */
2179:     if (treeStart->quadrants.elem_count) {
2180:       PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&(treeStart->quadrants), myFineStart, p4est_quadrant_disjoint));
2181:       if (overlapIndex < 0) {
2182:         firstCell = 0;
2183:       } else {
2184:         firstCell = treeStart->quadrants_offset + overlapIndex;
2185:       }
2186:     } else {
2187:       firstCell = 0;
2188:     }
2189:     if (treeEnd && treeEnd->quadrants.elem_count) {
2190:       PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&(treeEnd->quadrants), myFineEnd, p4est_quadrant_disjoint));
2191:       if (overlapIndex < 0) { /* all of this local section is overlapped */
2192:         lastCell = p4estC->local_num_quadrants;
2193:       } else {
2194:         p4est_quadrant_t *container = &(((p4est_quadrant_t *)treeEnd->quadrants.array)[overlapIndex]);
2195:         p4est_quadrant_t  first_desc;
2196:         int               equal;

2198:         PetscCallP4est(p4est_quadrant_first_descendant, (container, &first_desc, P4EST_QMAXLEVEL));
2199:         PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (myFineEnd, &first_desc));
2200:         if (equal) {
2201:           lastCell = treeEnd->quadrants_offset + overlapIndex;
2202:         } else {
2203:           lastCell = treeEnd->quadrants_offset + overlapIndex + 1;
2204:         }
2205:       }
2206:     } else {
2207:       lastCell = p4estC->local_num_quadrants;
2208:     }
2209:     send[2 * (p - startF)]     = firstCell;
2210:     send[2 * (p - startF) + 1] = lastCell - firstCell;
2211:     MPI_Isend(&send[2 * (p - startF)], 2, MPIU_INT, p, tag, comm, &sendReqs[p - startF]);
2212:   }
2213:   MPI_Waitall((PetscMPIInt)(endC - startC), recvReqs, MPI_STATUSES_IGNORE);
2214:   PetscSectionCreate(PETSC_COMM_SELF, &section);
2215:   PetscSectionSetChart(section, startC, endC);
2216:   for (p = startC; p < endC; p++) {
2217:     PetscInt numCells = recv[2 * (p - startC) + 1];
2218:     PetscSectionSetDof(section, p, numCells);
2219:   }
2220:   PetscSectionSetUp(section);
2221:   PetscSectionGetStorageSize(section, &nLeaves);
2222:   PetscMalloc1(nLeaves, &leaves);
2223:   for (p = startC; p < endC; p++) {
2224:     PetscInt firstCell = recv[2 * (p - startC)];
2225:     PetscInt numCells  = recv[2 * (p - startC) + 1];
2226:     PetscInt off, i;

2228:     PetscSectionGetOffset(section, p, &off);
2229:     for (i = 0; i < numCells; i++) {
2230:       leaves[off + i].rank  = p;
2231:       leaves[off + i].index = firstCell + i;
2232:     }
2233:   }
2234:   PetscSFCreate(comm, &sf);
2235:   PetscSFSetGraph(sf, cEnd - cStart, nLeaves, NULL, PETSC_OWN_POINTER, leaves, PETSC_OWN_POINTER);
2236:   PetscSectionDestroy(&section);
2237:   MPI_Waitall((PetscMPIInt)(endF - startF), sendReqs, MPI_STATUSES_IGNORE);
2238:   PetscFree2(send, sendReqs);
2239:   PetscFree2(recv, recvReqs);
2240:   *coveringSF = sf;
2241:   return 0;
2242: }

2244: /* closure points for locally-owned cells */
2245: static PetscErrorCode DMPforestGetCellSFNodes(DM dm, PetscInt numClosureIndices, PetscInt *numClosurePoints, PetscSFNode **closurePoints, PetscBool redirect)
2246: {
2247:   PetscInt           cStart, cEnd;
2248:   PetscInt           count, c;
2249:   PetscMPIInt        rank;
2250:   PetscInt           closureSize = -1;
2251:   PetscInt          *closure     = NULL;
2252:   PetscSF            pointSF;
2253:   PetscInt           nleaves, nroots;
2254:   const PetscInt    *ilocal;
2255:   const PetscSFNode *iremote;
2256:   DM                 plex;
2257:   DM_Forest         *forest;
2258:   DM_Forest_pforest *pforest;

2260:   forest  = (DM_Forest *)dm->data;
2261:   pforest = (DM_Forest_pforest *)forest->data;
2262:   cStart  = pforest->cLocalStart;
2263:   cEnd    = pforest->cLocalEnd;
2264:   DMPforestGetPlex(dm, &plex);
2265:   DMGetPointSF(dm, &pointSF);
2266:   PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote);
2267:   nleaves           = PetscMax(0, nleaves);
2268:   nroots            = PetscMax(0, nroots);
2269:   *numClosurePoints = numClosureIndices * (cEnd - cStart);
2270:   PetscMalloc1(*numClosurePoints, closurePoints);
2271:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
2272:   for (c = cStart, count = 0; c < cEnd; c++) {
2273:     PetscInt i;
2274:     DMPlexGetTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure);

2276:     for (i = 0; i < numClosureIndices; i++, count++) {
2277:       PetscInt p   = closure[2 * i];
2278:       PetscInt loc = -1;

2280:       PetscFindInt(p, nleaves, ilocal, &loc);
2281:       if (redirect && loc >= 0) {
2282:         (*closurePoints)[count].rank  = iremote[loc].rank;
2283:         (*closurePoints)[count].index = iremote[loc].index;
2284:       } else {
2285:         (*closurePoints)[count].rank  = rank;
2286:         (*closurePoints)[count].index = p;
2287:       }
2288:     }
2289:     DMPlexRestoreTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure);
2290:   }
2291:   return 0;
2292: }

2294: static void MPIAPI DMPforestMaxSFNode(void *a, void *b, PetscMPIInt *len, MPI_Datatype *type)
2295: {
2296:   PetscMPIInt i;

2298:   for (i = 0; i < *len; i++) {
2299:     PetscSFNode *A = (PetscSFNode *)a;
2300:     PetscSFNode *B = (PetscSFNode *)b;

2302:     if (B->rank < 0) *B = *A;
2303:   }
2304: }

2306: static PetscErrorCode DMPforestGetTransferSF_Point(DM coarse, DM fine, PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2307: {
2308:   MPI_Comm           comm;
2309:   PetscMPIInt        rank, size;
2310:   DM_Forest_pforest *pforestC, *pforestF;
2311:   p4est_t           *p4estC, *p4estF;
2312:   PetscInt           numClosureIndices;
2313:   PetscInt           numClosurePointsC, numClosurePointsF;
2314:   PetscSFNode       *closurePointsC, *closurePointsF;
2315:   p4est_quadrant_t  *coverQuads = NULL;
2316:   p4est_quadrant_t **treeQuads;
2317:   PetscInt          *treeQuadCounts;
2318:   MPI_Datatype       nodeType;
2319:   MPI_Datatype       nodeClosureType;
2320:   MPI_Op             sfNodeReduce;
2321:   p4est_topidx_t     fltF, lltF, t;
2322:   DM                 plexC, plexF;
2323:   PetscInt           pStartF, pEndF, pStartC, pEndC;
2324:   PetscBool          saveInCoarse = PETSC_FALSE;
2325:   PetscBool          saveInFine   = PETSC_FALSE;
2326:   PetscBool          formCids     = (childIds != NULL) ? PETSC_TRUE : PETSC_FALSE;
2327:   PetscInt          *cids         = NULL;

2329:   pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2330:   pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2331:   p4estC   = pforestC->forest;
2332:   p4estF   = pforestF->forest;
2334:   comm = PetscObjectComm((PetscObject)coarse);
2335:   MPI_Comm_rank(comm, &rank);
2336:   MPI_Comm_size(comm, &size);
2337:   DMPforestGetPlex(fine, &plexF);
2338:   DMPlexGetChart(plexF, &pStartF, &pEndF);
2339:   DMPforestGetPlex(coarse, &plexC);
2340:   DMPlexGetChart(plexC, &pStartC, &pEndC);
2341:   { /* check if the results have been cached */
2342:     DM adaptCoarse, adaptFine;

2344:     DMForestGetAdaptivityForest(coarse, &adaptCoarse);
2345:     DMForestGetAdaptivityForest(fine, &adaptFine);
2346:     if (adaptCoarse && adaptCoarse->data == fine->data) { /* coarse is adapted from fine */
2347:       if (pforestC->pointSelfToAdaptSF) {
2348:         PetscObjectReference((PetscObject)(pforestC->pointSelfToAdaptSF));
2349:         *sf = pforestC->pointSelfToAdaptSF;
2350:         if (childIds) {
2351:           PetscMalloc1(pEndF - pStartF, &cids);
2352:           PetscArraycpy(cids, pforestC->pointSelfToAdaptCids, pEndF - pStartF);
2353:           *childIds = cids;
2354:         }
2355:         return 0;
2356:       } else {
2357:         saveInCoarse = PETSC_TRUE;
2358:         formCids     = PETSC_TRUE;
2359:       }
2360:     } else if (adaptFine && adaptFine->data == coarse->data) { /* fine is adapted from coarse */
2361:       if (pforestF->pointAdaptToSelfSF) {
2362:         PetscObjectReference((PetscObject)(pforestF->pointAdaptToSelfSF));
2363:         *sf = pforestF->pointAdaptToSelfSF;
2364:         if (childIds) {
2365:           PetscMalloc1(pEndF - pStartF, &cids);
2366:           PetscArraycpy(cids, pforestF->pointAdaptToSelfCids, pEndF - pStartF);
2367:           *childIds = cids;
2368:         }
2369:         return 0;
2370:       } else {
2371:         saveInFine = PETSC_TRUE;
2372:         formCids   = PETSC_TRUE;
2373:       }
2374:     }
2375:   }

2377:   /* count the number of closure points that have dofs and create a list */
2378:   numClosureIndices = P4EST_INSUL;
2379:   /* create the datatype */
2380:   MPI_Type_contiguous(2, MPIU_INT, &nodeType);
2381:   MPI_Type_commit(&nodeType);
2382:   MPI_Op_create(DMPforestMaxSFNode, PETSC_FALSE, &sfNodeReduce);
2383:   MPI_Type_contiguous(numClosureIndices * 2, MPIU_INT, &nodeClosureType);
2384:   MPI_Type_commit(&nodeClosureType);
2385:   /* everything has to go through cells: for each cell, create a list of the sfnodes in its closure */
2386:   /* get lists of closure point SF nodes for every cell */
2387:   DMPforestGetCellSFNodes(coarse, numClosureIndices, &numClosurePointsC, &closurePointsC, PETSC_TRUE);
2388:   DMPforestGetCellSFNodes(fine, numClosureIndices, &numClosurePointsF, &closurePointsF, PETSC_FALSE);
2389:   /* create pointers for tree lists */
2390:   fltF = p4estF->first_local_tree;
2391:   lltF = p4estF->last_local_tree;
2392:   PetscCalloc2(lltF + 1 - fltF, &treeQuads, lltF + 1 - fltF, &treeQuadCounts);
2393:   /* if the partitions don't match, ship the coarse to cover the fine */
2394:   if (size > 1) {
2395:     PetscInt p;

2397:     for (p = 0; p < size; p++) {
2398:       int equal;

2400:       PetscCallP4estReturn(equal, p4est_quadrant_is_equal_piggy, (&p4estC->global_first_position[p], &p4estF->global_first_position[p]));
2401:       if (!equal) break;
2402:     }
2403:     if (p < size) { /* non-matching distribution: send the coarse to cover the fine */
2404:       PetscInt          cStartC, cEndC;
2405:       PetscSF           coveringSF;
2406:       PetscInt          nleaves;
2407:       PetscInt          count;
2408:       PetscSFNode      *newClosurePointsC;
2409:       p4est_quadrant_t *coverQuadsSend;
2410:       p4est_topidx_t    fltC = p4estC->first_local_tree;
2411:       p4est_topidx_t    lltC = p4estC->last_local_tree;
2412:       p4est_topidx_t    t;
2413:       PetscMPIInt       blockSizes[4]   = {P4EST_DIM, 2, 1, 1};
2414:       MPI_Aint          blockOffsets[4] = {offsetof(p4est_quadrant_t, x), offsetof(p4est_quadrant_t, level), offsetof(p4est_quadrant_t, pad16), offsetof(p4est_quadrant_t, p)};
2415:       MPI_Datatype      blockTypes[4]   = {MPI_INT32_T, MPI_INT8_T, MPI_INT16_T, MPI_INT32_T /* p.which_tree */};
2416:       MPI_Datatype      quadStruct, quadType;

2418:       DMPlexGetSimplexOrBoxCells(plexC, 0, &cStartC, &cEndC);
2419:       DMPforestGetCellCoveringSF(comm, p4estC, p4estF, pforestC->cLocalStart, pforestC->cLocalEnd, &coveringSF);
2420:       PetscSFGetGraph(coveringSF, NULL, &nleaves, NULL, NULL);
2421:       PetscMalloc1(numClosureIndices * nleaves, &newClosurePointsC);
2422:       PetscMalloc1(nleaves, &coverQuads);
2423:       PetscMalloc1(cEndC - cStartC, &coverQuadsSend);
2424:       count = 0;
2425:       for (t = fltC; t <= lltC; t++) { /* unfortunately, we need to pack a send array, since quads are not stored packed in p4est */
2426:         p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2427:         PetscInt      q;

2429:         PetscMemcpy(&coverQuadsSend[count], tree->quadrants.array, tree->quadrants.elem_count * sizeof(p4est_quadrant_t));
2430:         for (q = 0; (size_t)q < tree->quadrants.elem_count; q++) coverQuadsSend[count + q].p.which_tree = t;
2431:         count += tree->quadrants.elem_count;
2432:       }
2433:       /* p is of a union type p4est_quadrant_data, but only the p.which_tree field is active at this time. So, we
2434:          have a simple blockTypes[] to use. Note that quadStruct does not count potential padding in array of
2435:          p4est_quadrant_t. We have to call MPI_Type_create_resized() to change upper-bound of quadStruct.
2436:        */
2437:       MPI_Type_create_struct(4, blockSizes, blockOffsets, blockTypes, &quadStruct);
2438:       MPI_Type_create_resized(quadStruct, 0, sizeof(p4est_quadrant_t), &quadType);
2439:       MPI_Type_commit(&quadType);
2440:       PetscSFBcastBegin(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE);
2441:       PetscSFBcastBegin(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE);
2442:       PetscSFBcastEnd(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE);
2443:       PetscSFBcastEnd(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE);
2444:       MPI_Type_free(&quadStruct);
2445:       MPI_Type_free(&quadType);
2446:       PetscFree(coverQuadsSend);
2447:       PetscFree(closurePointsC);
2448:       PetscSFDestroy(&coveringSF);
2449:       closurePointsC = newClosurePointsC;

2451:       /* assign tree quads based on locations in coverQuads */
2452:       {
2453:         PetscInt q;
2454:         for (q = 0; q < nleaves; q++) {
2455:           p4est_locidx_t t = coverQuads[q].p.which_tree;
2456:           if (!treeQuadCounts[t - fltF]++) treeQuads[t - fltF] = &coverQuads[q];
2457:         }
2458:       }
2459:     }
2460:   }
2461:   if (!coverQuads) { /* matching partitions: assign tree quads based on locations in p4est native arrays */
2462:     for (t = fltF; t <= lltF; t++) {
2463:       p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);

2465:       treeQuadCounts[t - fltF] = tree->quadrants.elem_count;
2466:       treeQuads[t - fltF]      = (p4est_quadrant_t *)tree->quadrants.array;
2467:     }
2468:   }

2470:   {
2471:     PetscInt     p;
2472:     PetscInt     cLocalStartF;
2473:     PetscSF      pointSF;
2474:     PetscSFNode *roots;
2475:     PetscInt    *rootType;
2476:     DM           refTree = NULL;
2477:     DMLabel      canonical;
2478:     PetscInt    *childClosures[P4EST_CHILDREN] = {NULL};
2479:     PetscInt    *rootClosure                   = NULL;
2480:     PetscInt     coarseOffset;
2481:     PetscInt     numCoarseQuads;

2483:     PetscMalloc1(pEndF - pStartF, &roots);
2484:     PetscMalloc1(pEndF - pStartF, &rootType);
2485:     DMGetPointSF(fine, &pointSF);
2486:     for (p = pStartF; p < pEndF; p++) {
2487:       roots[p - pStartF].rank  = -1;
2488:       roots[p - pStartF].index = -1;
2489:       rootType[p - pStartF]    = -1;
2490:     }
2491:     if (formCids) {
2492:       PetscInt child;

2494:       PetscMalloc1(pEndF - pStartF, &cids);
2495:       for (p = pStartF; p < pEndF; p++) cids[p - pStartF] = -2;
2496:       DMPlexGetReferenceTree(plexF, &refTree);
2497:       DMPlexGetTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure);
2498:       for (child = 0; child < P4EST_CHILDREN; child++) { /* get the closures of the child cells in the reference tree */
2499:         DMPlexGetTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]);
2500:       }
2501:       DMGetLabel(refTree, "canonical", &canonical);
2502:     }
2503:     cLocalStartF = pforestF->cLocalStart;
2504:     for (t = fltF, coarseOffset = 0, numCoarseQuads = 0; t <= lltF; t++, coarseOffset += numCoarseQuads) {
2505:       p4est_tree_t     *tree         = &(((p4est_tree_t *)p4estF->trees->array)[t]);
2506:       PetscInt          numFineQuads = tree->quadrants.elem_count;
2507:       p4est_quadrant_t *coarseQuads  = treeQuads[t - fltF];
2508:       p4est_quadrant_t *fineQuads    = (p4est_quadrant_t *)tree->quadrants.array;
2509:       PetscInt          i, coarseCount = 0;
2510:       PetscInt          offset = tree->quadrants_offset;
2511:       sc_array_t        coarseQuadsArray;

2513:       numCoarseQuads = treeQuadCounts[t - fltF];
2514:       PetscCallP4est(sc_array_init_data, (&coarseQuadsArray, coarseQuads, sizeof(p4est_quadrant_t), (size_t)numCoarseQuads));
2515:       for (i = 0; i < numFineQuads; i++) {
2516:         PetscInt          c          = i + offset;
2517:         p4est_quadrant_t *quad       = &fineQuads[i];
2518:         p4est_quadrant_t *quadCoarse = NULL;
2519:         ssize_t           disjoint   = -1;

2521:         while (disjoint < 0 && coarseCount < numCoarseQuads) {
2522:           quadCoarse = &coarseQuads[coarseCount];
2523:           PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2524:           if (disjoint < 0) coarseCount++;
2525:         }
2527:         if (quadCoarse->level > quad->level || (quadCoarse->level == quad->level && !transferIdent)) { /* the "coarse" mesh is finer than the fine mesh at the point: continue */
2528:           if (transferIdent) {                                                                         /* find corners */
2529:             PetscInt j = 0;

2531:             do {
2532:               if (j < P4EST_CHILDREN) {
2533:                 p4est_quadrant_t cornerQuad;
2534:                 int              equal;

2536:                 PetscCallP4est(p4est_quadrant_corner_descendant, (quad, &cornerQuad, j, quadCoarse->level));
2537:                 PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (&cornerQuad, quadCoarse));
2538:                 if (equal) {
2539:                   PetscInt    petscJ = P4estVertToPetscVert[j];
2540:                   PetscInt    p      = closurePointsF[numClosureIndices * c + (P4EST_INSUL - P4EST_CHILDREN) + petscJ].index;
2541:                   PetscSFNode q      = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + (P4EST_INSUL - P4EST_CHILDREN) + petscJ];

2543:                   roots[p - pStartF]    = q;
2544:                   rootType[p - pStartF] = PETSC_MAX_INT;
2545:                   cids[p - pStartF]     = -1;
2546:                   j++;
2547:                 }
2548:               }
2549:               coarseCount++;
2550:               disjoint = 1;
2551:               if (coarseCount < numCoarseQuads) {
2552:                 quadCoarse = &coarseQuads[coarseCount];
2553:                 PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2554:               }
2555:             } while (!disjoint);
2556:           }
2557:           continue;
2558:         }
2559:         if (quadCoarse->level == quad->level) { /* same quad present in coarse and fine mesh */
2560:           PetscInt j;
2561:           for (j = 0; j < numClosureIndices; j++) {
2562:             PetscInt p = closurePointsF[numClosureIndices * c + j].index;

2564:             roots[p - pStartF]    = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + j];
2565:             rootType[p - pStartF] = PETSC_MAX_INT; /* unconditionally accept */
2566:             cids[p - pStartF]     = -1;
2567:           }
2568:         } else {
2569:           PetscInt levelDiff                 = quad->level - quadCoarse->level;
2570:           PetscInt proposedCids[P4EST_INSUL] = {0};

2572:           if (formCids) {
2573:             PetscInt  cl;
2574:             PetscInt *pointClosure = NULL;
2575:             int       cid;

2578:             PetscCallP4estReturn(cid, p4est_quadrant_child_id, (quad));
2579:             DMPlexGetTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure);
2580:             for (cl = 0; cl < P4EST_INSUL; cl++) {
2581:               PetscInt       p      = pointClosure[2 * cl];
2582:               PetscInt       point  = childClosures[cid][2 * cl];
2583:               PetscInt       ornt   = childClosures[cid][2 * cl + 1];
2584:               PetscInt       newcid = -1;
2585:               DMPolytopeType ct;

2587:               if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2588:               DMPlexGetCellType(refTree, point, &ct);
2589:               ornt = DMPolytopeConvertNewOrientation_Internal(ct, ornt);
2590:               if (!cl) {
2591:                 newcid = cid + 1;
2592:               } else {
2593:                 PetscInt rcl, parent, parentOrnt = 0;

2595:                 DMPlexGetTreeParent(refTree, point, &parent, NULL);
2596:                 if (parent == point) {
2597:                   newcid = -1;
2598:                 } else if (!parent) { /* in the root */
2599:                   newcid = point;
2600:                 } else {
2601:                   DMPolytopeType rct = DM_POLYTOPE_UNKNOWN;

2603:                   for (rcl = 1; rcl < P4EST_INSUL; rcl++) {
2604:                     if (rootClosure[2 * rcl] == parent) {
2605:                       DMPlexGetCellType(refTree, parent, &rct);
2606:                       parentOrnt = DMPolytopeConvertNewOrientation_Internal(rct, rootClosure[2 * rcl + 1]);
2607:                       break;
2608:                     }
2609:                   }
2611:                   DMPlexReferenceTreeGetChildSymmetry(refTree, parent, parentOrnt, ornt, point, DMPolytopeConvertNewOrientation_Internal(rct, pointClosure[2 * rcl + 1]), NULL, &newcid);
2612:                 }
2613:               }
2614:               if (newcid >= 0) {
2615:                 if (canonical) DMLabelGetValue(canonical, newcid, &newcid);
2616:                 proposedCids[cl] = newcid;
2617:               }
2618:             }
2619:             DMPlexRestoreTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure);
2620:           }
2621:           p4est_qcoord_t coarseBound[2][P4EST_DIM] = {
2622:             {quadCoarse->x, quadCoarse->y,
2623:   #if defined(P4_TO_P8)
2624:              quadCoarse->z
2625:   #endif
2626:             },
2627:             {0}
2628:           };
2629:           p4est_qcoord_t fineBound[2][P4EST_DIM] = {
2630:             {quad->x, quad->y,
2631:   #if defined(P4_TO_P8)
2632:              quad->z
2633:   #endif
2634:             },
2635:             {0}
2636:           };
2637:           PetscInt j;
2638:           for (j = 0; j < P4EST_DIM; j++) { /* get the coordinates of cell boundaries in each direction */
2639:             coarseBound[1][j] = coarseBound[0][j] + P4EST_QUADRANT_LEN(quadCoarse->level);
2640:             fineBound[1][j]   = fineBound[0][j] + P4EST_QUADRANT_LEN(quad->level);
2641:           }
2642:           for (j = 0; j < numClosureIndices; j++) {
2643:             PetscInt    l, p;
2644:             PetscSFNode q;

2646:             p = closurePointsF[numClosureIndices * c + j].index;
2647:             if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2648:             if (j == 0) { /* volume: ancestor is volume */
2649:               l = 0;
2650:             } else if (j < 1 + P4EST_FACES) { /* facet */
2651:               PetscInt face       = PetscFaceToP4estFace[j - 1];
2652:               PetscInt direction  = face / 2;
2653:               PetscInt coarseFace = -1;

2655:               if (coarseBound[face % 2][direction] == fineBound[face % 2][direction]) {
2656:                 coarseFace = face;
2657:                 l          = 1 + P4estFaceToPetscFace[coarseFace];
2658:               } else {
2659:                 l = 0;
2660:               }
2661:   #if defined(P4_TO_P8)
2662:             } else if (j < 1 + P4EST_FACES + P8EST_EDGES) {
2663:               PetscInt  edge       = PetscEdgeToP4estEdge[j - (1 + P4EST_FACES)];
2664:               PetscInt  direction  = edge / 4;
2665:               PetscInt  mod        = edge % 4;
2666:               PetscInt  coarseEdge = -1, coarseFace = -1;
2667:               PetscInt  minDir = PetscMin((direction + 1) % 3, (direction + 2) % 3);
2668:               PetscInt  maxDir = PetscMax((direction + 1) % 3, (direction + 2) % 3);
2669:               PetscBool dirTest[2];

2671:               dirTest[0] = (PetscBool)(coarseBound[mod % 2][minDir] == fineBound[mod % 2][minDir]);
2672:               dirTest[1] = (PetscBool)(coarseBound[mod / 2][maxDir] == fineBound[mod / 2][maxDir]);

2674:               if (dirTest[0] && dirTest[1]) { /* fine edge falls on coarse edge */
2675:                 coarseEdge = edge;
2676:                 l          = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2677:               } else if (dirTest[0]) { /* fine edge falls on a coarse face in the minDir direction */
2678:                 coarseFace = 2 * minDir + (mod % 2);
2679:                 l          = 1 + P4estFaceToPetscFace[coarseFace];
2680:               } else if (dirTest[1]) { /* fine edge falls on a coarse face in the maxDir direction */
2681:                 coarseFace = 2 * maxDir + (mod / 2);
2682:                 l          = 1 + P4estFaceToPetscFace[coarseFace];
2683:               } else {
2684:                 l = 0;
2685:               }
2686:   #endif
2687:             } else {
2688:               PetscInt  vertex = PetscVertToP4estVert[P4EST_CHILDREN - (P4EST_INSUL - j)];
2689:               PetscBool dirTest[P4EST_DIM];
2690:               PetscInt  m;
2691:               PetscInt  numMatch     = 0;
2692:               PetscInt  coarseVertex = -1, coarseFace = -1;
2693:   #if defined(P4_TO_P8)
2694:               PetscInt coarseEdge = -1;
2695:   #endif

2697:               for (m = 0; m < P4EST_DIM; m++) {
2698:                 dirTest[m] = (PetscBool)(coarseBound[(vertex >> m) & 1][m] == fineBound[(vertex >> m) & 1][m]);
2699:                 if (dirTest[m]) numMatch++;
2700:               }
2701:               if (numMatch == P4EST_DIM) { /* vertex on vertex */
2702:                 coarseVertex = vertex;
2703:                 l            = P4EST_INSUL - (P4EST_CHILDREN - P4estVertToPetscVert[coarseVertex]);
2704:               } else if (numMatch == 1) { /* vertex on face */
2705:                 for (m = 0; m < P4EST_DIM; m++) {
2706:                   if (dirTest[m]) {
2707:                     coarseFace = 2 * m + ((vertex >> m) & 1);
2708:                     break;
2709:                   }
2710:                 }
2711:                 l = 1 + P4estFaceToPetscFace[coarseFace];
2712:   #if defined(P4_TO_P8)
2713:               } else if (numMatch == 2) { /* vertex on edge */
2714:                 for (m = 0; m < P4EST_DIM; m++) {
2715:                   if (!dirTest[m]) {
2716:                     PetscInt otherDir1 = (m + 1) % 3;
2717:                     PetscInt otherDir2 = (m + 2) % 3;
2718:                     PetscInt minDir    = PetscMin(otherDir1, otherDir2);
2719:                     PetscInt maxDir    = PetscMax(otherDir1, otherDir2);

2721:                     coarseEdge = m * 4 + 2 * ((vertex >> maxDir) & 1) + ((vertex >> minDir) & 1);
2722:                     break;
2723:                   }
2724:                 }
2725:                 l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2726:   #endif
2727:               } else { /* volume */
2728:                 l = 0;
2729:               }
2730:             }
2731:             q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + l];
2732:             if (l > rootType[p - pStartF]) {
2733:               if (l >= P4EST_INSUL - P4EST_CHILDREN) { /* vertex on vertex: unconditional acceptance */
2734:                 if (transferIdent) {
2735:                   roots[p - pStartF]    = q;
2736:                   rootType[p - pStartF] = PETSC_MAX_INT;
2737:                   if (formCids) cids[p - pStartF] = -1;
2738:                 }
2739:               } else {
2740:                 PetscInt k, thisp = p, limit;

2742:                 roots[p - pStartF]    = q;
2743:                 rootType[p - pStartF] = l;
2744:                 if (formCids) cids[p - pStartF] = proposedCids[j];
2745:                 limit = transferIdent ? levelDiff : (levelDiff - 1);
2746:                 for (k = 0; k < limit; k++) {
2747:                   PetscInt parent;

2749:                   DMPlexGetTreeParent(plexF, thisp, &parent, NULL);
2750:                   if (parent == thisp) break;

2752:                   roots[parent - pStartF]    = q;
2753:                   rootType[parent - pStartF] = PETSC_MAX_INT;
2754:                   if (formCids) cids[parent - pStartF] = -1;
2755:                   thisp = parent;
2756:                 }
2757:               }
2758:             }
2759:           }
2760:         }
2761:       }
2762:     }

2764:     /* now every cell has labeled the points in its closure, so we first make sure everyone agrees by reducing to roots, and the broadcast the agreements */
2765:     if (size > 1) {
2766:       PetscInt *rootTypeCopy, p;

2768:       PetscMalloc1(pEndF - pStartF, &rootTypeCopy);
2769:       PetscArraycpy(rootTypeCopy, rootType, pEndF - pStartF);
2770:       PetscSFReduceBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX);
2771:       PetscSFReduceEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX);
2772:       PetscSFBcastBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE);
2773:       PetscSFBcastEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE);
2774:       for (p = pStartF; p < pEndF; p++) {
2775:         if (rootTypeCopy[p - pStartF] > rootType[p - pStartF]) { /* another process found a root of higher type (e.g. vertex instead of edge), which we want to accept, so nullify this */
2776:           roots[p - pStartF].rank  = -1;
2777:           roots[p - pStartF].index = -1;
2778:         }
2779:         if (formCids && rootTypeCopy[p - pStartF] == PETSC_MAX_INT) { cids[p - pStartF] = -1; /* we have found an antecedent that is the same: no child id */ }
2780:       }
2781:       PetscFree(rootTypeCopy);
2782:       PetscSFReduceBegin(pointSF, nodeType, roots, roots, sfNodeReduce);
2783:       PetscSFReduceEnd(pointSF, nodeType, roots, roots, sfNodeReduce);
2784:       PetscSFBcastBegin(pointSF, nodeType, roots, roots, MPI_REPLACE);
2785:       PetscSFBcastEnd(pointSF, nodeType, roots, roots, MPI_REPLACE);
2786:     }
2787:     PetscFree(rootType);

2789:     {
2790:       PetscInt     numRoots;
2791:       PetscInt     numLeaves;
2792:       PetscInt    *leaves;
2793:       PetscSFNode *iremote;
2794:       /* count leaves */

2796:       numRoots = pEndC - pStartC;

2798:       numLeaves = 0;
2799:       for (p = pStartF; p < pEndF; p++) {
2800:         if (roots[p - pStartF].index >= 0) numLeaves++;
2801:       }
2802:       PetscMalloc1(numLeaves, &leaves);
2803:       PetscMalloc1(numLeaves, &iremote);
2804:       numLeaves = 0;
2805:       for (p = pStartF; p < pEndF; p++) {
2806:         if (roots[p - pStartF].index >= 0) {
2807:           leaves[numLeaves]  = p - pStartF;
2808:           iremote[numLeaves] = roots[p - pStartF];
2809:           numLeaves++;
2810:         }
2811:       }
2812:       PetscFree(roots);
2813:       PetscSFCreate(comm, sf);
2814:       if (numLeaves == (pEndF - pStartF)) {
2815:         PetscFree(leaves);
2816:         PetscSFSetGraph(*sf, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
2817:       } else {
2818:         PetscSFSetGraph(*sf, numRoots, numLeaves, leaves, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
2819:       }
2820:     }
2821:     if (formCids) {
2822:       PetscSF  pointSF;
2823:       PetscInt child;

2825:       DMPlexGetReferenceTree(plexF, &refTree);
2826:       DMGetPointSF(plexF, &pointSF);
2827:       PetscSFReduceBegin(pointSF, MPIU_INT, cids, cids, MPI_MAX);
2828:       PetscSFReduceEnd(pointSF, MPIU_INT, cids, cids, MPI_MAX);
2829:       if (childIds) *childIds = cids;
2830:       for (child = 0; child < P4EST_CHILDREN; child++) DMPlexRestoreTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]);
2831:       DMPlexRestoreTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure);
2832:     }
2833:   }
2834:   if (saveInCoarse) { /* cache results */
2835:     PetscObjectReference((PetscObject)*sf);
2836:     pforestC->pointSelfToAdaptSF = *sf;
2837:     if (!childIds) {
2838:       pforestC->pointSelfToAdaptCids = cids;
2839:     } else {
2840:       PetscMalloc1(pEndF - pStartF, &pforestC->pointSelfToAdaptCids);
2841:       PetscArraycpy(pforestC->pointSelfToAdaptCids, cids, pEndF - pStartF);
2842:     }
2843:   } else if (saveInFine) {
2844:     PetscObjectReference((PetscObject)*sf);
2845:     pforestF->pointAdaptToSelfSF = *sf;
2846:     if (!childIds) {
2847:       pforestF->pointAdaptToSelfCids = cids;
2848:     } else {
2849:       PetscMalloc1(pEndF - pStartF, &pforestF->pointAdaptToSelfCids);
2850:       PetscArraycpy(pforestF->pointAdaptToSelfCids, cids, pEndF - pStartF);
2851:     }
2852:   }
2853:   PetscFree2(treeQuads, treeQuadCounts);
2854:   PetscFree(coverQuads);
2855:   PetscFree(closurePointsC);
2856:   PetscFree(closurePointsF);
2857:   MPI_Type_free(&nodeClosureType);
2858:   MPI_Op_free(&sfNodeReduce);
2859:   MPI_Type_free(&nodeType);
2860:   return 0;
2861: }

2863: /* children are sf leaves of parents */
2864: static PetscErrorCode DMPforestGetTransferSF_Internal(DM coarse, DM fine, const PetscInt dofPerDim[], PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2865: {
2866:   MPI_Comm           comm;
2867:   PetscMPIInt        rank;
2868:   DM_Forest_pforest *pforestC, *pforestF;
2869:   DM                 plexC, plexF;
2870:   PetscInt           pStartC, pEndC, pStartF, pEndF;
2871:   PetscSF            pointTransferSF;
2872:   PetscBool          allOnes = PETSC_TRUE;

2874:   pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2875:   pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2877:   comm = PetscObjectComm((PetscObject)coarse);
2878:   MPI_Comm_rank(comm, &rank);

2880:   {
2881:     PetscInt i;
2882:     for (i = 0; i <= P4EST_DIM; i++) {
2883:       if (dofPerDim[i] != 1) {
2884:         allOnes = PETSC_FALSE;
2885:         break;
2886:       }
2887:     }
2888:   }
2889:   DMPforestGetTransferSF_Point(coarse, fine, &pointTransferSF, transferIdent, childIds);
2890:   if (allOnes) {
2891:     *sf = pointTransferSF;
2892:     return 0;
2893:   }

2895:   DMPforestGetPlex(fine, &plexF);
2896:   DMPlexGetChart(plexF, &pStartF, &pEndF);
2897:   DMPforestGetPlex(coarse, &plexC);
2898:   DMPlexGetChart(plexC, &pStartC, &pEndC);
2899:   {
2900:     PetscInt           numRoots;
2901:     PetscInt           numLeaves;
2902:     const PetscInt    *leaves;
2903:     const PetscSFNode *iremote;
2904:     PetscInt           d;
2905:     PetscSection       leafSection, rootSection;
2906:     /* count leaves */

2908:     PetscSFGetGraph(pointTransferSF, &numRoots, &numLeaves, &leaves, &iremote);
2909:     PetscSectionCreate(PETSC_COMM_SELF, &rootSection);
2910:     PetscSectionCreate(PETSC_COMM_SELF, &leafSection);
2911:     PetscSectionSetChart(rootSection, pStartC, pEndC);
2912:     PetscSectionSetChart(leafSection, pStartF, pEndF);

2914:     for (d = 0; d <= P4EST_DIM; d++) {
2915:       PetscInt startC, endC, e;

2917:       DMPlexGetSimplexOrBoxCells(plexC, P4EST_DIM - d, &startC, &endC);
2918:       for (e = startC; e < endC; e++) PetscSectionSetDof(rootSection, e, dofPerDim[d]);
2919:     }

2921:     for (d = 0; d <= P4EST_DIM; d++) {
2922:       PetscInt startF, endF, e;

2924:       DMPlexGetSimplexOrBoxCells(plexF, P4EST_DIM - d, &startF, &endF);
2925:       for (e = startF; e < endF; e++) PetscSectionSetDof(leafSection, e, dofPerDim[d]);
2926:     }

2928:     PetscSectionSetUp(rootSection);
2929:     PetscSectionSetUp(leafSection);
2930:     {
2931:       PetscInt     nroots, nleaves;
2932:       PetscInt    *mine, i, p;
2933:       PetscInt    *offsets, *offsetsRoot;
2934:       PetscSFNode *remote;

2936:       PetscMalloc1(pEndF - pStartF, &offsets);
2937:       PetscMalloc1(pEndC - pStartC, &offsetsRoot);
2938:       for (p = pStartC; p < pEndC; p++) PetscSectionGetOffset(rootSection, p, &offsetsRoot[p - pStartC]);
2939:       PetscSFBcastBegin(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE);
2940:       PetscSFBcastEnd(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE);
2941:       PetscSectionGetStorageSize(rootSection, &nroots);
2942:       nleaves = 0;
2943:       for (i = 0; i < numLeaves; i++) {
2944:         PetscInt leaf = leaves ? leaves[i] : i;
2945:         PetscInt dof;

2947:         PetscSectionGetDof(leafSection, leaf, &dof);
2948:         nleaves += dof;
2949:       }
2950:       PetscMalloc1(nleaves, &mine);
2951:       PetscMalloc1(nleaves, &remote);
2952:       nleaves = 0;
2953:       for (i = 0; i < numLeaves; i++) {
2954:         PetscInt leaf = leaves ? leaves[i] : i;
2955:         PetscInt dof;
2956:         PetscInt off, j;

2958:         PetscSectionGetDof(leafSection, leaf, &dof);
2959:         PetscSectionGetOffset(leafSection, leaf, &off);
2960:         for (j = 0; j < dof; j++) {
2961:           remote[nleaves].rank  = iremote[i].rank;
2962:           remote[nleaves].index = offsets[leaf] + j;
2963:           mine[nleaves++]       = off + j;
2964:         }
2965:       }
2966:       PetscFree(offsetsRoot);
2967:       PetscFree(offsets);
2968:       PetscSFCreate(comm, sf);
2969:       PetscSFSetGraph(*sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
2970:     }
2971:     PetscSectionDestroy(&leafSection);
2972:     PetscSectionDestroy(&rootSection);
2973:     PetscSFDestroy(&pointTransferSF);
2974:   }
2975:   return 0;
2976: }

2978: static PetscErrorCode DMPforestGetTransferSF(DM dmA, DM dmB, const PetscInt dofPerDim[], PetscSF *sfAtoB, PetscSF *sfBtoA)
2979: {
2980:   DM          adaptA, adaptB;
2981:   DMAdaptFlag purpose;

2983:   DMForestGetAdaptivityForest(dmA, &adaptA);
2984:   DMForestGetAdaptivityForest(dmB, &adaptB);
2985:   /* it is more efficient when the coarser mesh is the first argument: reorder if we know one is coarser than the other */
2986:   if (adaptA && adaptA->data == dmB->data) { /* dmA was adapted from dmB */
2987:     DMForestGetAdaptivityPurpose(dmA, &purpose);
2988:     if (purpose == DM_ADAPT_REFINE) {
2989:       DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB);
2990:       return 0;
2991:     }
2992:   } else if (adaptB && adaptB->data == dmA->data) { /* dmB was adapted from dmA */
2993:     DMForestGetAdaptivityPurpose(dmB, &purpose);
2994:     if (purpose == DM_ADAPT_COARSEN) {
2995:       DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB);
2996:       return 0;
2997:     }
2998:   }
2999:   if (sfAtoB) DMPforestGetTransferSF_Internal(dmA, dmB, dofPerDim, sfAtoB, PETSC_TRUE, NULL);
3000:   if (sfBtoA) DMPforestGetTransferSF_Internal(dmB, dmA, dofPerDim, sfBtoA, (PetscBool)(sfAtoB == NULL), NULL);
3001:   return 0;
3002: }

3004: static PetscErrorCode DMPforestLabelsInitialize(DM dm, DM plex)
3005: {
3006:   DM_Forest         *forest  = (DM_Forest *)dm->data;
3007:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
3008:   PetscInt           cLocalStart, cLocalEnd, cStart, cEnd, fStart, fEnd, eStart, eEnd, vStart, vEnd;
3009:   PetscInt           cStartBase, cEndBase, fStartBase, fEndBase, vStartBase, vEndBase, eStartBase, eEndBase;
3010:   PetscInt           pStart, pEnd, pStartBase, pEndBase, p;
3011:   DM                 base;
3012:   PetscInt          *star      = NULL, starSize;
3013:   DMLabelLink        next      = dm->labels;
3014:   PetscInt           guess     = 0;
3015:   p4est_topidx_t     num_trees = pforest->topo->conn->num_trees;

3017:   pforest->labelsFinalized = PETSC_TRUE;
3018:   cLocalStart              = pforest->cLocalStart;
3019:   cLocalEnd                = pforest->cLocalEnd;
3020:   DMForestGetBaseDM(dm, &base);
3021:   if (!base) {
3022:     if (pforest->ghostName) { /* insert a label to make the boundaries, with stratum values denoting which face of the element touches the boundary */
3023:       p4est_connectivity_t *conn  = pforest->topo->conn;
3024:       p4est_t              *p4est = pforest->forest;
3025:       p4est_tree_t         *trees = (p4est_tree_t *)p4est->trees->array;
3026:       p4est_topidx_t        t, flt = p4est->first_local_tree;
3027:       p4est_topidx_t        llt = pforest->forest->last_local_tree;
3028:       DMLabel               ghostLabel;
3029:       PetscInt              c;

3031:       DMCreateLabel(plex, pforest->ghostName);
3032:       DMGetLabel(plex, pforest->ghostName, &ghostLabel);
3033:       for (c = cLocalStart, t = flt; t <= llt; t++) {
3034:         p4est_tree_t     *tree     = &trees[t];
3035:         p4est_quadrant_t *quads    = (p4est_quadrant_t *)tree->quadrants.array;
3036:         PetscInt          numQuads = (PetscInt)tree->quadrants.elem_count;
3037:         PetscInt          q;

3039:         for (q = 0; q < numQuads; q++, c++) {
3040:           p4est_quadrant_t *quad = &quads[q];
3041:           PetscInt          f;

3043:           for (f = 0; f < P4EST_FACES; f++) {
3044:             p4est_quadrant_t neigh;
3045:             int              isOutside;

3047:             PetscCallP4est(p4est_quadrant_face_neighbor, (quad, f, &neigh));
3048:             PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&neigh));
3049:             if (isOutside) {
3050:               p4est_topidx_t nt;
3051:               PetscInt       nf;

3053:               nt = conn->tree_to_tree[t * P4EST_FACES + f];
3054:               nf = (PetscInt)conn->tree_to_face[t * P4EST_FACES + f];
3055:               nf = nf % P4EST_FACES;
3056:               if (nt == t && nf == f) {
3057:                 PetscInt        plexF = P4estFaceToPetscFace[f];
3058:                 const PetscInt *cone;

3060:                 DMPlexGetCone(plex, c, &cone);
3061:                 DMLabelSetValue(ghostLabel, cone[plexF], plexF + 1);
3062:               }
3063:             }
3064:           }
3065:         }
3066:       }
3067:     }
3068:     return 0;
3069:   }
3070:   DMPlexGetSimplexOrBoxCells(base, 0, &cStartBase, &cEndBase);
3071:   DMPlexGetSimplexOrBoxCells(base, 1, &fStartBase, &fEndBase);
3072:   DMPlexGetSimplexOrBoxCells(base, P4EST_DIM - 1, &eStartBase, &eEndBase);
3073:   DMPlexGetDepthStratum(base, 0, &vStartBase, &vEndBase);

3075:   DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);
3076:   DMPlexGetSimplexOrBoxCells(plex, 1, &fStart, &fEnd);
3077:   DMPlexGetSimplexOrBoxCells(plex, P4EST_DIM - 1, &eStart, &eEnd);
3078:   DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);

3080:   DMPlexGetChart(plex, &pStart, &pEnd);
3081:   DMPlexGetChart(base, &pStartBase, &pEndBase);
3082:   /* go through the mesh: use star to find a quadrant that borders a point.  Use the closure to determine the
3083:    * orientation of the quadrant relative to that point.  Use that to relate the point to the numbering in the base
3084:    * mesh, and extract a label value (since the base mesh is redundantly distributed, can be found locally). */
3085:   while (next) {
3086:     DMLabel     baseLabel;
3087:     DMLabel     label = next->label;
3088:     PetscBool   isDepth, isCellType, isGhost, isVTK, isSpmap;
3089:     const char *name;

3091:     PetscObjectGetName((PetscObject)label, &name);
3092:     PetscStrcmp(name, "depth", &isDepth);
3093:     if (isDepth) {
3094:       next = next->next;
3095:       continue;
3096:     }
3097:     PetscStrcmp(name, "celltype", &isCellType);
3098:     if (isCellType) {
3099:       next = next->next;
3100:       continue;
3101:     }
3102:     PetscStrcmp(name, "ghost", &isGhost);
3103:     if (isGhost) {
3104:       next = next->next;
3105:       continue;
3106:     }
3107:     PetscStrcmp(name, "vtk", &isVTK);
3108:     if (isVTK) {
3109:       next = next->next;
3110:       continue;
3111:     }
3112:     PetscStrcmp(name, "_forest_base_subpoint_map", &isSpmap);
3113:     if (!isSpmap) {
3114:       DMGetLabel(base, name, &baseLabel);
3115:       if (!baseLabel) {
3116:         next = next->next;
3117:         continue;
3118:       }
3119:       DMLabelCreateIndex(baseLabel, pStartBase, pEndBase);
3120:     } else baseLabel = NULL;

3122:     for (p = pStart; p < pEnd; p++) {
3123:       PetscInt          s, c = -1, l;
3124:       PetscInt         *closure = NULL, closureSize;
3125:       p4est_quadrant_t *ghosts  = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3126:       p4est_tree_t     *trees   = (p4est_tree_t *)pforest->forest->trees->array;
3127:       p4est_quadrant_t *q;
3128:       PetscInt          t, val;
3129:       PetscBool         zerosupportpoint = PETSC_FALSE;

3131:       DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star);
3132:       for (s = 0; s < starSize; s++) {
3133:         PetscInt point = star[2 * s];

3135:         if (cStart <= point && point < cEnd) {
3136:           DMPlexGetTransitiveClosure(plex, point, PETSC_TRUE, &closureSize, &closure);
3137:           for (l = 0; l < closureSize; l++) {
3138:             PetscInt qParent = closure[2 * l], q, pp = p, pParent = p;
3139:             do { /* check parents of q */
3140:               q = qParent;
3141:               if (q == p) {
3142:                 c = point;
3143:                 break;
3144:               }
3145:               DMPlexGetTreeParent(plex, q, &qParent, NULL);
3146:             } while (qParent != q);
3147:             if (c != -1) break;
3148:             DMPlexGetTreeParent(plex, pp, &pParent, NULL);
3149:             q = closure[2 * l];
3150:             while (pParent != pp) { /* check parents of p */
3151:               pp = pParent;
3152:               if (pp == q) {
3153:                 c = point;
3154:                 break;
3155:               }
3156:               DMPlexGetTreeParent(plex, pp, &pParent, NULL);
3157:             }
3158:             if (c != -1) break;
3159:           }
3160:           DMPlexRestoreTransitiveClosure(plex, point, PETSC_TRUE, NULL, &closure);
3161:           if (l < closureSize) break;
3162:         } else {
3163:           PetscInt supportSize;

3165:           DMPlexGetSupportSize(plex, point, &supportSize);
3166:           zerosupportpoint = (PetscBool)(zerosupportpoint || !supportSize);
3167:         }
3168:       }
3169:       if (c < 0) {
3170:         const char *prefix;
3171:         PetscBool   print = PETSC_FALSE;

3173:         PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
3174:         PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_forest_print_label_error", &print, NULL);
3175:         if (print) {
3176:           PetscInt i;

3178:           PetscPrintf(PETSC_COMM_SELF, "[%d] Failed to find cell with point %" PetscInt_FMT " in its closure for label %s (starSize %" PetscInt_FMT ")\n", PetscGlobalRank, p, baseLabel ? ((PetscObject)baseLabel)->name : "_forest_base_subpoint_map", starSize);
3179:           for (i = 0; i < starSize; i++) PetscPrintf(PETSC_COMM_SELF, "  star[%" PetscInt_FMT "] = %" PetscInt_FMT ",%" PetscInt_FMT "\n", i, star[2 * i], star[2 * i + 1]);
3180:         }
3181:         DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star);
3182:         if (zerosupportpoint) continue;
3183:         else
3184:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Failed to find cell with point %" PetscInt_FMT " in its closure for label %s. Rerun with -dm_forest_print_label_error for more information", p, baseLabel ? ((PetscObject)baseLabel)->name : "_forest_base_subpoint_map");
3185:       }
3186:       DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star);

3188:       if (c < cLocalStart) {
3189:         /* get from the beginning of the ghost layer */
3190:         q = &(ghosts[c]);
3191:         t = (PetscInt)q->p.which_tree;
3192:       } else if (c < cLocalEnd) {
3193:         PetscInt lo = 0, hi = num_trees;
3194:         /* get from local quadrants: have to find the right tree */

3196:         c -= cLocalStart;

3198:         do {
3199:           p4est_tree_t *tree;

3202:           tree = &trees[guess];
3203:           if (c < tree->quadrants_offset) {
3204:             hi = guess;
3205:           } else if (c < tree->quadrants_offset + (PetscInt)tree->quadrants.elem_count) {
3206:             q = &((p4est_quadrant_t *)tree->quadrants.array)[c - (PetscInt)tree->quadrants_offset];
3207:             t = guess;
3208:             break;
3209:           } else {
3210:             lo = guess + 1;
3211:           }
3212:           guess = lo + (hi - lo) / 2;
3213:         } while (1);
3214:       } else {
3215:         /* get from the end of the ghost layer */
3216:         c -= (cLocalEnd - cLocalStart);

3218:         q = &(ghosts[c]);
3219:         t = (PetscInt)q->p.which_tree;
3220:       }

3222:       if (l == 0) { /* cell */
3223:         if (baseLabel) {
3224:           DMLabelGetValue(baseLabel, t + cStartBase, &val);
3225:         } else {
3226:           val = t + cStartBase;
3227:         }
3228:         DMLabelSetValue(label, p, val);
3229:       } else if (l >= 1 && l < 1 + P4EST_FACES) { /* facet */
3230:         p4est_quadrant_t nq;
3231:         int              isInside;

3233:         l = PetscFaceToP4estFace[l - 1];
3234:         PetscCallP4est(p4est_quadrant_face_neighbor, (q, l, &nq));
3235:         PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3236:         if (isInside) {
3237:           /* this facet is in the interior of a tree, so it inherits the label of the tree */
3238:           if (baseLabel) {
3239:             DMLabelGetValue(baseLabel, t + cStartBase, &val);
3240:           } else {
3241:             val = t + cStartBase;
3242:           }
3243:           DMLabelSetValue(label, p, val);
3244:         } else {
3245:           PetscInt f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + l];

3247:           if (baseLabel) {
3248:             DMLabelGetValue(baseLabel, f + fStartBase, &val);
3249:           } else {
3250:             val = f + fStartBase;
3251:           }
3252:           DMLabelSetValue(label, p, val);
3253:         }
3254:   #if defined(P4_TO_P8)
3255:       } else if (l >= 1 + P4EST_FACES && l < 1 + P4EST_FACES + P8EST_EDGES) { /* edge */
3256:         p4est_quadrant_t nq;
3257:         int              isInside;

3259:         l = PetscEdgeToP4estEdge[l - (1 + P4EST_FACES)];
3260:         PetscCallP4est(p8est_quadrant_edge_neighbor, (q, l, &nq));
3261:         PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3262:         if (isInside) {
3263:           /* this edge is in the interior of a tree, so it inherits the label of the tree */
3264:           if (baseLabel) {
3265:             DMLabelGetValue(baseLabel, t + cStartBase, &val);
3266:           } else {
3267:             val = t + cStartBase;
3268:           }
3269:           DMLabelSetValue(label, p, val);
3270:         } else {
3271:           int isOutsideFace;

3273:           PetscCallP4estReturn(isOutsideFace, p4est_quadrant_is_outside_face, (&nq));
3274:           if (isOutsideFace) {
3275:             PetscInt f;

3277:             if (nq.x < 0) {
3278:               f = 0;
3279:             } else if (nq.x >= P4EST_../../../..4est/.._LEN) {
3280:               f = 1;
3281:             } else if (nq.y < 0) {
3282:               f = 2;
3283:             } else if (nq.y >= P4EST_../../../..4est/.._LEN) {
3284:               f = 3;
3285:             } else if (nq.z < 0) {
3286:               f = 4;
3287:             } else {
3288:               f = 5;
3289:             }
3290:             f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3291:             if (baseLabel) {
3292:               DMLabelGetValue(baseLabel, f + fStartBase, &val);
3293:             } else {
3294:               val = f + fStartBase;
3295:             }
3296:             DMLabelSetValue(label, p, val);
3297:           } else { /* the quadrant edge corresponds to the tree edge */
3298:             PetscInt e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + l];

3300:             if (baseLabel) {
3301:               DMLabelGetValue(baseLabel, e + eStartBase, &val);
3302:             } else {
3303:               val = e + eStartBase;
3304:             }
3305:             DMLabelSetValue(label, p, val);
3306:           }
3307:         }
3308:   #endif
3309:       } else { /* vertex */
3310:         p4est_quadrant_t nq;
3311:         int              isInside;

3313:   #if defined(P4_TO_P8)
3314:         l = PetscVertToP4estVert[l - (1 + P4EST_FACES + P8EST_EDGES)];
3315:   #else
3316:         l = PetscVertToP4estVert[l - (1 + P4EST_FACES)];
3317:   #endif
3318:         PetscCallP4est(p4est_quadrant_corner_neighbor, (q, l, &nq));
3319:         PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3320:         if (isInside) {
3321:           if (baseLabel) {
3322:             DMLabelGetValue(baseLabel, t + cStartBase, &val);
3323:           } else {
3324:             val = t + cStartBase;
3325:           }
3326:           DMLabelSetValue(label, p, val);
3327:         } else {
3328:           int isOutside;

3330:           PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&nq));
3331:           if (isOutside) {
3332:             PetscInt f = -1;

3334:             if (nq.x < 0) {
3335:               f = 0;
3336:             } else if (nq.x >= P4EST_../../../..4est/.._LEN) {
3337:               f = 1;
3338:             } else if (nq.y < 0) {
3339:               f = 2;
3340:             } else if (nq.y >= P4EST_../../../..4est/.._LEN) {
3341:               f = 3;
3342:   #if defined(P4_TO_P8)
3343:             } else if (nq.z < 0) {
3344:               f = 4;
3345:             } else {
3346:               f = 5;
3347:   #endif
3348:             }
3349:             f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3350:             if (baseLabel) {
3351:               DMLabelGetValue(baseLabel, f + fStartBase, &val);
3352:             } else {
3353:               val = f + fStartBase;
3354:             }
3355:             DMLabelSetValue(label, p, val);
3356:             continue;
3357:           }
3358:   #if defined(P4_TO_P8)
3359:           PetscCallP4estReturn(isOutside, p8est_quadrant_is_outside_edge, (&nq));
3360:           if (isOutside) {
3361:             /* outside edge */
3362:             PetscInt e = -1;

3364:             if (nq.x >= 0 && nq.x < P4EST_../../../..4est/.._LEN) {
3365:               if (nq.z < 0) {
3366:                 if (nq.y < 0) {
3367:                   e = 0;
3368:                 } else {
3369:                   e = 1;
3370:                 }
3371:               } else {
3372:                 if (nq.y < 0) {
3373:                   e = 2;
3374:                 } else {
3375:                   e = 3;
3376:                 }
3377:               }
3378:             } else if (nq.y >= 0 && nq.y < P4EST_../../../..4est/.._LEN) {
3379:               if (nq.z < 0) {
3380:                 if (nq.x < 0) {
3381:                   e = 4;
3382:                 } else {
3383:                   e = 5;
3384:                 }
3385:               } else {
3386:                 if (nq.x < 0) {
3387:                   e = 6;
3388:                 } else {
3389:                   e = 7;
3390:                 }
3391:               }
3392:             } else {
3393:               if (nq.y < 0) {
3394:                 if (nq.x < 0) {
3395:                   e = 8;
3396:                 } else {
3397:                   e = 9;
3398:                 }
3399:               } else {
3400:                 if (nq.x < 0) {
3401:                   e = 10;
3402:                 } else {
3403:                   e = 11;
3404:                 }
3405:               }
3406:             }

3408:             e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + e];
3409:             if (baseLabel) {
3410:               DMLabelGetValue(baseLabel, e + eStartBase, &val);
3411:             } else {
3412:               val = e + eStartBase;
3413:             }
3414:             DMLabelSetValue(label, p, val);
3415:             continue;
3416:           }
3417:   #endif
3418:           {
3419:             /* outside vertex: same corner as quadrant corner */
3420:             PetscInt v = pforest->topo->conn->tree_to_corner[P4EST_CHILDREN * t + l];

3422:             if (baseLabel) {
3423:               DMLabelGetValue(baseLabel, v + vStartBase, &val);
3424:             } else {
3425:               val = v + vStartBase;
3426:             }
3427:             DMLabelSetValue(label, p, val);
3428:           }
3429:         }
3430:       }
3431:     }
3432:     next = next->next;
3433:   }
3434:   return 0;
3435: }

3437: static PetscErrorCode DMPforestLabelsFinalize(DM dm, DM plex)
3438: {
3439:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
3440:   DM                 adapt;

3442:   if (pforest->labelsFinalized) return 0;
3443:   pforest->labelsFinalized = PETSC_TRUE;
3444:   DMForestGetAdaptivityForest(dm, &adapt);
3445:   if (!adapt) {
3446:     /* Initialize labels from the base dm */
3447:     DMPforestLabelsInitialize(dm, plex);
3448:   } else {
3449:     PetscInt    dofPerDim[4] = {1, 1, 1, 1};
3450:     PetscSF     transferForward, transferBackward, pointSF;
3451:     PetscInt    pStart, pEnd, pStartA, pEndA;
3452:     PetscInt   *values, *adaptValues;
3453:     DMLabelLink next = adapt->labels;
3454:     DMLabel     adaptLabel;
3455:     DM          adaptPlex;

3457:     DMForestGetAdaptivityLabel(dm, &adaptLabel);
3458:     DMPforestGetPlex(adapt, &adaptPlex);
3459:     DMPforestGetTransferSF(adapt, dm, dofPerDim, &transferForward, &transferBackward);
3460:     DMPlexGetChart(plex, &pStart, &pEnd);
3461:     DMPlexGetChart(adaptPlex, &pStartA, &pEndA);
3462:     PetscMalloc2(pEnd - pStart, &values, pEndA - pStartA, &adaptValues);
3463:     DMGetPointSF(plex, &pointSF);
3464:     if (PetscDefined(USE_DEBUG)) {
3465:       PetscInt p;
3466:       for (p = pStartA; p < pEndA; p++) adaptValues[p - pStartA] = -1;
3467:       for (p = pStart; p < pEnd; p++) values[p - pStart] = -2;
3468:       if (transferForward) {
3469:         PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE);
3470:         PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE);
3471:       }
3472:       if (transferBackward) {
3473:         PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX);
3474:         PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX);
3475:       }
3476:       for (p = pStart; p < pEnd; p++) {
3477:         PetscInt q = p, parent;

3479:         DMPlexGetTreeParent(plex, q, &parent, NULL);
3480:         while (parent != q) {
3481:           if (values[parent] == -2) values[parent] = values[q];
3482:           q = parent;
3483:           DMPlexGetTreeParent(plex, q, &parent, NULL);
3484:         }
3485:       }
3486:       PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX);
3487:       PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX);
3488:       PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE);
3489:       PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE);
3491:     }
3492:     while (next) {
3493:       DMLabel     nextLabel = next->label;
3494:       const char *name;
3495:       PetscBool   isDepth, isCellType, isGhost, isVTK;
3496:       DMLabel     label;
3497:       PetscInt    p;

3499:       PetscObjectGetName((PetscObject)nextLabel, &name);
3500:       PetscStrcmp(name, "depth", &isDepth);
3501:       if (isDepth) {
3502:         next = next->next;
3503:         continue;
3504:       }
3505:       PetscStrcmp(name, "celltype", &isCellType);
3506:       if (isCellType) {
3507:         next = next->next;
3508:         continue;
3509:       }
3510:       PetscStrcmp(name, "ghost", &isGhost);
3511:       if (isGhost) {
3512:         next = next->next;
3513:         continue;
3514:       }
3515:       PetscStrcmp(name, "vtk", &isVTK);
3516:       if (isVTK) {
3517:         next = next->next;
3518:         continue;
3519:       }
3520:       if (nextLabel == adaptLabel) {
3521:         next = next->next;
3522:         continue;
3523:       }
3524:       /* label was created earlier */
3525:       DMGetLabel(dm, name, &label);
3526:       for (p = pStartA; p < pEndA; p++) DMLabelGetValue(nextLabel, p, &adaptValues[p]);
3527:       for (p = pStart; p < pEnd; p++) values[p] = PETSC_MIN_INT;

3529:       if (transferForward) PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE);
3530:       if (transferBackward) PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX);
3531:       if (transferForward) PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE);
3532:       if (transferBackward) PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX);
3533:       for (p = pStart; p < pEnd; p++) {
3534:         PetscInt q = p, parent;

3536:         DMPlexGetTreeParent(plex, q, &parent, NULL);
3537:         while (parent != q) {
3538:           if (values[parent] == PETSC_MIN_INT) values[parent] = values[q];
3539:           q = parent;
3540:           DMPlexGetTreeParent(plex, q, &parent, NULL);
3541:         }
3542:       }
3543:       PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX);
3544:       PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX);
3545:       PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE);
3546:       PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE);

3548:       for (p = pStart; p < pEnd; p++) DMLabelSetValue(label, p, values[p]);
3549:       next = next->next;
3550:     }
3551:     PetscFree2(values, adaptValues);
3552:     PetscSFDestroy(&transferForward);
3553:     PetscSFDestroy(&transferBackward);
3554:     pforest->labelsFinalized = PETSC_TRUE;
3555:   }
3556:   return 0;
3557: }

3559: static PetscErrorCode DMPforestMapCoordinates_Cell(DM plex, p4est_geometry_t *geom, PetscInt cell, p4est_quadrant_t *q, p4est_topidx_t t, p4est_connectivity_t *conn, PetscScalar *coords)
3560: {
3561:   PetscInt     closureSize, c, coordStart, coordEnd, coordDim;
3562:   PetscInt    *closure = NULL;
3563:   PetscSection coordSec;

3565:   DMGetCoordinateSection(plex, &coordSec);
3566:   PetscSectionGetChart(coordSec, &coordStart, &coordEnd);
3567:   DMGetCoordinateDim(plex, &coordDim);
3568:   DMPlexGetTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure);
3569:   for (c = 0; c < closureSize; c++) {
3570:     PetscInt point = closure[2 * c];

3572:     if (point >= coordStart && point < coordEnd) {
3573:       PetscInt dof, off;
3574:       PetscInt nCoords, i;
3575:       PetscSectionGetDof(coordSec, point, &dof);
3577:       nCoords = dof / coordDim;
3578:       PetscSectionGetOffset(coordSec, point, &off);
3579:       for (i = 0; i < nCoords; i++) {
3580:         PetscScalar *coord               = &coords[off + i * coordDim];
3581:         double       coordP4est[3]       = {0.};
3582:         double       coordP4estMapped[3] = {0.};
3583:         PetscInt     j;
3584:         PetscReal    treeCoords[P4EST_CHILDREN][3] = {{0.}};
3585:         PetscReal    eta[3]                        = {0.};
3586:         PetscInt     numRounds                     = 10;
3587:         PetscReal    coordGuess[3]                 = {0.};

3589:         eta[0] = (PetscReal)q->x / (PetscReal)P4EST_../../../..4est/.._LEN;
3590:         eta[1] = (PetscReal)q->y / (PetscReal)P4EST_../../../..4est/.._LEN;
3591:   #if defined(P4_TO_P8)
3592:         eta[2] = (PetscReal)q->z / (PetscReal)P4EST_../../../..4est/.._LEN;
3593:   #endif

3595:         for (j = 0; j < P4EST_CHILDREN; j++) {
3596:           PetscInt k;

3598:           for (k = 0; k < 3; k++) treeCoords[j][k] = conn->vertices[3 * conn->tree_to_vertex[P4EST_CHILDREN * t + j] + k];
3599:         }

3601:         for (j = 0; j < P4EST_CHILDREN; j++) {
3602:           PetscInt  k;
3603:           PetscReal prod = 1.;

3605:           for (k = 0; k < P4EST_DIM; k++) prod *= (j & (1 << k)) ? eta[k] : (1. - eta[k]);
3606:           for (k = 0; k < 3; k++) coordGuess[k] += prod * treeCoords[j][k];
3607:         }

3609:         for (j = 0; j < numRounds; j++) {
3610:           PetscInt dir;

3612:           for (dir = 0; dir < P4EST_DIM; dir++) {
3613:             PetscInt  k;
3614:             PetscReal diff[3];
3615:             PetscReal dXdeta[3] = {0.};
3616:             PetscReal rhs, scale, update;

3618:             for (k = 0; k < 3; k++) diff[k] = coordP4est[k] - coordGuess[k];
3619:             for (k = 0; k < P4EST_CHILDREN; k++) {
3620:               PetscInt  l;
3621:               PetscReal prod = 1.;

3623:               for (l = 0; l < P4EST_DIM; l++) {
3624:                 if (l == dir) {
3625:                   prod *= (k & (1 << l)) ? 1. : -1.;
3626:                 } else {
3627:                   prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3628:                 }
3629:               }
3630:               for (l = 0; l < 3; l++) dXdeta[l] += prod * treeCoords[k][l];
3631:             }
3632:             rhs   = 0.;
3633:             scale = 0;
3634:             for (k = 0; k < 3; k++) {
3635:               rhs += diff[k] * dXdeta[k];
3636:               scale += dXdeta[k] * dXdeta[k];
3637:             }
3638:             update = rhs / scale;
3639:             eta[dir] += update;
3640:             eta[dir] = PetscMin(eta[dir], 1.);
3641:             eta[dir] = PetscMax(eta[dir], 0.);

3643:             coordGuess[0] = coordGuess[1] = coordGuess[2] = 0.;
3644:             for (k = 0; k < P4EST_CHILDREN; k++) {
3645:               PetscInt  l;
3646:               PetscReal prod = 1.;

3648:               for (l = 0; l < P4EST_DIM; l++) prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3649:               for (l = 0; l < 3; l++) coordGuess[l] += prod * treeCoords[k][l];
3650:             }
3651:           }
3652:         }
3653:         for (j = 0; j < 3; j++) coordP4est[j] = (double)eta[j];

3655:         if (geom) {
3656:           (geom->X)(geom, t, coordP4est, coordP4estMapped);
3657:           for (j = 0; j < coordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3658:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded");
3659:       }
3660:     }
3661:   }
3662:   DMPlexRestoreTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure);
3663:   return 0;
3664: }

3666: static PetscErrorCode DMPforestMapCoordinates(DM dm, DM plex)
3667: {
3668:   DM_Forest         *forest;
3669:   DM_Forest_pforest *pforest;
3670:   p4est_geometry_t  *geom;
3671:   PetscInt           cLocalStart, cLocalEnd;
3672:   Vec                coordLocalVec;
3673:   PetscScalar       *coords;
3674:   p4est_topidx_t     flt, llt, t;
3675:   p4est_tree_t      *trees;
3676:   PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
3677:   void *mapCtx;

3679:   forest  = (DM_Forest *)dm->data;
3680:   pforest = (DM_Forest_pforest *)forest->data;
3681:   geom    = pforest->topo->geom;
3682:   DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx);
3683:   if (!geom && !map) return 0;
3684:   DMGetCoordinatesLocal(plex, &coordLocalVec);
3685:   VecGetArray(coordLocalVec, &coords);
3686:   cLocalStart = pforest->cLocalStart;
3687:   cLocalEnd   = pforest->cLocalEnd;
3688:   flt         = pforest->forest->first_local_tree;
3689:   llt         = pforest->forest->last_local_tree;
3690:   trees       = (p4est_tree_t *)pforest->forest->trees->array;
3691:   if (map) { /* apply the map directly to the existing coordinates */
3692:     PetscSection coordSec;
3693:     PetscInt     coordStart, coordEnd, p, coordDim, p4estCoordDim, cStart, cEnd, cEndInterior;
3694:     DM           base;

3696:     DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
3697:     DMPlexGetGhostCellStratum(plex, &cEndInterior, NULL);
3698:     cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3699:     DMForestGetBaseDM(dm, &base);
3700:     DMGetCoordinateSection(plex, &coordSec);
3701:     PetscSectionGetChart(coordSec, &coordStart, &coordEnd);
3702:     DMGetCoordinateDim(plex, &coordDim);
3703:     p4estCoordDim = PetscMin(coordDim, 3);
3704:     for (p = coordStart; p < coordEnd; p++) {
3705:       PetscInt *star = NULL, starSize;
3706:       PetscInt  dof, off, cell = -1, coarsePoint = -1;
3707:       PetscInt  nCoords, i;
3708:       PetscSectionGetDof(coordSec, p, &dof);
3710:       nCoords = dof / coordDim;
3711:       PetscSectionGetOffset(coordSec, p, &off);
3712:       DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star);
3713:       for (i = 0; i < starSize; i++) {
3714:         PetscInt point = star[2 * i];

3716:         if (cStart <= point && point < cEnd) {
3717:           cell = point;
3718:           break;
3719:         }
3720:       }
3721:       DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star);
3722:       if (cell >= 0) {
3723:         if (cell < cLocalStart) {
3724:           p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;

3726:           coarsePoint = ghosts[cell].p.which_tree;
3727:         } else if (cell < cLocalEnd) {
3728:           cell -= cLocalStart;
3729:           for (t = flt; t <= llt; t++) {
3730:             p4est_tree_t *tree = &(trees[t]);

3732:             if (cell >= tree->quadrants_offset && (size_t)cell < tree->quadrants_offset + tree->quadrants.elem_count) {
3733:               coarsePoint = t;
3734:               break;
3735:             }
3736:           }
3737:         } else {
3738:           p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;

3740:           coarsePoint = ghosts[cell - cLocalEnd].p.which_tree;
3741:         }
3742:       }
3743:       for (i = 0; i < nCoords; i++) {
3744:         PetscScalar *coord               = &coords[off + i * coordDim];
3745:         PetscReal    coordP4est[3]       = {0.};
3746:         PetscReal    coordP4estMapped[3] = {0.};
3747:         PetscInt     j;

3749:         for (j = 0; j < p4estCoordDim; j++) coordP4est[j] = PetscRealPart(coord[j]);
3750:         (map)(base, coarsePoint, p4estCoordDim, coordP4est, coordP4estMapped, mapCtx);
3751:         for (j = 0; j < p4estCoordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3752:       }
3753:     }
3754:   } else { /* we have to transform coordinates back to the unit cube (where geom is defined), and then apply geom */
3755:     PetscInt cStart, cEnd, cEndInterior;

3757:     DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
3758:     DMPlexGetGhostCellStratum(plex, &cEndInterior, NULL);
3759:     cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3760:     if (cLocalStart > 0) {
3761:       p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3762:       PetscInt          count;

3764:       for (count = 0; count < cLocalStart; count++) {
3765:         p4est_quadrant_t *quad = &ghosts[count];
3766:         p4est_topidx_t    t    = quad->p.which_tree;

3768:         DMPforestMapCoordinates_Cell(plex, geom, count, quad, t, pforest->topo->conn, coords);
3769:       }
3770:     }
3771:     for (t = flt; t <= llt; t++) {
3772:       p4est_tree_t     *tree     = &(trees[t]);
3773:       PetscInt          offset   = cLocalStart + tree->quadrants_offset, i;
3774:       PetscInt          numQuads = (PetscInt)tree->quadrants.elem_count;
3775:       p4est_quadrant_t *quads    = (p4est_quadrant_t *)tree->quadrants.array;

3777:       for (i = 0; i < numQuads; i++) {
3778:         PetscInt count = i + offset;

3780:         DMPforestMapCoordinates_Cell(plex, geom, count, &quads[i], t, pforest->topo->conn, coords);
3781:       }
3782:     }
3783:     if (cLocalEnd - cLocalStart < cEnd - cStart) {
3784:       p4est_quadrant_t *ghosts    = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3785:       PetscInt          numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
3786:       PetscInt          count;

3788:       for (count = 0; count < numGhosts - cLocalStart; count++) {
3789:         p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
3790:         p4est_topidx_t    t    = quad->p.which_tree;

3792:         DMPforestMapCoordinates_Cell(plex, geom, count + cLocalEnd, quad, t, pforest->topo->conn, coords);
3793:       }
3794:     }
3795:   }
3796:   VecRestoreArray(coordLocalVec, &coords);
3797:   return 0;
3798: }

3800: static PetscErrorCode PforestQuadrantIsInterior(p4est_quadrant_t *quad, PetscBool *is_interior)
3801: {
3802:   p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3803:   if ((quad->x > 0) && (quad->x + h < P4EST_../../../..4est/.._LEN)
3804:   #ifdef P4_TO_P8
3805:       && (quad->z > 0) && (quad->z + h < P4EST_../../../..4est/.._LEN)
3806:   #endif
3807:       && (quad->y > 0) && (quad->y + h < P4EST_../../../..4est/.._LEN)) {
3808:     *is_interior = PETSC_TRUE;
3809:   } else {
3810:     *is_interior = PETSC_FALSE;
3811:   }
3812:   return 0;
3813: }

3815: /* We always use DG coordinates with p4est: if they do not match the vertex
3816:    coordinates, add space for them in the section */
3817: static PetscErrorCode PforestCheckLocalizeCell(DM plex, PetscInt cDim, Vec cVecOld, DM_Forest_pforest *pforest, PetscSection oldSection, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad)
3818: {
3819:   PetscBool is_interior;

3821:   PforestQuadrantIsInterior(quad, &is_interior);
3822:   if (is_interior) { // quads in the interior of a coarse cell can't touch periodic interfaces
3823:     PetscSectionSetDof(newSection, cell, 0);
3824:     PetscSectionSetFieldDof(newSection, cell, 0, 0);
3825:   } else {
3826:     PetscInt     cSize;
3827:     PetscScalar *values      = NULL;
3828:     PetscBool    same_coords = PETSC_TRUE;

3830:     DMPlexVecGetClosure(plex, oldSection, cVecOld, cell, &cSize, &values);
3831:     PetscAssert(cSize == cDim * P4EST_CHILDREN, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected closure size");
3832:     for (int c = 0; c < P4EST_CHILDREN; c++) {
3833:       p4est_qcoord_t quad_coords[3];
3834:       p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3835:       double         corner_coords[3];
3836:       double         vert_coords[3];
3837:       PetscInt       corner = PetscVertToP4estVert[c];

3839:       for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) vert_coords[d] = PetscRealPart(values[c * cDim + d]);

3841:       quad_coords[0] = quad->x;
3842:       quad_coords[1] = quad->y;
3843:   #ifdef P4_TO_P8
3844:       quad_coords[2] = quad->z;
3845:   #endif
3846:       for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3847:   #ifndef P4_TO_P8
3848:       PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3849:   #else
3850:       PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3851:   #endif
3852:       for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) {
3853:         if (fabs(vert_coords[d] - corner_coords[d]) > PETSC_SMALL) {
3854:           same_coords = PETSC_FALSE;
3855:           break;
3856:         }
3857:       }
3858:     }
3859:     if (same_coords) {
3860:       PetscSectionSetDof(newSection, cell, 0);
3861:       PetscSectionSetFieldDof(newSection, cell, 0, 0);
3862:     } else {
3863:       PetscSectionSetDof(newSection, cell, cSize);
3864:       PetscSectionSetFieldDof(newSection, cell, 0, cSize);
3865:     }
3866:     DMPlexVecRestoreClosure(plex, oldSection, cVecOld, cell, &cSize, &values);
3867:   }
3868:   return 0;
3869: }

3871: static PetscErrorCode PforestLocalizeCell(DM plex, PetscInt cDim, DM_Forest_pforest *pforest, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad, PetscScalar coords[])
3872: {
3873:   PetscInt cdof, off;

3875:   PetscSectionGetDof(newSection, cell, &cdof);
3876:   if (!cdof) return 0;

3878:   PetscSectionGetOffset(newSection, cell, &off);
3879:   for (PetscInt c = 0, pos = off; c < P4EST_CHILDREN; c++) {
3880:     p4est_qcoord_t quad_coords[3];
3881:     p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3882:     double         corner_coords[3];
3883:     PetscInt       corner = PetscVertToP4estVert[c];

3885:     quad_coords[0] = quad->x;
3886:     quad_coords[1] = quad->y;
3887:   #ifdef P4_TO_P8
3888:     quad_coords[2] = quad->z;
3889:   #endif
3890:     for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3891:   #ifndef P4_TO_P8
3892:     PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3893:   #else
3894:     PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3895:   #endif
3896:     for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) coords[pos++] = corner_coords[d];
3897:     for (PetscInt d = PetscMin(cDim, 3); d < cDim; d++) coords[pos++] = 0.;
3898:   }
3899:   return 0;
3900: }

3902: static PetscErrorCode DMPforestLocalizeCoordinates(DM dm, DM plex)
3903: {
3904:   DM_Forest         *forest;
3905:   DM_Forest_pforest *pforest;
3906:   DM                 base, cdm, cdmCell;
3907:   Vec                cVec, cVecOld;
3908:   PetscSection       oldSection, newSection;
3909:   PetscScalar       *coords2;
3910:   const PetscReal   *L;
3911:   PetscInt           cLocalStart, cLocalEnd, coarsePoint;
3912:   PetscInt           cDim, newStart, newEnd;
3913:   PetscInt           v, vStart, vEnd, cp, cStart, cEnd, cEndInterior;
3914:   p4est_topidx_t     flt, llt, t;
3915:   p4est_tree_t      *trees;
3916:   PetscBool          baseLocalized = PETSC_FALSE;

3918:   DMGetPeriodicity(dm, NULL, NULL, &L);
3919:   /* we localize on all cells if we don't have a base DM or the base DM coordinates have not been localized */
3920:   DMGetCoordinateDim(dm, &cDim);
3921:   DMForestGetBaseDM(dm, &base);
3922:   if (base) DMGetCoordinatesLocalized(base, &baseLocalized);
3923:   if (!baseLocalized) base = NULL;
3924:   if (!baseLocalized && !L) return 0;
3925:   DMPlexGetChart(plex, &newStart, &newEnd);

3927:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newSection);
3928:   PetscSectionSetNumFields(newSection, 1);
3929:   PetscSectionSetFieldComponents(newSection, 0, cDim);
3930:   PetscSectionSetChart(newSection, newStart, newEnd);

3932:   DMGetCoordinateSection(plex, &oldSection);
3933:   DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);
3934:   DMGetCoordinatesLocal(plex, &cVecOld);

3936:   forest      = (DM_Forest *)dm->data;
3937:   pforest     = (DM_Forest_pforest *)forest->data;
3938:   cLocalStart = pforest->cLocalStart;
3939:   cLocalEnd   = pforest->cLocalEnd;
3940:   flt         = pforest->forest->first_local_tree;
3941:   llt         = pforest->forest->last_local_tree;
3942:   trees       = (p4est_tree_t *)pforest->forest->trees->array;

3944:   DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
3945:   DMPlexGetGhostCellStratum(plex, &cEndInterior, NULL);
3946:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3947:   cp   = 0;
3948:   if (cLocalStart > 0) {
3949:     p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3950:     PetscInt          cell;

3952:     for (cell = 0; cell < cLocalStart; ++cell, cp++) {
3953:       p4est_quadrant_t *quad = &ghosts[cell];

3955:       coarsePoint = quad->p.which_tree;
3956:       PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad);
3957:     }
3958:   }
3959:   for (t = flt; t <= llt; t++) {
3960:     p4est_tree_t     *tree     = &(trees[t]);
3961:     PetscInt          offset   = cLocalStart + tree->quadrants_offset;
3962:     PetscInt          numQuads = (PetscInt)tree->quadrants.elem_count;
3963:     p4est_quadrant_t *quads    = (p4est_quadrant_t *)tree->quadrants.array;
3964:     PetscInt          i;

3966:     if (!numQuads) continue;
3967:     coarsePoint = t;
3968:     for (i = 0; i < numQuads; i++, cp++) {
3969:       PetscInt          cell = i + offset;
3970:       p4est_quadrant_t *quad = &quads[i];

3972:       PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad);
3973:     }
3974:   }
3975:   if (cLocalEnd - cLocalStart < cEnd - cStart) {
3976:     p4est_quadrant_t *ghosts    = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3977:     PetscInt          numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
3978:     PetscInt          count;

3980:     for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
3981:       p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
3982:       coarsePoint            = quad->p.which_tree;
3983:       PetscInt cell          = count + cLocalEnd;

3985:       PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad);
3986:     }
3987:   }
3988:   PetscAssert(cp == cEnd - cStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of fine cells %" PetscInt_FMT " != %" PetscInt_FMT, cp, cEnd - cStart);

3990:   PetscSectionSetUp(newSection);
3991:   DMGetCoordinateDM(plex, &cdm);
3992:   DMClone(cdm, &cdmCell);
3993:   DMSetCellCoordinateDM(plex, cdmCell);
3994:   DMDestroy(&cdmCell);
3995:   DMSetCellCoordinateSection(plex, cDim, newSection);
3996:   PetscSectionGetStorageSize(newSection, &v);
3997:   VecCreate(PETSC_COMM_SELF, &cVec);
3998:   PetscObjectSetName((PetscObject)cVec, "coordinates");
3999:   VecSetBlockSize(cVec, cDim);
4000:   VecSetSizes(cVec, v, PETSC_DETERMINE);
4001:   VecSetType(cVec, VECSTANDARD);
4002:   VecSet(cVec, PETSC_MIN_REAL);

4004:   /* Localize coordinates on cells if needed */
4005:   VecGetArray(cVec, &coords2);
4006:   cp = 0;
4007:   if (cLocalStart > 0) {
4008:     p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4009:     PetscInt          cell;

4011:     for (cell = 0; cell < cLocalStart; ++cell, cp++) {
4012:       p4est_quadrant_t *quad = &ghosts[cell];

4014:       coarsePoint = quad->p.which_tree;
4015:       PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2);
4016:     }
4017:   }
4018:   for (t = flt; t <= llt; t++) {
4019:     p4est_tree_t     *tree     = &(trees[t]);
4020:     PetscInt          offset   = cLocalStart + tree->quadrants_offset;
4021:     PetscInt          numQuads = (PetscInt)tree->quadrants.elem_count;
4022:     p4est_quadrant_t *quads    = (p4est_quadrant_t *)tree->quadrants.array;
4023:     PetscInt          i;

4025:     if (!numQuads) continue;
4026:     coarsePoint = t;
4027:     for (i = 0; i < numQuads; i++, cp++) {
4028:       PetscInt          cell = i + offset;
4029:       p4est_quadrant_t *quad = &quads[i];

4031:       PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2);
4032:     }
4033:   }
4034:   if (cLocalEnd - cLocalStart < cEnd - cStart) {
4035:     p4est_quadrant_t *ghosts    = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4036:     PetscInt          numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4037:     PetscInt          count;

4039:     for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4040:       p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4041:       coarsePoint            = quad->p.which_tree;
4042:       PetscInt cell          = count + cLocalEnd;

4044:       PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2);
4045:     }
4046:   }
4047:   VecRestoreArray(cVec, &coords2);
4048:   DMSetCellCoordinatesLocal(plex, cVec);
4049:   VecDestroy(&cVec);
4050:   PetscSectionDestroy(&newSection);
4051:   return 0;
4052: }

4054:   #define DMForestClearAdaptivityForest_pforest _append_pforest(DMForestClearAdaptivityForest)
4055: static PetscErrorCode DMForestClearAdaptivityForest_pforest(DM dm)
4056: {
4057:   DM_Forest         *forest;
4058:   DM_Forest_pforest *pforest;

4060:   forest  = (DM_Forest *)dm->data;
4061:   pforest = (DM_Forest_pforest *)forest->data;
4062:   PetscSFDestroy(&(pforest->pointAdaptToSelfSF));
4063:   PetscSFDestroy(&(pforest->pointSelfToAdaptSF));
4064:   PetscFree(pforest->pointAdaptToSelfCids);
4065:   PetscFree(pforest->pointSelfToAdaptCids);
4066:   return 0;
4067: }

4069: static PetscErrorCode DMConvert_pforest_plex(DM dm, DMType newtype, DM *plex)
4070: {
4071:   DM_Forest           *forest;
4072:   DM_Forest_pforest   *pforest;
4073:   DM                   refTree, newPlex, base;
4074:   PetscInt             adjDim, adjCodim, coordDim;
4075:   MPI_Comm             comm;
4076:   PetscBool            isPforest;
4077:   PetscInt             dim;
4078:   PetscInt             overlap;
4079:   p4est_connect_type_t ctype;
4080:   p4est_locidx_t       first_local_quad = -1;
4081:   sc_array_t          *points_per_dim, *cone_sizes, *cones, *cone_orientations, *coords, *children, *parents, *childids, *leaves, *remotes;
4082:   PetscSection         parentSection;
4083:   PetscSF              pointSF;
4084:   size_t               zz, count;
4085:   PetscInt             pStart, pEnd;
4086:   DMLabel              ghostLabelBase = NULL;


4090:   comm = PetscObjectComm((PetscObject)dm);
4091:   PetscObjectTypeCompare((PetscObject)dm, DMPFOREST, &isPforest);
4093:   DMGetDimension(dm, &dim);
4095:   forest  = (DM_Forest *)dm->data;
4096:   pforest = (DM_Forest_pforest *)forest->data;
4097:   DMForestGetBaseDM(dm, &base);
4098:   if (base) DMGetLabel(base, "ghost", &ghostLabelBase);
4099:   if (!pforest->plex) {
4100:     PetscMPIInt size;

4102:     MPI_Comm_size(comm, &size);
4103:     DMCreate(comm, &newPlex);
4104:     DMSetType(newPlex, DMPLEX);
4105:     DMSetMatType(newPlex, dm->mattype);
4106:     /* share labels */
4107:     DMCopyLabels(dm, newPlex, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL);
4108:     DMForestGetAdjacencyDimension(dm, &adjDim);
4109:     DMForestGetAdjacencyCodimension(dm, &adjCodim);
4110:     DMGetCoordinateDim(dm, &coordDim);
4111:     if (adjDim == 0) {
4112:       ctype = P4EST_CONNECT_FULL;
4113:     } else if (adjCodim == 1) {
4114:       ctype = P4EST_CONNECT_FACE;
4115:   #if defined(P4_TO_P8)
4116:     } else if (adjDim == 1) {
4117:       ctype = P8EST_CONNECT_EDGE;
4118:   #endif
4119:     } else {
4120:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid adjacency dimension %" PetscInt_FMT, adjDim);
4121:     }
4123:     DMForestGetPartitionOverlap(dm, &overlap);
4124:     DMPlexSetOverlap_Plex(newPlex, NULL, overlap);

4126:     points_per_dim    = sc_array_new(sizeof(p4est_locidx_t));
4127:     cone_sizes        = sc_array_new(sizeof(p4est_locidx_t));
4128:     cones             = sc_array_new(sizeof(p4est_locidx_t));
4129:     cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
4130:     coords            = sc_array_new(3 * sizeof(double));
4131:     children          = sc_array_new(sizeof(p4est_locidx_t));
4132:     parents           = sc_array_new(sizeof(p4est_locidx_t));
4133:     childids          = sc_array_new(sizeof(p4est_locidx_t));
4134:     leaves            = sc_array_new(sizeof(p4est_locidx_t));
4135:     remotes           = sc_array_new(2 * sizeof(p4est_locidx_t));

4137:     PetscCallP4est(p4est_get_plex_data_ext, (pforest->forest, &pforest->ghost, &pforest->lnodes, ctype, (int)((size > 1) ? overlap : 0), &first_local_quad, points_per_dim, cone_sizes, cones, cone_orientations, coords, children, parents, childids, leaves, remotes, 1));

4139:     pforest->cLocalStart = (PetscInt)first_local_quad;
4140:     pforest->cLocalEnd   = pforest->cLocalStart + (PetscInt)pforest->forest->local_num_quadrants;
4141:     locidx_to_PetscInt(points_per_dim);
4142:     locidx_to_PetscInt(cone_sizes);
4143:     locidx_to_PetscInt(cones);
4144:     locidx_to_PetscInt(cone_orientations);
4145:     coords_double_to_PetscScalar(coords, coordDim);
4146:     locidx_to_PetscInt(children);
4147:     locidx_to_PetscInt(parents);
4148:     locidx_to_PetscInt(childids);
4149:     locidx_to_PetscInt(leaves);
4150:     locidx_pair_to_PetscSFNode(remotes);

4152:     DMSetDimension(newPlex, P4EST_DIM);
4153:     DMSetCoordinateDim(newPlex, coordDim);
4154:     DMPlexSetMaxProjectionHeight(newPlex, P4EST_DIM - 1);
4155:     DMPlexCreateFromDAG(newPlex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array);
4156:     DMPlexConvertOldOrientations_Internal(newPlex);
4157:     DMCreateReferenceTree_pforest(comm, &refTree);
4158:     DMPlexSetReferenceTree(newPlex, refTree);
4159:     PetscSectionCreate(comm, &parentSection);
4160:     DMPlexGetChart(newPlex, &pStart, &pEnd);
4161:     PetscSectionSetChart(parentSection, pStart, pEnd);
4162:     count = children->elem_count;
4163:     for (zz = 0; zz < count; zz++) {
4164:       PetscInt child = *((PetscInt *)sc_array_index(children, zz));

4166:       PetscSectionSetDof(parentSection, child, 1);
4167:     }
4168:     PetscSectionSetUp(parentSection);
4169:     DMPlexSetTree(newPlex, parentSection, (PetscInt *)parents->array, (PetscInt *)childids->array);
4170:     PetscSectionDestroy(&parentSection);
4171:     PetscSFCreate(comm, &pointSF);
4172:     /*
4173:        These arrays defining the sf are from the p4est library, but the code there shows the leaves being populated in increasing order.
4174:        https://gitlab.com/petsc/petsc/merge_requests/2248#note_240186391
4175:     */
4176:     PetscSFSetGraph(pointSF, pEnd - pStart, (PetscInt)leaves->elem_count, (PetscInt *)leaves->array, PETSC_COPY_VALUES, (PetscSFNode *)remotes->array, PETSC_COPY_VALUES);
4177:     DMSetPointSF(newPlex, pointSF);
4178:     DMSetPointSF(dm, pointSF);
4179:     {
4180:       DM coordDM;

4182:       DMGetCoordinateDM(newPlex, &coordDM);
4183:       DMSetPointSF(coordDM, pointSF);
4184:     }
4185:     PetscSFDestroy(&pointSF);
4186:     sc_array_destroy(points_per_dim);
4187:     sc_array_destroy(cone_sizes);
4188:     sc_array_destroy(cones);
4189:     sc_array_destroy(cone_orientations);
4190:     sc_array_destroy(coords);
4191:     sc_array_destroy(children);
4192:     sc_array_destroy(parents);
4193:     sc_array_destroy(childids);
4194:     sc_array_destroy(leaves);
4195:     sc_array_destroy(remotes);

4197:     {
4198:       const PetscReal *maxCell, *Lstart, *L;

4200:       DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
4201:       DMSetPeriodicity(newPlex, maxCell, Lstart, L);
4202:       DMPforestLocalizeCoordinates(dm, newPlex);
4203:     }

4205:     if (overlap > 0) { /* the p4est routine can't set all of the coordinates in its routine if there is overlap */
4206:       Vec                coordsGlobal, coordsLocal;
4207:       const PetscScalar *globalArray;
4208:       PetscScalar       *localArray;
4209:       PetscSF            coordSF;
4210:       DM                 coordDM;

4212:       DMGetCoordinateDM(newPlex, &coordDM);
4213:       DMGetSectionSF(coordDM, &coordSF);
4214:       DMGetCoordinates(newPlex, &coordsGlobal);
4215:       DMGetCoordinatesLocal(newPlex, &coordsLocal);
4216:       VecGetArrayRead(coordsGlobal, &globalArray);
4217:       VecGetArray(coordsLocal, &localArray);
4218:       PetscSFBcastBegin(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE);
4219:       PetscSFBcastEnd(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE);
4220:       VecRestoreArray(coordsLocal, &localArray);
4221:       VecRestoreArrayRead(coordsGlobal, &globalArray);
4222:       DMSetCoordinatesLocal(newPlex, coordsLocal);
4223:     }
4224:     DMPforestMapCoordinates(dm, newPlex);

4226:     pforest->plex = newPlex;

4228:     /* copy labels */
4229:     DMPforestLabelsFinalize(dm, newPlex);

4231:     if (ghostLabelBase || pforest->ghostName) { /* we have to do this after copying labels because the labels drive the construction of ghost cells */
4232:       PetscInt numAdded;
4233:       DM       newPlexGhosted;
4234:       void    *ctx;

4236:       DMPlexConstructGhostCells(newPlex, pforest->ghostName, &numAdded, &newPlexGhosted);
4237:       DMGetApplicationContext(newPlex, &ctx);
4238:       DMSetApplicationContext(newPlexGhosted, ctx);
4239:       /* we want the sf for the ghost dm to be the one for the p4est dm as well */
4240:       DMGetPointSF(newPlexGhosted, &pointSF);
4241:       DMSetPointSF(dm, pointSF);
4242:       DMDestroy(&newPlex);
4243:       DMPlexSetReferenceTree(newPlexGhosted, refTree);
4244:       DMForestClearAdaptivityForest_pforest(dm);
4245:       newPlex = newPlexGhosted;

4247:       /* share the labels back */
4248:       DMDestroyLabelLinkList_Internal(dm);
4249:       DMCopyLabels(newPlex, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL);
4250:       pforest->plex = newPlex;
4251:     }
4252:     DMDestroy(&refTree);
4253:     if (dm->setfromoptionscalled) {
4254:       PetscObjectOptionsBegin((PetscObject)newPlex);
4255:       DMSetFromOptions_NonRefinement_Plex(newPlex, PetscOptionsObject);
4256:       PetscObjectProcessOptionsHandlers((PetscObject)newPlex, PetscOptionsObject);
4257:       PetscOptionsEnd();
4258:     }
4259:     DMViewFromOptions(newPlex, NULL, "-dm_p4est_plex_view");
4260:     {
4261:       DM           cdm;
4262:       PetscSection coordsSec;
4263:       Vec          coords;
4264:       PetscInt     cDim;

4266:       DMGetCoordinateDim(newPlex, &cDim);
4267:       DMGetCoordinateSection(newPlex, &coordsSec);
4268:       DMSetCoordinateSection(dm, cDim, coordsSec);
4269:       DMGetCoordinatesLocal(newPlex, &coords);
4270:       DMSetCoordinatesLocal(dm, coords);
4271:       DMGetCellCoordinateDM(newPlex, &cdm);
4272:       if (cdm) DMSetCellCoordinateDM(dm, cdm);
4273:       DMGetCellCoordinateSection(newPlex, &coordsSec);
4274:       if (coordsSec) DMSetCellCoordinateSection(dm, cDim, coordsSec);
4275:       DMGetCellCoordinatesLocal(newPlex, &coords);
4276:       if (coords) DMSetCellCoordinatesLocal(dm, coords);
4277:     }
4278:   }
4279:   newPlex = pforest->plex;
4280:   if (plex) {
4281:     DMClone(newPlex, plex);
4282:   #if 0
4283:     DMGetCoordinateDM(newPlex,&coordDM);
4284:     DMSetCoordinateDM(*plex,coordDM);
4285:     DMGetCellCoordinateDM(newPlex,&coordDM);
4286:     DMSetCellCoordinateDM(*plex,coordDM);
4287:   #endif
4288:     DMShareDiscretization(dm, *plex);
4289:   }
4290:   return 0;
4291: }

4293: static PetscErrorCode DMSetFromOptions_pforest(DM dm, PetscOptionItems *PetscOptionsObject)
4294: {
4295:   DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4296:   char               stringBuffer[256];
4297:   PetscBool          flg;

4299:   DMSetFromOptions_Forest(dm, PetscOptionsObject);
4300:   PetscOptionsHeadBegin(PetscOptionsObject, "DM" P4EST_STRING " options");
4301:   PetscOptionsBool("-dm_p4est_partition_for_coarsening", "partition forest to allow for coarsening", "DMP4estSetPartitionForCoarsening", pforest->partition_for_coarsening, &(pforest->partition_for_coarsening), NULL);
4302:   PetscOptionsString("-dm_p4est_ghost_label_name", "the name of the ghost label when converting from a DMPlex", NULL, NULL, stringBuffer, sizeof(stringBuffer), &flg);
4303:   PetscOptionsHeadEnd();
4304:   if (flg) {
4305:     PetscFree(pforest->ghostName);
4306:     PetscStrallocpy(stringBuffer, &pforest->ghostName);
4307:   }
4308:   return 0;
4309: }

4311:   #if !defined(P4_TO_P8)
4312:     #define DMPforestGetPartitionForCoarsening DMP4estGetPartitionForCoarsening
4313:     #define DMPforestSetPartitionForCoarsening DMP4estSetPartitionForCoarsening
4314:   #else
4315:     #define DMPforestGetPartitionForCoarsening DMP8estGetPartitionForCoarsening
4316:     #define DMPforestSetPartitionForCoarsening DMP8estSetPartitionForCoarsening
4317:   #endif

4319: PETSC_EXTERN PetscErrorCode DMPforestGetPartitionForCoarsening(DM dm, PetscBool *flg)
4320: {
4321:   DM_Forest_pforest *pforest;

4324:   pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4325:   *flg    = pforest->partition_for_coarsening;
4326:   return 0;
4327: }

4329: PETSC_EXTERN PetscErrorCode DMPforestSetPartitionForCoarsening(DM dm, PetscBool flg)
4330: {
4331:   DM_Forest_pforest *pforest;

4334:   pforest                           = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4335:   pforest->partition_for_coarsening = flg;
4336:   return 0;
4337: }

4339: static PetscErrorCode DMPforestGetPlex(DM dm, DM *plex)
4340: {
4341:   DM_Forest_pforest *pforest;

4343:   if (plex) *plex = NULL;
4344:   DMSetUp(dm);
4345:   pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4346:   if (!pforest->plex) DMConvert_pforest_plex(dm, DMPLEX, NULL);
4347:   DMShareDiscretization(dm, pforest->plex);
4348:   if (plex) *plex = pforest->plex;
4349:   return 0;
4350: }

4352:   #define DMCreateInterpolation_pforest _append_pforest(DMCreateInterpolation)
4353: static PetscErrorCode DMCreateInterpolation_pforest(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
4354: {
4355:   PetscSection gsc, gsf;
4356:   PetscInt     m, n;
4357:   DM           cdm;

4359:   DMGetGlobalSection(dmFine, &gsf);
4360:   PetscSectionGetConstrainedStorageSize(gsf, &m);
4361:   DMGetGlobalSection(dmCoarse, &gsc);
4362:   PetscSectionGetConstrainedStorageSize(gsc, &n);

4364:   MatCreate(PetscObjectComm((PetscObject)dmFine), interpolation);
4365:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
4366:   MatSetType(*interpolation, MATAIJ);

4368:   DMGetCoarseDM(dmFine, &cdm);

4371:   {
4372:     DM        plexF, plexC;
4373:     PetscSF   sf;
4374:     PetscInt *cids;
4375:     PetscInt  dofPerDim[4] = {1, 1, 1, 1};

4377:     DMPforestGetPlex(dmCoarse, &plexC);
4378:     DMPforestGetPlex(dmFine, &plexF);
4379:     DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids);
4380:     PetscSFSetUp(sf);
4381:     DMPlexComputeInterpolatorTree(plexC, plexF, sf, cids, *interpolation);
4382:     PetscSFDestroy(&sf);
4383:     PetscFree(cids);
4384:   }
4385:   MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
4386:   /* Use naive scaling */
4387:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
4388:   return 0;
4389: }

4391:   #define DMCreateInjection_pforest _append_pforest(DMCreateInjection)
4392: static PetscErrorCode DMCreateInjection_pforest(DM dmCoarse, DM dmFine, Mat *injection)
4393: {
4394:   PetscSection gsc, gsf;
4395:   PetscInt     m, n;
4396:   DM           cdm;

4398:   DMGetGlobalSection(dmFine, &gsf);
4399:   PetscSectionGetConstrainedStorageSize(gsf, &n);
4400:   DMGetGlobalSection(dmCoarse, &gsc);
4401:   PetscSectionGetConstrainedStorageSize(gsc, &m);

4403:   MatCreate(PetscObjectComm((PetscObject)dmFine), injection);
4404:   MatSetSizes(*injection, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
4405:   MatSetType(*injection, MATAIJ);

4407:   DMGetCoarseDM(dmFine, &cdm);

4410:   {
4411:     DM        plexF, plexC;
4412:     PetscSF   sf;
4413:     PetscInt *cids;
4414:     PetscInt  dofPerDim[4] = {1, 1, 1, 1};

4416:     DMPforestGetPlex(dmCoarse, &plexC);
4417:     DMPforestGetPlex(dmFine, &plexF);
4418:     DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids);
4419:     PetscSFSetUp(sf);
4420:     DMPlexComputeInjectorTree(plexC, plexF, sf, cids, *injection);
4421:     PetscSFDestroy(&sf);
4422:     PetscFree(cids);
4423:   }
4424:   MatViewFromOptions(*injection, NULL, "-inject_mat_view");
4425:   /* Use naive scaling */
4426:   return 0;
4427: }

4429:   #define DMForestTransferVecFromBase_pforest _append_pforest(DMForestTransferVecFromBase)
4430: static PetscErrorCode DMForestTransferVecFromBase_pforest(DM dm, Vec vecIn, Vec vecOut)
4431: {
4432:   DM        dmIn, dmVecIn, base, basec, plex, coarseDM;
4433:   DM       *hierarchy;
4434:   PetscSF   sfRed = NULL;
4435:   PetscDS   ds;
4436:   Vec       vecInLocal, vecOutLocal;
4437:   DMLabel   subpointMap;
4438:   PetscInt  minLevel, mh, n_hi, i;
4439:   PetscBool hiforest, *hierarchy_forest;

4441:   VecGetDM(vecIn, &dmVecIn);
4442:   DMGetDS(dmVecIn, &ds);
4444:   { /* we cannot stick user contexts into function callbacks for DMProjectFieldLocal! */
4445:     PetscSection section;
4446:     PetscInt     Nf;

4448:     DMGetLocalSection(dmVecIn, &section);
4449:     PetscSectionGetNumFields(section, &Nf);
4451:   }
4452:   DMForestGetMinimumRefinement(dm, &minLevel);
4454:   DMForestGetBaseDM(dm, &base);

4457:   VecSet(vecOut, 0.0);
4458:   if (dmVecIn == base) { /* sequential runs */
4459:     PetscObjectReference((PetscObject)vecIn);
4460:   } else {
4461:     PetscSection secIn, secInRed;
4462:     Vec          vecInRed, vecInLocal;

4464:     PetscObjectQuery((PetscObject)base, "_base_migration_sf", (PetscObject *)&sfRed);
4466:     PetscSectionCreate(PetscObjectComm((PetscObject)dmVecIn), &secInRed);
4467:     VecCreate(PETSC_COMM_SELF, &vecInRed);
4468:     DMGetLocalSection(dmVecIn, &secIn);
4469:     DMGetLocalVector(dmVecIn, &vecInLocal);
4470:     DMGlobalToLocalBegin(dmVecIn, vecIn, INSERT_VALUES, vecInLocal);
4471:     DMGlobalToLocalEnd(dmVecIn, vecIn, INSERT_VALUES, vecInLocal);
4472:     DMPlexDistributeField(dmVecIn, sfRed, secIn, vecInLocal, secInRed, vecInRed);
4473:     DMRestoreLocalVector(dmVecIn, &vecInLocal);
4474:     PetscSectionDestroy(&secInRed);
4475:     vecIn = vecInRed;
4476:   }

4478:   /* we first search through the AdaptivityForest hierarchy
4479:      once we found the first disconnected forest, we upsweep the DM hierarchy */
4480:   hiforest = PETSC_TRUE;

4482:   /* upsweep to the coarsest DM */
4483:   n_hi     = 0;
4484:   coarseDM = dm;
4485:   do {
4486:     PetscBool isforest;

4488:     dmIn = coarseDM;
4489:     /* need to call DMSetUp to have the hierarchy recursively setup */
4490:     DMSetUp(dmIn);
4491:     DMIsForest(dmIn, &isforest);
4493:     coarseDM = NULL;
4494:     if (hiforest) DMForestGetAdaptivityForest(dmIn, &coarseDM);
4495:     if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4496:       hiforest = PETSC_FALSE;
4497:       DMGetCoarseDM(dmIn, &coarseDM);
4498:     }
4499:     n_hi++;
4500:   } while (coarseDM);

4502:   PetscMalloc2(n_hi, &hierarchy, n_hi, &hierarchy_forest);

4504:   i        = 0;
4505:   hiforest = PETSC_TRUE;
4506:   coarseDM = dm;
4507:   do {
4508:     dmIn     = coarseDM;
4509:     coarseDM = NULL;
4510:     if (hiforest) DMForestGetAdaptivityForest(dmIn, &coarseDM);
4511:     if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4512:       hiforest = PETSC_FALSE;
4513:       DMGetCoarseDM(dmIn, &coarseDM);
4514:     }
4515:     i++;
4516:     hierarchy[n_hi - i] = dmIn;
4517:   } while (coarseDM);

4519:   /* project base vector on the coarsest forest (minimum refinement = 0) */
4520:   DMPforestGetPlex(dmIn, &plex);

4522:   /* Check this plex is compatible with the base */
4523:   {
4524:     IS       gnum[2];
4525:     PetscInt ncells[2], gncells[2];

4527:     DMPlexGetCellNumbering(base, &gnum[0]);
4528:     DMPlexGetCellNumbering(plex, &gnum[1]);
4529:     ISGetMinMax(gnum[0], NULL, &ncells[0]);
4530:     ISGetMinMax(gnum[1], NULL, &ncells[1]);
4531:     MPIU_Allreduce(ncells, gncells, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
4533:   }

4535:   DMGetLabel(dmIn, "_forest_base_subpoint_map", &subpointMap);

4538:   DMPlexGetMaxProjectionHeight(base, &mh);
4539:   DMPlexSetMaxProjectionHeight(plex, mh);

4541:   DMClone(base, &basec);
4542:   DMCopyDisc(dmVecIn, basec);
4543:   if (sfRed) {
4544:     PetscObjectReference((PetscObject)vecIn);
4545:     vecInLocal = vecIn;
4546:   } else {
4547:     DMCreateLocalVector(basec, &vecInLocal);
4548:     DMGlobalToLocalBegin(basec, vecIn, INSERT_VALUES, vecInLocal);
4549:     DMGlobalToLocalEnd(basec, vecIn, INSERT_VALUES, vecInLocal);
4550:   }

4552:   DMGetLocalVector(dmIn, &vecOutLocal);
4553:   { /* get degrees of freedom ordered onto dmIn */
4554:     PetscSF            basetocoarse;
4555:     PetscInt           bStart, bEnd, nroots;
4556:     PetscInt           iStart, iEnd, nleaves, leaf;
4557:     PetscMPIInt        rank;
4558:     PetscSFNode       *remotes;
4559:     PetscSection       secIn, secOut;
4560:     PetscInt          *remoteOffsets;
4561:     PetscSF            transferSF;
4562:     const PetscScalar *inArray;
4563:     PetscScalar       *outArray;

4565:     MPI_Comm_rank(PetscObjectComm((PetscObject)basec), &rank);
4566:     DMPlexGetChart(basec, &bStart, &bEnd);
4567:     nroots = PetscMax(bEnd - bStart, 0);
4568:     DMPlexGetChart(plex, &iStart, &iEnd);
4569:     nleaves = PetscMax(iEnd - iStart, 0);

4571:     PetscMalloc1(nleaves, &remotes);
4572:     for (leaf = iStart; leaf < iEnd; leaf++) {
4573:       PetscInt index;

4575:       remotes[leaf - iStart].rank = rank;
4576:       DMLabelGetValue(subpointMap, leaf, &index);
4577:       remotes[leaf - iStart].index = index;
4578:     }

4580:     PetscSFCreate(PetscObjectComm((PetscObject)basec), &basetocoarse);
4581:     PetscSFSetGraph(basetocoarse, nroots, nleaves, NULL, PETSC_OWN_POINTER, remotes, PETSC_OWN_POINTER);
4582:     PetscSFSetUp(basetocoarse);
4583:     DMGetLocalSection(basec, &secIn);
4584:     PetscSectionCreate(PetscObjectComm((PetscObject)dmIn), &secOut);
4585:     PetscSFDistributeSection(basetocoarse, secIn, &remoteOffsets, secOut);
4586:     PetscSFCreateSectionSF(basetocoarse, secIn, remoteOffsets, secOut, &transferSF);
4587:     PetscFree(remoteOffsets);
4588:     VecGetArrayWrite(vecOutLocal, &outArray);
4589:     VecGetArrayRead(vecInLocal, &inArray);
4590:     PetscSFBcastBegin(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE);
4591:     PetscSFBcastEnd(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE);
4592:     VecRestoreArrayRead(vecInLocal, &inArray);
4593:     VecRestoreArrayWrite(vecOutLocal, &outArray);
4594:     PetscSFDestroy(&transferSF);
4595:     PetscSectionDestroy(&secOut);
4596:     PetscSFDestroy(&basetocoarse);
4597:   }
4598:   VecDestroy(&vecInLocal);
4599:   DMDestroy(&basec);
4600:   VecDestroy(&vecIn);

4602:   /* output */
4603:   if (n_hi > 1) { /* downsweep the stored hierarchy */
4604:     Vec vecOut1, vecOut2;
4605:     DM  fineDM;

4607:     DMGetGlobalVector(dmIn, &vecOut1);
4608:     DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut1);
4609:     DMRestoreLocalVector(dmIn, &vecOutLocal);
4610:     for (i = 1; i < n_hi - 1; i++) {
4611:       fineDM = hierarchy[i];
4612:       DMGetGlobalVector(fineDM, &vecOut2);
4613:       DMForestTransferVec(dmIn, vecOut1, fineDM, vecOut2, PETSC_TRUE, 0.0);
4614:       DMRestoreGlobalVector(dmIn, &vecOut1);
4615:       vecOut1 = vecOut2;
4616:       dmIn    = fineDM;
4617:     }
4618:     DMForestTransferVec(dmIn, vecOut1, dm, vecOut, PETSC_TRUE, 0.0);
4619:     DMRestoreGlobalVector(dmIn, &vecOut1);
4620:   } else {
4621:     DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut);
4622:     DMRestoreLocalVector(dmIn, &vecOutLocal);
4623:   }
4624:   PetscFree2(hierarchy, hierarchy_forest);
4625:   return 0;
4626: }

4628:   #define DMForestTransferVec_pforest _append_pforest(DMForestTransferVec)
4629: static PetscErrorCode DMForestTransferVec_pforest(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscBool useBCs, PetscReal time)
4630: {
4631:   DM          adaptIn, adaptOut, plexIn, plexOut;
4632:   DM_Forest  *forestIn, *forestOut, *forestAdaptIn, *forestAdaptOut;
4633:   PetscInt    dofPerDim[] = {1, 1, 1, 1};
4634:   PetscSF     inSF = NULL, outSF = NULL;
4635:   PetscInt   *inCids = NULL, *outCids = NULL;
4636:   DMAdaptFlag purposeIn, purposeOut;

4638:   forestOut = (DM_Forest *)dmOut->data;
4639:   forestIn  = (DM_Forest *)dmIn->data;

4641:   DMForestGetAdaptivityForest(dmOut, &adaptOut);
4642:   DMForestGetAdaptivityPurpose(dmOut, &purposeOut);
4643:   forestAdaptOut = adaptOut ? (DM_Forest *)adaptOut->data : NULL;

4645:   DMForestGetAdaptivityForest(dmIn, &adaptIn);
4646:   DMForestGetAdaptivityPurpose(dmIn, &purposeIn);
4647:   forestAdaptIn = adaptIn ? (DM_Forest *)adaptIn->data : NULL;

4649:   if (forestAdaptOut == forestIn) {
4650:     switch (purposeOut) {
4651:     case DM_ADAPT_REFINE:
4652:       DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids);
4653:       PetscSFSetUp(inSF);
4654:       break;
4655:     case DM_ADAPT_COARSEN:
4656:     case DM_ADAPT_COARSEN_LAST:
4657:       DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &outCids);
4658:       PetscSFSetUp(outSF);
4659:       break;
4660:     default:
4661:       DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids);
4662:       DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids);
4663:       PetscSFSetUp(inSF);
4664:       PetscSFSetUp(outSF);
4665:     }
4666:   } else if (forestAdaptIn == forestOut) {
4667:     switch (purposeIn) {
4668:     case DM_ADAPT_REFINE:
4669:       DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &inCids);
4670:       PetscSFSetUp(outSF);
4671:       break;
4672:     case DM_ADAPT_COARSEN:
4673:     case DM_ADAPT_COARSEN_LAST:
4674:       DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids);
4675:       PetscSFSetUp(inSF);
4676:       break;
4677:     default:
4678:       DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids);
4679:       DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids);
4680:       PetscSFSetUp(inSF);
4681:       PetscSFSetUp(outSF);
4682:     }
4683:   } else SETERRQ(PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Only support transfer from pre-adaptivity to post-adaptivity right now");
4684:   DMPforestGetPlex(dmIn, &plexIn);
4685:   DMPforestGetPlex(dmOut, &plexOut);

4687:   DMPlexTransferVecTree(plexIn, vecIn, plexOut, vecOut, inSF, outSF, inCids, outCids, useBCs, time);
4688:   PetscFree(inCids);
4689:   PetscFree(outCids);
4690:   PetscSFDestroy(&inSF);
4691:   PetscSFDestroy(&outSF);
4692:   PetscFree(inCids);
4693:   PetscFree(outCids);
4694:   return 0;
4695: }

4697:   #define DMCreateCoordinateDM_pforest _append_pforest(DMCreateCoordinateDM)
4698: static PetscErrorCode DMCreateCoordinateDM_pforest(DM dm, DM *cdm)
4699: {
4700:   DM plex;

4703:   DMPforestGetPlex(dm, &plex);
4704:   DMGetCoordinateDM(plex, cdm);
4705:   PetscObjectReference((PetscObject)*cdm);
4706:   return 0;
4707: }

4709:   #define VecViewLocal_pforest _append_pforest(VecViewLocal)
4710: static PetscErrorCode VecViewLocal_pforest(Vec vec, PetscViewer viewer)
4711: {
4712:   DM dm, plex;

4714:   VecGetDM(vec, &dm);
4715:   DMPforestGetPlex(dm, &plex);
4716:   VecSetDM(vec, plex);
4717:   VecView_Plex_Local(vec, viewer);
4718:   VecSetDM(vec, dm);
4719:   return 0;
4720: }

4722:   #define VecView_pforest _append_pforest(VecView)
4723: static PetscErrorCode VecView_pforest(Vec vec, PetscViewer viewer)
4724: {
4725:   DM dm, plex;

4727:   VecGetDM(vec, &dm);
4728:   DMPforestGetPlex(dm, &plex);
4729:   VecSetDM(vec, plex);
4730:   VecView_Plex(vec, viewer);
4731:   VecSetDM(vec, dm);
4732:   return 0;
4733: }

4735:   #define VecView_pforest_Native _infix_pforest(VecView, _Native)
4736: static PetscErrorCode VecView_pforest_Native(Vec vec, PetscViewer viewer)
4737: {
4738:   DM dm, plex;

4740:   VecGetDM(vec, &dm);
4741:   DMPforestGetPlex(dm, &plex);
4742:   VecSetDM(vec, plex);
4743:   VecView_Plex_Native(vec, viewer);
4744:   VecSetDM(vec, dm);
4745:   return 0;
4746: }

4748:   #define VecLoad_pforest _append_pforest(VecLoad)
4749: static PetscErrorCode VecLoad_pforest(Vec vec, PetscViewer viewer)
4750: {
4751:   DM dm, plex;

4753:   VecGetDM(vec, &dm);
4754:   DMPforestGetPlex(dm, &plex);
4755:   VecSetDM(vec, plex);
4756:   VecLoad_Plex(vec, viewer);
4757:   VecSetDM(vec, dm);
4758:   return 0;
4759: }

4761:   #define VecLoad_pforest_Native _infix_pforest(VecLoad, _Native)
4762: static PetscErrorCode VecLoad_pforest_Native(Vec vec, PetscViewer viewer)
4763: {
4764:   DM dm, plex;

4766:   VecGetDM(vec, &dm);
4767:   DMPforestGetPlex(dm, &plex);
4768:   VecSetDM(vec, plex);
4769:   VecLoad_Plex_Native(vec, viewer);
4770:   VecSetDM(vec, dm);
4771:   return 0;
4772: }

4774:   #define DMCreateGlobalVector_pforest _append_pforest(DMCreateGlobalVector)
4775: static PetscErrorCode DMCreateGlobalVector_pforest(DM dm, Vec *vec)
4776: {
4777:   DMCreateGlobalVector_Section_Private(dm, vec);
4778:   /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
4779:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_pforest);
4780:   VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void))VecView_pforest_Native);
4781:   VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_pforest);
4782:   VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void))VecLoad_pforest_Native);
4783:   return 0;
4784: }

4786:   #define DMCreateLocalVector_pforest _append_pforest(DMCreateLocalVector)
4787: static PetscErrorCode DMCreateLocalVector_pforest(DM dm, Vec *vec)
4788: {
4789:   DMCreateLocalVector_Section_Private(dm, vec);
4790:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecViewLocal_pforest);
4791:   return 0;
4792: }

4794:   #define DMCreateMatrix_pforest _append_pforest(DMCreateMatrix)
4795: static PetscErrorCode DMCreateMatrix_pforest(DM dm, Mat *mat)
4796: {
4797:   DM plex;

4800:   DMPforestGetPlex(dm, &plex);
4801:   if (plex->prealloc_only != dm->prealloc_only) plex->prealloc_only = dm->prealloc_only; /* maybe this should go into forest->plex */
4802:   DMCreateMatrix(plex, mat);
4803:   MatSetDM(*mat, dm);
4804:   return 0;
4805: }

4807:   #define DMProjectFunctionLocal_pforest _append_pforest(DMProjectFunctionLocal)
4808: static PetscErrorCode DMProjectFunctionLocal_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4809: {
4810:   DM plex;

4813:   DMPforestGetPlex(dm, &plex);
4814:   DMProjectFunctionLocal(plex, time, funcs, ctxs, mode, localX);
4815:   return 0;
4816: }

4818:   #define DMProjectFunctionLabelLocal_pforest _append_pforest(DMProjectFunctionLabelLocal)
4819: static PetscErrorCode DMProjectFunctionLabelLocal_pforest(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4820: {
4821:   DM plex;

4824:   DMPforestGetPlex(dm, &plex);
4825:   DMProjectFunctionLabelLocal(plex, time, label, numIds, ids, Ncc, comps, funcs, ctxs, mode, localX);
4826:   return 0;
4827: }

4829:   #define DMProjectFieldLocal_pforest _append_pforest(DMProjectFieldLocal)
4830: PetscErrorCode DMProjectFieldLocal_pforest(DM dm, PetscReal time, Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
4831: {
4832:   DM plex;

4835:   DMPforestGetPlex(dm, &plex);
4836:   DMProjectFieldLocal(plex, time, localU, funcs, mode, localX);
4837:   return 0;
4838: }

4840:   #define DMComputeL2Diff_pforest _append_pforest(DMComputeL2Diff)
4841: PetscErrorCode DMComputeL2Diff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
4842: {
4843:   DM plex;

4846:   DMPforestGetPlex(dm, &plex);
4847:   DMComputeL2Diff(plex, time, funcs, ctxs, X, diff);
4848:   return 0;
4849: }

4851:   #define DMComputeL2FieldDiff_pforest _append_pforest(DMComputeL2FieldDiff)
4852: PetscErrorCode DMComputeL2FieldDiff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
4853: {
4854:   DM plex;

4857:   DMPforestGetPlex(dm, &plex);
4858:   DMComputeL2FieldDiff(plex, time, funcs, ctxs, X, diff);
4859:   return 0;
4860: }

4862:   #define DMCreatelocalsection_pforest _append_pforest(DMCreatelocalsection)
4863: static PetscErrorCode DMCreatelocalsection_pforest(DM dm)
4864: {
4865:   DM           plex;
4866:   PetscSection section;

4869:   DMPforestGetPlex(dm, &plex);
4870:   DMGetLocalSection(plex, &section);
4871:   DMSetLocalSection(dm, section);
4872:   return 0;
4873: }

4875:   #define DMCreateDefaultConstraints_pforest _append_pforest(DMCreateDefaultConstraints)
4876: static PetscErrorCode DMCreateDefaultConstraints_pforest(DM dm)
4877: {
4878:   DM           plex;
4879:   Mat          mat;
4880:   Vec          bias;
4881:   PetscSection section;

4884:   DMPforestGetPlex(dm, &plex);
4885:   DMGetDefaultConstraints(plex, &section, &mat, &bias);
4886:   DMSetDefaultConstraints(dm, section, mat, bias);
4887:   return 0;
4888: }

4890:   #define DMGetDimPoints_pforest _append_pforest(DMGetDimPoints)
4891: static PetscErrorCode DMGetDimPoints_pforest(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
4892: {
4893:   DM plex;

4896:   DMPforestGetPlex(dm, &plex);
4897:   DMGetDimPoints(plex, dim, cStart, cEnd);
4898:   return 0;
4899: }

4901:   /* Need to forward declare */
4902:   #define DMInitialize_pforest _append_pforest(DMInitialize)
4903: static PetscErrorCode DMInitialize_pforest(DM dm);

4905:   #define DMClone_pforest _append_pforest(DMClone)
4906: static PetscErrorCode DMClone_pforest(DM dm, DM *newdm)
4907: {
4908:   DMClone_Forest(dm, newdm);
4909:   DMInitialize_pforest(*newdm);
4910:   return 0;
4911: }

4913:   #define DMForestCreateCellChart_pforest _append_pforest(DMForestCreateCellChart)
4914: static PetscErrorCode DMForestCreateCellChart_pforest(DM dm, PetscInt *cStart, PetscInt *cEnd)
4915: {
4916:   DM_Forest         *forest;
4917:   DM_Forest_pforest *pforest;
4918:   PetscInt           overlap;

4920:   DMSetUp(dm);
4921:   forest  = (DM_Forest *)dm->data;
4922:   pforest = (DM_Forest_pforest *)forest->data;
4923:   *cStart = 0;
4924:   DMForestGetPartitionOverlap(dm, &overlap);
4925:   if (overlap && pforest->ghost) {
4926:     *cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[pforest->forest->mpisize];
4927:   } else {
4928:     *cEnd = pforest->forest->local_num_quadrants;
4929:   }
4930:   return 0;
4931: }

4933:   #define DMForestCreateCellSF_pforest _append_pforest(DMForestCreateCellSF)
4934: static PetscErrorCode DMForestCreateCellSF_pforest(DM dm, PetscSF *cellSF)
4935: {
4936:   DM_Forest         *forest;
4937:   DM_Forest_pforest *pforest;
4938:   PetscMPIInt        rank;
4939:   PetscInt           overlap;
4940:   PetscInt           cStart, cEnd, cLocalStart, cLocalEnd;
4941:   PetscInt           nRoots, nLeaves, *mine = NULL;
4942:   PetscSFNode       *remote = NULL;
4943:   PetscSF            sf;

4945:   DMForestGetCellChart(dm, &cStart, &cEnd);
4946:   forest      = (DM_Forest *)dm->data;
4947:   pforest     = (DM_Forest_pforest *)forest->data;
4948:   nRoots      = cEnd - cStart;
4949:   cLocalStart = pforest->cLocalStart;
4950:   cLocalEnd   = pforest->cLocalEnd;
4951:   nLeaves     = 0;
4952:   DMForestGetPartitionOverlap(dm, &overlap);
4953:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
4954:   if (overlap && pforest->ghost) {
4955:     PetscSFNode      *mirror;
4956:     p4est_quadrant_t *mirror_array;
4957:     PetscInt          nMirror, nGhostPre, nSelf, q;
4958:     void            **mirrorPtrs;

4960:     nMirror   = (PetscInt)pforest->ghost->mirrors.elem_count;
4961:     nSelf     = cLocalEnd - cLocalStart;
4962:     nLeaves   = nRoots - nSelf;
4963:     nGhostPre = (PetscInt)pforest->ghost->proc_offsets[rank];
4964:     PetscMalloc1(nLeaves, &mine);
4965:     PetscMalloc1(nLeaves, &remote);
4966:     PetscMalloc2(nMirror, &mirror, nMirror, &mirrorPtrs);
4967:     mirror_array = (p4est_quadrant_t *)pforest->ghost->mirrors.array;
4968:     for (q = 0; q < nMirror; q++) {
4969:       p4est_quadrant_t *mir = &(mirror_array[q]);

4971:       mirror[q].rank  = rank;
4972:       mirror[q].index = (PetscInt)mir->p.piggy3.local_num + cLocalStart;
4973:       mirrorPtrs[q]   = (void *)&(mirror[q]);
4974:     }
4975:     PetscCallP4est(p4est_ghost_exchange_custom, (pforest->forest, pforest->ghost, sizeof(PetscSFNode), mirrorPtrs, remote));
4976:     PetscFree2(mirror, mirrorPtrs);
4977:     for (q = 0; q < nGhostPre; q++) mine[q] = q;
4978:     for (; q < nLeaves; q++) mine[q] = (q - nGhostPre) + cLocalEnd;
4979:   }
4980:   PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
4981:   PetscSFSetGraph(sf, nRoots, nLeaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4982:   *cellSF = sf;
4983:   return 0;
4984: }

4986: static PetscErrorCode DMCreateNeumannOverlap_pforest(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx)
4987: {
4988:   DM plex;

4990:   DMPforestGetPlex(dm, &plex);
4991:   DMCreateNeumannOverlap_Plex(plex, ovl, J, setup, setup_ctx);
4992:   if (!*setup) {
4993:     PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup);
4994:     if (*setup) PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm);
4995:   }
4996:   return 0;
4997: }

4999: static PetscErrorCode DMInitialize_pforest(DM dm)
5000: {
5001:   dm->ops->setup                     = DMSetUp_pforest;
5002:   dm->ops->view                      = DMView_pforest;
5003:   dm->ops->clone                     = DMClone_pforest;
5004:   dm->ops->createinterpolation       = DMCreateInterpolation_pforest;
5005:   dm->ops->createinjection           = DMCreateInjection_pforest;
5006:   dm->ops->setfromoptions            = DMSetFromOptions_pforest;
5007:   dm->ops->createcoordinatedm        = DMCreateCoordinateDM_pforest;
5008:   dm->ops->createglobalvector        = DMCreateGlobalVector_pforest;
5009:   dm->ops->createlocalvector         = DMCreateLocalVector_pforest;
5010:   dm->ops->creatematrix              = DMCreateMatrix_pforest;
5011:   dm->ops->projectfunctionlocal      = DMProjectFunctionLocal_pforest;
5012:   dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_pforest;
5013:   dm->ops->projectfieldlocal         = DMProjectFieldLocal_pforest;
5014:   dm->ops->createlocalsection        = DMCreatelocalsection_pforest;
5015:   dm->ops->createdefaultconstraints  = DMCreateDefaultConstraints_pforest;
5016:   dm->ops->computel2diff             = DMComputeL2Diff_pforest;
5017:   dm->ops->computel2fielddiff        = DMComputeL2FieldDiff_pforest;
5018:   dm->ops->getdimpoints              = DMGetDimPoints_pforest;

5020:   PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", DMConvert_plex_pforest);
5021:   PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", DMConvert_pforest_plex);
5022:   PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_pforest);
5023:   PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMForestGetPartitionOverlap);
5024:   return 0;
5025: }

5027:   #define DMCreate_pforest _append_pforest(DMCreate)
5028: PETSC_EXTERN PetscErrorCode DMCreate_pforest(DM dm)
5029: {
5030:   DM_Forest         *forest;
5031:   DM_Forest_pforest *pforest;

5033:   PetscP4estInitialize();
5034:   DMCreate_Forest(dm);
5035:   DMInitialize_pforest(dm);
5036:   DMSetDimension(dm, P4EST_DIM);

5038:   /* set forest defaults */
5039:   DMForestSetTopology(dm, "unit");
5040:   DMForestSetMinimumRefinement(dm, 0);
5041:   DMForestSetInitialRefinement(dm, 0);
5042:   DMForestSetMaximumRefinement(dm, P4EST_QMAXLEVEL);
5043:   DMForestSetGradeFactor(dm, 2);
5044:   DMForestSetAdjacencyDimension(dm, 0);
5045:   DMForestSetPartitionOverlap(dm, 0);

5047:   /* create p4est data */
5048:   PetscNew(&pforest);

5050:   forest                            = (DM_Forest *)dm->data;
5051:   forest->data                      = pforest;
5052:   forest->destroy                   = DMForestDestroy_pforest;
5053:   forest->ftemplate                 = DMForestTemplate_pforest;
5054:   forest->transfervec               = DMForestTransferVec_pforest;
5055:   forest->transfervecfrombase       = DMForestTransferVecFromBase_pforest;
5056:   forest->createcellchart           = DMForestCreateCellChart_pforest;
5057:   forest->createcellsf              = DMForestCreateCellSF_pforest;
5058:   forest->clearadaptivityforest     = DMForestClearAdaptivityForest_pforest;
5059:   forest->getadaptivitysuccess      = DMForestGetAdaptivitySuccess_pforest;
5060:   pforest->topo                     = NULL;
5061:   pforest->forest                   = NULL;
5062:   pforest->ghost                    = NULL;
5063:   pforest->lnodes                   = NULL;
5064:   pforest->partition_for_coarsening = PETSC_TRUE;
5065:   pforest->coarsen_hierarchy        = PETSC_FALSE;
5066:   pforest->cLocalStart              = -1;
5067:   pforest->cLocalEnd                = -1;
5068:   pforest->labelsFinalized          = PETSC_FALSE;
5069:   pforest->ghostName                = NULL;
5070:   return 0;
5071: }

5073: #endif /* defined(PETSC_HAVE_P4EST) */