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, §ion);
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(§ion);
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, §ion);
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, §ion);
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, §ion, &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) */