Actual source code: pforest.h
1: #pragma once
3: #include <petscds.h>
4: #include <petsc/private/dmimpl.h>
5: #include <petsc/private/dmforestimpl.h>
6: #include <petsc/private/dmpleximpl.h>
7: #include <petsc/private/dmlabelimpl.h>
8: #include <petsc/private/viewerimpl.h>
9: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
10: #include "petsc_p4est_package.h"
12: #if defined(PETSC_HAVE_P4EST)
14: #if !defined(P4_TO_P8)
15: #include <p4est.h>
16: #include <p4est_extended.h>
17: #include <p4est_geometry.h>
18: #include <p4est_ghost.h>
19: #include <p4est_lnodes.h>
20: #include <p4est_vtk.h>
21: #include <p4est_plex.h>
22: #include <p4est_bits.h>
23: #include <p4est_algorithms.h>
24: #else
25: #include <p8est.h>
26: #include <p8est_extended.h>
27: #include <p8est_geometry.h>
28: #include <p8est_ghost.h>
29: #include <p8est_lnodes.h>
30: #include <p8est_vtk.h>
31: #include <p8est_plex.h>
32: #include <p8est_bits.h>
33: #include <p8est_algorithms.h>
34: #endif
36: typedef enum {
37: PATTERN_HASH,
38: PATTERN_FRACTAL,
39: PATTERN_CORNER,
40: PATTERN_CENTER,
41: PATTERN_COUNT
42: } DMRefinePattern;
43: static const char *DMRefinePatternName[PATTERN_COUNT] = {"hash", "fractal", "corner", "center"};
45: typedef struct _DMRefinePatternCtx {
46: PetscInt corner;
47: PetscBool fractal[P4EST_CHILDREN];
48: PetscReal hashLikelihood;
49: PetscInt maxLevel;
50: p4est_refine_t refine_fn;
51: } DMRefinePatternCtx;
53: static int DMRefinePattern_Corner(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
54: {
55: p4est_quadrant_t root, rootcorner;
56: DMRefinePatternCtx *ctx;
58: ctx = (DMRefinePatternCtx *)p4est->user_pointer;
59: if (quadrant->level >= ctx->maxLevel) return 0;
61: root.x = root.y = 0;
62: #if defined(P4_TO_P8)
63: root.z = 0;
64: #endif
65: root.level = 0;
66: p4est_quadrant_corner_descendant(&root, &rootcorner, ctx->corner, quadrant->level);
67: if (p4est_quadrant_is_equal(quadrant, &rootcorner)) return 1;
68: return 0;
69: }
71: static int DMRefinePattern_Center(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
72: {
73: int cid;
74: p4est_quadrant_t ancestor, ancestorcorner;
75: DMRefinePatternCtx *ctx;
77: ctx = (DMRefinePatternCtx *)p4est->user_pointer;
78: if (quadrant->level >= ctx->maxLevel) return 0;
79: if (quadrant->level <= 1) return 1;
81: p4est_quadrant_ancestor(quadrant, 1, &ancestor);
82: cid = p4est_quadrant_child_id(&ancestor);
83: p4est_quadrant_corner_descendant(&ancestor, &ancestorcorner, P4EST_CHILDREN - 1 - cid, quadrant->level);
84: if (p4est_quadrant_is_equal(quadrant, &ancestorcorner)) return 1;
85: return 0;
86: }
88: static int DMRefinePattern_Fractal(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
89: {
90: int cid;
91: DMRefinePatternCtx *ctx;
93: ctx = (DMRefinePatternCtx *)p4est->user_pointer;
94: if (quadrant->level >= ctx->maxLevel) return 0;
95: if (!quadrant->level) return 1;
96: cid = p4est_quadrant_child_id(quadrant);
97: if (ctx->fractal[cid ^ ((int)(quadrant->level % P4EST_CHILDREN))]) return 1;
98: return 0;
99: }
101: /* simplified from MurmurHash3 by Austin Appleby */
102: #define DMPROT32(x, y) ((x << y) | (x >> (32 - y)))
103: static uint32_t DMPforestHash(const uint32_t *blocks, uint32_t nblocks)
104: {
105: uint32_t c1 = 0xcc9e2d51;
106: uint32_t c2 = 0x1b873593;
107: uint32_t r1 = 15;
108: uint32_t r2 = 13;
109: uint32_t m = 5;
110: uint32_t n = 0xe6546b64;
111: uint32_t hash = 0;
112: int len = nblocks * 4;
113: uint32_t i;
115: for (i = 0; i < nblocks; i++) {
116: uint32_t k;
118: k = blocks[i];
119: k *= c1;
120: k = DMPROT32(k, r1);
121: k *= c2;
123: hash ^= k;
124: hash = DMPROT32(hash, r2) * m + n;
125: }
127: hash ^= len;
128: hash ^= (hash >> 16);
129: hash *= 0x85ebca6b;
130: hash ^= (hash >> 13);
131: hash *= 0xc2b2ae35;
132: hash ^= (hash >> 16);
134: return hash;
135: }
137: #if defined(UINT32_MAX)
138: #define DMP4EST_HASH_MAX UINT32_MAX
139: #else
140: #define DMP4EST_HASH_MAX ((uint32_t)0xffffffff)
141: #endif
143: static int DMRefinePattern_Hash(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
144: {
145: uint32_t data[5];
146: uint32_t result;
147: DMRefinePatternCtx *ctx;
149: ctx = (DMRefinePatternCtx *)p4est->user_pointer;
150: if (quadrant->level >= ctx->maxLevel) return 0;
151: data[0] = ((uint32_t)quadrant->level) << 24;
152: data[1] = (uint32_t)which_tree;
153: data[2] = (uint32_t)quadrant->x;
154: data[3] = (uint32_t)quadrant->y;
155: #if defined(P4_TO_P8)
156: data[4] = (uint32_t)quadrant->z;
157: #endif
159: result = DMPforestHash(data, 2 + P4EST_DIM);
160: if (((double)result / (double)DMP4EST_HASH_MAX) < ctx->hashLikelihood) return 1;
161: return 0;
162: }
164: #define DMConvert_pforest_plex _infix_pforest(DMConvert, _plex)
165: static PetscErrorCode DMConvert_pforest_plex(DM, DMType, DM *);
167: #define DMFTopology_pforest _append_pforest(DMFTopology)
168: typedef struct {
169: PetscInt refct;
170: p4est_connectivity_t *conn;
171: p4est_geometry_t *geom;
172: PetscInt *tree_face_to_uniq; /* p4est does not explicitly enumerate facets, but we must to keep track of labels */
173: } DMFTopology_pforest;
175: #define DM_Forest_pforest _append_pforest(DM_Forest)
176: typedef struct {
177: DMFTopology_pforest *topo;
178: p4est_t *forest;
179: p4est_ghost_t *ghost;
180: p4est_lnodes_t *lnodes;
181: PetscBool partition_for_coarsening;
182: PetscBool coarsen_hierarchy;
183: PetscBool labelsFinalized;
184: PetscBool adaptivitySuccess;
185: PetscInt cLocalStart;
186: PetscInt cLocalEnd;
187: DM plex;
188: char *ghostName;
189: PetscSF pointAdaptToSelfSF;
190: PetscSF pointSelfToAdaptSF;
191: PetscInt *pointAdaptToSelfCids;
192: PetscInt *pointSelfToAdaptCids;
193: } DM_Forest_pforest;
195: #define DM_Forest_geometry_pforest _append_pforest(DM_Forest_geometry)
196: typedef struct {
197: DM base;
198: PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
199: void *mapCtx;
200: PetscInt coordDim;
201: p4est_geometry_t *inner;
202: } DM_Forest_geometry_pforest;
204: #define GeometryMapping_pforest _append_pforest(GeometryMapping)
205: static void GeometryMapping_pforest(p4est_geometry_t *geom, p4est_topidx_t which_tree, const double abc[3], double xyz[3])
206: {
207: DM_Forest_geometry_pforest *geom_pforest = (DM_Forest_geometry_pforest *)geom->user;
208: PetscReal PetscABC[3] = {0.};
209: PetscReal PetscXYZ[3] = {0.};
210: PetscInt i, d = PetscMin(3, geom_pforest->coordDim);
211: double ABC[3];
212: PetscErrorCode ierr;
214: (geom_pforest->inner->X)(geom_pforest->inner, which_tree, abc, ABC);
216: for (i = 0; i < d; i++) PetscABC[i] = ABC[i];
217: ierr = (geom_pforest->map)(geom_pforest->base, (PetscInt)which_tree, geom_pforest->coordDim, PetscABC, PetscXYZ, geom_pforest->mapCtx);
218: PETSC_P4EST_ASSERT(!ierr);
219: for (i = 0; i < d; i++) xyz[i] = PetscXYZ[i];
220: }
222: #define GeometryDestroy_pforest _append_pforest(GeometryDestroy)
223: static void GeometryDestroy_pforest(p4est_geometry_t *geom)
224: {
225: DM_Forest_geometry_pforest *geom_pforest = (DM_Forest_geometry_pforest *)geom->user;
226: PetscErrorCode ierr;
228: p4est_geometry_destroy(geom_pforest->inner);
229: ierr = PetscFree(geom->user);
230: PETSC_P4EST_ASSERT(!ierr);
231: ierr = PetscFree(geom);
232: PETSC_P4EST_ASSERT(!ierr);
233: }
235: #define DMFTopologyDestroy_pforest _append_pforest(DMFTopologyDestroy)
236: static PetscErrorCode DMFTopologyDestroy_pforest(DMFTopology_pforest **topo)
237: {
238: PetscFunctionBegin;
239: if (!(*topo)) PetscFunctionReturn(PETSC_SUCCESS);
240: if (--((*topo)->refct) > 0) {
241: *topo = NULL;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
244: if ((*topo)->geom) PetscCallP4est(p4est_geometry_destroy, ((*topo)->geom));
245: PetscCallP4est(p4est_connectivity_destroy, ((*topo)->conn));
246: PetscCall(PetscFree((*topo)->tree_face_to_uniq));
247: PetscCall(PetscFree(*topo));
248: *topo = NULL;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *, PetscInt **);
254: #define DMFTopologyCreateBrick_pforest _append_pforest(DMFTopologyCreateBrick)
255: static PetscErrorCode DMFTopologyCreateBrick_pforest(DM dm, PetscInt N[], PetscInt P[], PetscReal B[], DMFTopology_pforest **topo, PetscBool useMorton)
256: {
257: double *vertices;
258: PetscInt i, numVerts;
260: PetscFunctionBegin;
261: PetscCheck(useMorton, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Lexicographic ordering not implemented yet");
262: PetscCall(PetscNew(topo));
264: (*topo)->refct = 1;
265: #if !defined(P4_TO_P8)
266: 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));
267: #else
268: 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));
269: #endif
270: numVerts = (*topo)->conn->num_vertices;
271: vertices = (*topo)->conn->vertices;
272: for (i = 0; i < 3 * numVerts; i++) {
273: PetscInt j = i % 3;
275: vertices[i] = B[2 * j] + (vertices[i] / N[j]) * (B[2 * j + 1] - B[2 * j]);
276: }
277: (*topo)->geom = NULL;
278: PetscCall(PforestConnectivityEnumerateFacets((*topo)->conn, &(*topo)->tree_face_to_uniq));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: #define DMFTopologyCreate_pforest _append_pforest(DMFTopologyCreate)
283: static PetscErrorCode DMFTopologyCreate_pforest(DM dm, DMForestTopology topologyName, DMFTopology_pforest **topo)
284: {
285: const char *name = (const char *)topologyName;
286: const char *prefix;
287: PetscBool isBrick, isShell, isSphere, isMoebius;
289: PetscFunctionBegin;
291: PetscAssertPointer(name, 2);
292: PetscAssertPointer(topo, 3);
293: PetscCall(PetscStrcmp(name, "brick", &isBrick));
294: PetscCall(PetscStrcmp(name, "shell", &isShell));
295: PetscCall(PetscStrcmp(name, "sphere", &isSphere));
296: PetscCall(PetscStrcmp(name, "moebius", &isMoebius));
297: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
298: if (isBrick) {
299: PetscBool flgN, flgP, flgM, flgB, useMorton = PETSC_TRUE, periodic = PETSC_FALSE;
300: PetscInt N[3] = {2, 2, 2}, P[3] = {0, 0, 0}, nretN = P4EST_DIM, nretP = P4EST_DIM, nretB = 2 * P4EST_DIM, i;
301: 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};
303: if (dm->setfromoptionscalled) {
304: PetscCall(PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_size", N, &nretN, &flgN));
305: PetscCall(PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_periodicity", P, &nretP, &flgP));
306: PetscCall(PetscOptionsGetRealArray(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_bounds", B, &nretB, &flgB));
307: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_p4est_brick_use_morton_curve", &useMorton, &flgM));
308: PetscCheck(!flgN || nretN == P4EST_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Need to give %d sizes in -dm_p4est_brick_size, gave %" PetscInt_FMT, P4EST_DIM, nretN);
309: PetscCheck(!flgP || nretP == P4EST_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Need to give %d periodicities in -dm_p4est_brick_periodicity, gave %" PetscInt_FMT, P4EST_DIM, nretP);
310: PetscCheck(!flgB || nretB == 2 * P4EST_DIM, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Need to give %d bounds in -dm_p4est_brick_bounds, gave %" PetscInt_FMT, P4EST_DIM, nretP);
311: }
312: for (i = 0; i < P4EST_DIM; i++) {
313: P[i] = (P[i] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE);
314: periodic = (PetscBool)(P[i] || periodic);
315: if (!flgB) B[2 * i + 1] = N[i];
316: if (P[i]) {
317: Lstart[i] = B[2 * i + 0];
318: L[i] = B[2 * i + 1] - B[2 * i + 0];
319: maxCell[i] = 1.1 * (L[i] / N[i]);
320: }
321: }
322: PetscCall(DMFTopologyCreateBrick_pforest(dm, N, P, B, topo, useMorton));
323: if (periodic) PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
324: } else {
325: PetscCall(PetscNew(topo));
327: (*topo)->refct = 1;
328: PetscCallP4estReturn((*topo)->conn, p4est_connectivity_new_byname, (name));
329: (*topo)->geom = NULL;
330: if (isMoebius) PetscCall(DMSetCoordinateDim(dm, 3));
331: #if defined(P4_TO_P8)
332: if (isShell) {
333: PetscReal R2 = 1., R1 = .55;
335: if (dm->setfromoptionscalled) {
336: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_shell_outer_radius", &R2, NULL));
337: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_shell_inner_radius", &R1, NULL));
338: }
339: PetscCallP4estReturn((*topo)->geom, p8est_geometry_new_shell, ((*topo)->conn, R2, R1));
340: } else if (isSphere) {
341: PetscReal R2 = 1., R1 = 0.191728, R0 = 0.039856;
343: if (dm->setfromoptionscalled) {
344: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_outer_radius", &R2, NULL));
345: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_inner_radius", &R1, NULL));
346: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_sphere_core_radius", &R0, NULL));
347: }
348: PetscCallP4estReturn((*topo)->geom, p8est_geometry_new_sphere, ((*topo)->conn, R2, R1, R0));
349: }
350: #endif
351: PetscCall(PforestConnectivityEnumerateFacets((*topo)->conn, &(*topo)->tree_face_to_uniq));
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: #define DMConvert_plex_pforest _append_pforest(DMConvert_plex)
357: static PetscErrorCode DMConvert_plex_pforest(DM dm, DMType newtype, DM *pforest)
358: {
359: MPI_Comm comm;
360: PetscBool isPlex;
361: PetscInt dim;
362: void *ctx;
364: PetscFunctionBegin;
367: comm = PetscObjectComm((PetscObject)dm);
368: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
369: PetscCheck(isPlex, comm, PETSC_ERR_ARG_WRONG, "Expected DM type %s, got %s", DMPLEX, ((PetscObject)dm)->type_name);
370: PetscCall(DMGetDimension(dm, &dim));
371: PetscCheck(dim == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Expected DM dimension %d, got %" PetscInt_FMT, P4EST_DIM, dim);
372: PetscCall(DMCreate(comm, pforest));
373: PetscCall(DMSetType(*pforest, DMPFOREST));
374: PetscCall(DMForestSetBaseDM(*pforest, dm));
375: PetscCall(DMGetApplicationContext(dm, &ctx));
376: PetscCall(DMSetApplicationContext(*pforest, ctx));
377: PetscCall(DMCopyDisc(dm, *pforest));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: #define DMForestDestroy_pforest _append_pforest(DMForestDestroy)
382: static PetscErrorCode DMForestDestroy_pforest(DM dm)
383: {
384: DM_Forest *forest = (DM_Forest *)dm->data;
385: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
387: PetscFunctionBegin;
389: if (pforest->lnodes) PetscCallP4est(p4est_lnodes_destroy, (pforest->lnodes));
390: pforest->lnodes = NULL;
391: if (pforest->ghost) PetscCallP4est(p4est_ghost_destroy, (pforest->ghost));
392: pforest->ghost = NULL;
393: if (pforest->forest) PetscCallP4est(p4est_destroy, (pforest->forest));
394: pforest->forest = NULL;
395: PetscCall(DMFTopologyDestroy_pforest(&pforest->topo));
396: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", NULL));
397: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", NULL));
398: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
399: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", NULL));
400: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL));
401: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", NULL));
402: PetscCall(PetscFree(pforest->ghostName));
403: PetscCall(DMDestroy(&pforest->plex));
404: PetscCall(PetscSFDestroy(&pforest->pointAdaptToSelfSF));
405: PetscCall(PetscSFDestroy(&pforest->pointSelfToAdaptSF));
406: PetscCall(PetscFree(pforest->pointAdaptToSelfCids));
407: PetscCall(PetscFree(pforest->pointSelfToAdaptCids));
408: PetscCall(PetscFree(forest->data));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: #define DMForestTemplate_pforest _append_pforest(DMForestTemplate)
413: static PetscErrorCode DMForestTemplate_pforest(DM dm, DM tdm)
414: {
415: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
416: DM_Forest_pforest *tpforest = (DM_Forest_pforest *)((DM_Forest *)tdm->data)->data;
418: PetscFunctionBegin;
419: if (pforest->topo) pforest->topo->refct++;
420: PetscCall(DMFTopologyDestroy_pforest(&(tpforest->topo)));
421: tpforest->topo = pforest->topo;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: #define DMPlexCreateConnectivity_pforest _append_pforest(DMPlexCreateConnectivity)
426: static PetscErrorCode DMPlexCreateConnectivity_pforest(DM, p4est_connectivity_t **, PetscInt **);
428: typedef struct _PforestAdaptCtx {
429: PetscInt maxLevel;
430: PetscInt minLevel;
431: PetscInt currLevel;
432: PetscBool anyChange;
433: } PforestAdaptCtx;
435: static int pforest_coarsen_currlevel(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
436: {
437: PforestAdaptCtx *ctx = (PforestAdaptCtx *)p4est->user_pointer;
438: PetscInt minLevel = ctx->minLevel;
439: PetscInt currLevel = ctx->currLevel;
441: if (quadrants[0]->level <= minLevel) return 0;
442: return (int)((PetscInt)quadrants[0]->level == currLevel);
443: }
445: static int pforest_coarsen_uniform(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
446: {
447: PforestAdaptCtx *ctx = (PforestAdaptCtx *)p4est->user_pointer;
448: PetscInt minLevel = ctx->minLevel;
450: return (int)((PetscInt)quadrants[0]->level > minLevel);
451: }
453: static int pforest_coarsen_flag_any(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
454: {
455: PetscInt i;
456: PetscBool any = PETSC_FALSE;
457: PforestAdaptCtx *ctx = (PforestAdaptCtx *)p4est->user_pointer;
458: PetscInt minLevel = ctx->minLevel;
460: if (quadrants[0]->level <= minLevel) return 0;
461: for (i = 0; i < P4EST_CHILDREN; i++) {
462: if (quadrants[i]->p.user_int == DM_ADAPT_KEEP) {
463: any = PETSC_FALSE;
464: break;
465: }
466: if (quadrants[i]->p.user_int == DM_ADAPT_COARSEN) {
467: any = PETSC_TRUE;
468: break;
469: }
470: }
471: return any ? 1 : 0;
472: }
474: static int pforest_coarsen_flag_all(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrants[])
475: {
476: PetscInt i;
477: PetscBool all = PETSC_TRUE;
478: PforestAdaptCtx *ctx = (PforestAdaptCtx *)p4est->user_pointer;
479: PetscInt minLevel = ctx->minLevel;
481: if (quadrants[0]->level <= minLevel) return 0;
482: for (i = 0; i < P4EST_CHILDREN; i++) {
483: if (quadrants[i]->p.user_int != DM_ADAPT_COARSEN) {
484: all = PETSC_FALSE;
485: break;
486: }
487: }
488: return all ? 1 : 0;
489: }
491: static void pforest_init_determine(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
492: {
493: quadrant->p.user_int = DM_ADAPT_DETERMINE;
494: }
496: static int pforest_refine_uniform(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: return ((PetscInt)quadrant->level < maxLevel);
502: }
504: static int pforest_refine_flag(p4est_t *p4est, p4est_topidx_t which_tree, p4est_quadrant_t *quadrant)
505: {
506: PforestAdaptCtx *ctx = (PforestAdaptCtx *)p4est->user_pointer;
507: PetscInt maxLevel = ctx->maxLevel;
509: if ((PetscInt)quadrant->level >= maxLevel) return 0;
511: return (quadrant->p.user_int == DM_ADAPT_REFINE);
512: }
514: 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)
515: {
516: PetscMPIInt rank = p4estFrom->mpirank;
517: p4est_topidx_t t;
518: PetscInt toFineLeaves = 0, fromFineLeaves = 0;
520: PetscFunctionBegin;
521: for (t = flt; t <= llt; t++) { /* count roots and leaves */
522: p4est_tree_t *treeFrom = &(((p4est_tree_t *)p4estFrom->trees->array)[t]);
523: p4est_tree_t *treeTo = &(((p4est_tree_t *)p4estTo->trees->array)[t]);
524: p4est_quadrant_t *firstFrom = &treeFrom->first_desc;
525: p4est_quadrant_t *firstTo = &treeTo->first_desc;
526: PetscInt numFrom = (PetscInt)treeFrom->quadrants.elem_count;
527: PetscInt numTo = (PetscInt)treeTo->quadrants.elem_count;
528: p4est_quadrant_t *quadsFrom = (p4est_quadrant_t *)treeFrom->quadrants.array;
529: p4est_quadrant_t *quadsTo = (p4est_quadrant_t *)treeTo->quadrants.array;
530: PetscInt currentFrom, currentTo;
531: PetscInt treeOffsetFrom = (PetscInt)treeFrom->quadrants_offset;
532: PetscInt treeOffsetTo = (PetscInt)treeTo->quadrants_offset;
533: int comp;
535: PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (firstFrom, firstTo));
536: PetscCheck(comp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "non-matching partitions");
538: for (currentFrom = 0, currentTo = 0; currentFrom < numFrom && currentTo < numTo;) {
539: p4est_quadrant_t *quadFrom = &quadsFrom[currentFrom];
540: p4est_quadrant_t *quadTo = &quadsTo[currentTo];
542: if (quadFrom->level == quadTo->level) {
543: if (toLeaves) {
544: toLeaves[toFineLeaves] = currentTo + treeOffsetTo + ToOffset;
545: fromRoots[toFineLeaves].rank = rank;
546: fromRoots[toFineLeaves].index = currentFrom + treeOffsetFrom + FromOffset;
547: }
548: toFineLeaves++;
549: currentFrom++;
550: currentTo++;
551: } else {
552: int fromIsAncestor;
554: PetscCallP4estReturn(fromIsAncestor, p4est_quadrant_is_ancestor, (quadFrom, quadTo));
555: if (fromIsAncestor) {
556: p4est_quadrant_t lastDesc;
558: if (toLeaves) {
559: toLeaves[toFineLeaves] = currentTo + treeOffsetTo + ToOffset;
560: fromRoots[toFineLeaves].rank = rank;
561: fromRoots[toFineLeaves].index = currentFrom + treeOffsetFrom + FromOffset;
562: }
563: toFineLeaves++;
564: currentTo++;
565: PetscCallP4est(p4est_quadrant_last_descendant, (quadFrom, &lastDesc, quadTo->level));
566: PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (quadTo, &lastDesc));
567: if (comp) currentFrom++;
568: } else {
569: p4est_quadrant_t lastDesc;
571: if (fromLeaves) {
572: fromLeaves[fromFineLeaves] = currentFrom + treeOffsetFrom + FromOffset;
573: toRoots[fromFineLeaves].rank = rank;
574: toRoots[fromFineLeaves].index = currentTo + treeOffsetTo + ToOffset;
575: }
576: fromFineLeaves++;
577: currentFrom++;
578: PetscCallP4est(p4est_quadrant_last_descendant, (quadTo, &lastDesc, quadFrom->level));
579: PetscCallP4estReturn(comp, p4est_quadrant_is_equal, (quadFrom, &lastDesc));
580: if (comp) currentTo++;
581: }
582: }
583: }
584: }
585: *toFineLeavesCount = toFineLeaves;
586: *fromFineLeavesCount = fromFineLeaves;
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /* Compute the maximum level across all the trees */
591: static PetscErrorCode DMPforestGetRefinementLevel(DM dm, PetscInt *lev)
592: {
593: p4est_topidx_t t, flt, llt;
594: DM_Forest *forest = (DM_Forest *)dm->data;
595: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
596: PetscInt maxlevelloc = 0;
597: p4est_t *p4est;
599: PetscFunctionBegin;
600: PetscCheck(pforest, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing DM_Forest_pforest");
601: PetscCheck(pforest->forest, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing p4est_t");
602: p4est = pforest->forest;
603: flt = p4est->first_local_tree;
604: llt = p4est->last_local_tree;
605: for (t = flt; t <= llt; t++) {
606: p4est_tree_t *tree = &(((p4est_tree_t *)p4est->trees->array)[t]);
607: maxlevelloc = PetscMax((PetscInt)tree->maxlevel, maxlevelloc);
608: }
609: PetscCall(MPIU_Allreduce(&maxlevelloc, lev, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: /* Puts identity in coarseToFine */
614: /* assumes a matching partition */
615: static PetscErrorCode DMPforestComputeLocalCellTransferSF(MPI_Comm comm, p4est_t *p4estFrom, PetscInt FromOffset, p4est_t *p4estTo, PetscInt ToOffset, PetscSF *fromCoarseToFine, PetscSF *toCoarseFromFine)
616: {
617: p4est_topidx_t flt, llt;
618: PetscSF fromCoarse, toCoarse;
619: PetscInt numRootsFrom, numRootsTo, numLeavesFrom, numLeavesTo;
620: PetscInt *fromLeaves = NULL, *toLeaves = NULL;
621: PetscSFNode *fromRoots = NULL, *toRoots = NULL;
623: PetscFunctionBegin;
624: flt = p4estFrom->first_local_tree;
625: llt = p4estFrom->last_local_tree;
626: PetscCall(PetscSFCreate(comm, &fromCoarse));
627: if (toCoarseFromFine) PetscCall(PetscSFCreate(comm, &toCoarse));
628: numRootsFrom = p4estFrom->local_num_quadrants + FromOffset;
629: numRootsTo = p4estTo->local_num_quadrants + ToOffset;
630: PetscCall(DMPforestComputeLocalCellTransferSF_loop(p4estFrom, FromOffset, p4estTo, ToOffset, flt, llt, &numLeavesTo, NULL, NULL, &numLeavesFrom, NULL, NULL));
631: PetscCall(PetscMalloc1(numLeavesTo, &toLeaves));
632: PetscCall(PetscMalloc1(numLeavesTo, &fromRoots));
633: if (toCoarseFromFine) {
634: PetscCall(PetscMalloc1(numLeavesFrom, &fromLeaves));
635: PetscCall(PetscMalloc1(numLeavesFrom, &fromRoots));
636: }
637: PetscCall(DMPforestComputeLocalCellTransferSF_loop(p4estFrom, FromOffset, p4estTo, ToOffset, flt, llt, &numLeavesTo, toLeaves, fromRoots, &numLeavesFrom, fromLeaves, toRoots));
638: if (!ToOffset && (numLeavesTo == numRootsTo)) { /* compress */
639: PetscCall(PetscFree(toLeaves));
640: PetscCall(PetscSFSetGraph(fromCoarse, numRootsFrom, numLeavesTo, NULL, PETSC_OWN_POINTER, fromRoots, PETSC_OWN_POINTER));
641: } else PetscCall(PetscSFSetGraph(fromCoarse, numRootsFrom, numLeavesTo, toLeaves, PETSC_OWN_POINTER, fromRoots, PETSC_OWN_POINTER));
642: *fromCoarseToFine = fromCoarse;
643: if (toCoarseFromFine) {
644: PetscCall(PetscSFSetGraph(toCoarse, numRootsTo, numLeavesFrom, fromLeaves, PETSC_OWN_POINTER, toRoots, PETSC_OWN_POINTER));
645: *toCoarseFromFine = toCoarse;
646: }
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: /* range of processes whose B sections overlap this ranks A section */
651: static PetscErrorCode DMPforestComputeOverlappingRanks(PetscMPIInt size, PetscMPIInt rank, p4est_t *p4estA, p4est_t *p4estB, PetscInt *startB, PetscInt *endB)
652: {
653: p4est_quadrant_t *myCoarseStart = &(p4estA->global_first_position[rank]);
654: p4est_quadrant_t *myCoarseEnd = &(p4estA->global_first_position[rank + 1]);
655: p4est_quadrant_t *globalFirstB = p4estB->global_first_position;
657: PetscFunctionBegin;
658: *startB = -1;
659: *endB = -1;
660: if (p4estA->local_num_quadrants) {
661: PetscInt lo, hi, guess;
662: /* binary search to find interval containing myCoarseStart */
663: lo = 0;
664: hi = size;
665: guess = rank;
666: while (1) {
667: int startCompMy, myCompEnd;
669: PetscCallP4estReturn(startCompMy, p4est_quadrant_compare_piggy, (&globalFirstB[guess], myCoarseStart));
670: PetscCallP4estReturn(myCompEnd, p4est_quadrant_compare_piggy, (myCoarseStart, &globalFirstB[guess + 1]));
671: if (startCompMy <= 0 && myCompEnd < 0) {
672: *startB = guess;
673: break;
674: } else if (startCompMy > 0) { /* guess is to high */
675: hi = guess;
676: } else { /* guess is to low */
677: lo = guess + 1;
678: }
679: guess = lo + (hi - lo) / 2;
680: }
681: /* reset bounds, but not guess */
682: lo = 0;
683: hi = size;
684: while (1) {
685: int startCompMy, myCompEnd;
687: PetscCallP4estReturn(startCompMy, p4est_quadrant_compare_piggy, (&globalFirstB[guess], myCoarseEnd));
688: PetscCallP4estReturn(myCompEnd, p4est_quadrant_compare_piggy, (myCoarseEnd, &globalFirstB[guess + 1]));
689: if (startCompMy < 0 && myCompEnd <= 0) { /* notice that the comparison operators are different from above */
690: *endB = guess + 1;
691: break;
692: } else if (startCompMy >= 0) { /* guess is to high */
693: hi = guess;
694: } else { /* guess is to low */
695: lo = guess + 1;
696: }
697: guess = lo + (hi - lo) / 2;
698: }
699: }
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode DMPforestGetPlex(DM, DM *);
705: #define DMSetUp_pforest _append_pforest(DMSetUp)
706: static PetscErrorCode DMSetUp_pforest(DM dm)
707: {
708: DM_Forest *forest = (DM_Forest *)dm->data;
709: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
710: DM base, adaptFrom;
711: DMForestTopology topoName;
712: PetscSF preCoarseToFine = NULL, coarseToPreFine = NULL;
713: PforestAdaptCtx ctx;
715: PetscFunctionBegin;
716: ctx.minLevel = PETSC_MAX_INT;
717: ctx.maxLevel = 0;
718: ctx.currLevel = 0;
719: ctx.anyChange = PETSC_FALSE;
720: /* sanity check */
721: PetscCall(DMForestGetAdaptivityForest(dm, &adaptFrom));
722: PetscCall(DMForestGetBaseDM(dm, &base));
723: PetscCall(DMForestGetTopology(dm, &topoName));
724: PetscCheck(adaptFrom || base || topoName, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "A forest needs a topology, a base DM, or a DM to adapt from");
726: /* === Step 1: DMFTopology === */
727: if (adaptFrom) { /* reference already created topology */
728: PetscBool ispforest;
729: DM_Forest *aforest = (DM_Forest *)adaptFrom->data;
730: DM_Forest_pforest *apforest = (DM_Forest_pforest *)aforest->data;
732: PetscCall(PetscObjectTypeCompare((PetscObject)adaptFrom, DMPFOREST, &ispforest));
733: PetscCheck(ispforest, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NOTSAMETYPE, "Trying to adapt from %s, which is not %s", ((PetscObject)adaptFrom)->type_name, DMPFOREST);
734: PetscCheck(apforest->topo, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The pre-adaptation forest must have a topology");
735: PetscCall(DMSetUp(adaptFrom));
736: PetscCall(DMForestGetBaseDM(dm, &base));
737: PetscCall(DMForestGetTopology(dm, &topoName));
738: } else if (base) { /* construct a connectivity from base */
739: PetscBool isPlex, isDA;
741: PetscCall(PetscObjectGetName((PetscObject)base, &topoName));
742: PetscCall(DMForestSetTopology(dm, topoName));
743: PetscCall(PetscObjectTypeCompare((PetscObject)base, DMPLEX, &isPlex));
744: PetscCall(PetscObjectTypeCompare((PetscObject)base, DMDA, &isDA));
745: if (isPlex) {
746: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
747: PetscInt depth;
748: PetscMPIInt size;
749: p4est_connectivity_t *conn = NULL;
750: DMFTopology_pforest *topo;
751: PetscInt *tree_face_to_uniq = NULL;
753: PetscCall(DMPlexGetDepth(base, &depth));
754: if (depth == 1) {
755: DM connDM;
757: PetscCall(DMPlexInterpolate(base, &connDM));
758: base = connDM;
759: PetscCall(DMForestSetBaseDM(dm, base));
760: PetscCall(DMDestroy(&connDM));
761: } else PetscCheck(depth == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Base plex is neither interpolated nor uninterpolated? depth %" PetscInt_FMT ", expected 2 or %d", depth, P4EST_DIM + 1);
762: PetscCallMPI(MPI_Comm_size(comm, &size));
763: if (size > 1) {
764: DM dmRedundant;
765: PetscSF sf;
767: PetscCall(DMPlexGetRedundantDM(base, &sf, &dmRedundant));
768: PetscCheck(dmRedundant, comm, PETSC_ERR_PLIB, "Could not create redundant DM");
769: PetscCall(PetscObjectCompose((PetscObject)dmRedundant, "_base_migration_sf", (PetscObject)sf));
770: PetscCall(PetscSFDestroy(&sf));
771: base = dmRedundant;
772: PetscCall(DMForestSetBaseDM(dm, base));
773: PetscCall(DMDestroy(&dmRedundant));
774: }
775: PetscCall(DMViewFromOptions(base, NULL, "-dm_p4est_base_view"));
776: PetscCall(DMPlexCreateConnectivity_pforest(base, &conn, &tree_face_to_uniq));
777: PetscCall(PetscNew(&topo));
778: topo->refct = 1;
779: topo->conn = conn;
780: topo->geom = NULL;
781: {
782: PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
783: void *mapCtx;
785: PetscCall(DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx));
786: if (map) {
787: DM_Forest_geometry_pforest *geom_pforest;
788: p4est_geometry_t *geom;
790: PetscCall(PetscNew(&geom_pforest));
791: PetscCall(DMGetCoordinateDim(dm, &geom_pforest->coordDim));
792: geom_pforest->map = map;
793: geom_pforest->mapCtx = mapCtx;
794: PetscCallP4estReturn(geom_pforest->inner, p4est_geometry_new_connectivity, (conn));
795: PetscCall(PetscNew(&geom));
796: geom->name = topoName;
797: geom->user = geom_pforest;
798: geom->X = GeometryMapping_pforest;
799: geom->destroy = GeometryDestroy_pforest;
800: topo->geom = geom;
801: }
802: }
803: topo->tree_face_to_uniq = tree_face_to_uniq;
804: pforest->topo = topo;
805: } else PetscCheck(!isDA, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Not implemented yet");
806: #if 0
807: PetscInt N[3], P[3];
809: /* get the sizes, periodicities */
810: /* ... */
811: /* don't use Morton order */
812: PetscCall(DMFTopologyCreateBrick_pforest(dm,N,P,&pforest->topo,PETSC_FALSE));
813: #endif
814: {
815: PetscInt numLabels, l;
817: PetscCall(DMGetNumLabels(base, &numLabels));
818: for (l = 0; l < numLabels; l++) {
819: PetscBool isDepth, isGhost, isVTK, isDim, isCellType;
820: DMLabel label, labelNew;
821: PetscInt defVal;
822: const char *name;
824: PetscCall(DMGetLabelName(base, l, &name));
825: PetscCall(DMGetLabelByNum(base, l, &label));
826: PetscCall(PetscStrcmp(name, "depth", &isDepth));
827: if (isDepth) continue;
828: PetscCall(PetscStrcmp(name, "dim", &isDim));
829: if (isDim) continue;
830: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
831: if (isCellType) continue;
832: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
833: if (isGhost) continue;
834: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
835: if (isVTK) continue;
836: PetscCall(DMCreateLabel(dm, name));
837: PetscCall(DMGetLabel(dm, name, &labelNew));
838: PetscCall(DMLabelGetDefaultValue(label, &defVal));
839: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
840: }
841: /* map dm points (internal plex) to base
842: we currently create the subpoint_map for the entire hierarchy, starting from the finest forest
843: and propagating back to the coarsest
844: This is not an optimal approach, since we need the map only on the coarsest level
845: during DMForestTransferVecFromBase */
846: PetscCall(DMForestGetMinimumRefinement(dm, &l));
847: if (!l) PetscCall(DMCreateLabel(dm, "_forest_base_subpoint_map"));
848: }
849: } else { /* construct from topology name */
850: DMFTopology_pforest *topo;
852: PetscCall(DMFTopologyCreate_pforest(dm, topoName, &topo));
853: pforest->topo = topo;
854: /* TODO: construct base? */
855: }
857: /* === Step 2: get the leaves of the forest === */
858: if (adaptFrom) { /* start with the old forest */
859: DMLabel adaptLabel;
860: PetscInt defaultValue;
861: PetscInt numValues, numValuesGlobal, cLocalStart, count;
862: DM_Forest *aforest = (DM_Forest *)adaptFrom->data;
863: DM_Forest_pforest *apforest = (DM_Forest_pforest *)aforest->data;
864: PetscBool computeAdaptSF;
865: p4est_topidx_t flt, llt, t;
867: flt = apforest->forest->first_local_tree;
868: llt = apforest->forest->last_local_tree;
869: cLocalStart = apforest->cLocalStart;
870: PetscCall(DMForestGetComputeAdaptivitySF(dm, &computeAdaptSF));
871: PetscCallP4estReturn(pforest->forest, p4est_copy, (apforest->forest, 0)); /* 0 indicates no data copying */
872: PetscCall(DMForestGetAdaptivityLabel(dm, &adaptLabel));
873: if (adaptLabel) {
874: /* apply the refinement/coarsening by flags, plus minimum/maximum refinement */
875: PetscCall(DMLabelGetNumValues(adaptLabel, &numValues));
876: PetscCall(MPIU_Allreduce(&numValues, &numValuesGlobal, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)adaptFrom)));
877: PetscCall(DMLabelGetDefaultValue(adaptLabel, &defaultValue));
878: if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN_LAST) { /* uniform coarsen of the last level only (equivalent to DM_ADAPT_COARSEN for conforming grids) */
879: PetscCall(DMForestGetMinimumRefinement(dm, &ctx.minLevel));
880: PetscCall(DMPforestGetRefinementLevel(dm, &ctx.currLevel));
881: pforest->forest->user_pointer = (void *)&ctx;
882: PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_currlevel, NULL));
883: pforest->forest->user_pointer = (void *)dm;
884: PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
885: /* we will have to change the offset after we compute the overlap */
886: if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL));
887: } else if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN) { /* uniform coarsen */
888: PetscCall(DMForestGetMinimumRefinement(dm, &ctx.minLevel));
889: pforest->forest->user_pointer = (void *)&ctx;
890: PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_uniform, NULL));
891: pforest->forest->user_pointer = (void *)dm;
892: PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
893: /* we will have to change the offset after we compute the overlap */
894: if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL));
895: } else if (!numValuesGlobal && defaultValue == DM_ADAPT_REFINE) { /* uniform refine */
896: PetscCall(DMForestGetMaximumRefinement(dm, &ctx.maxLevel));
897: pforest->forest->user_pointer = (void *)&ctx;
898: PetscCallP4est(p4est_refine, (pforest->forest, 0, pforest_refine_uniform, NULL));
899: pforest->forest->user_pointer = (void *)dm;
900: PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
901: /* we will have to change the offset after we compute the overlap */
902: if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), apforest->forest, apforest->cLocalStart, pforest->forest, 0, &preCoarseToFine, NULL));
903: } else if (numValuesGlobal) {
904: p4est_t *p4est = pforest->forest;
905: PetscInt *cellFlags;
906: DMForestAdaptivityStrategy strategy;
907: PetscSF cellSF;
908: PetscInt c, cStart, cEnd;
909: PetscBool adaptAny;
911: PetscCall(DMForestGetMaximumRefinement(dm, &ctx.maxLevel));
912: PetscCall(DMForestGetMinimumRefinement(dm, &ctx.minLevel));
913: PetscCall(DMForestGetAdaptivityStrategy(dm, &strategy));
914: PetscCall(PetscStrncmp(strategy, "any", 3, &adaptAny));
915: PetscCall(DMForestGetCellChart(adaptFrom, &cStart, &cEnd));
916: PetscCall(DMForestGetCellSF(adaptFrom, &cellSF));
917: PetscCall(PetscMalloc1(cEnd - cStart, &cellFlags));
918: for (c = cStart; c < cEnd; c++) PetscCall(DMLabelGetValue(adaptLabel, c, &cellFlags[c - cStart]));
919: if (cellSF) {
920: if (adaptAny) {
921: PetscCall(PetscSFReduceBegin(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MAX));
922: PetscCall(PetscSFReduceEnd(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MAX));
923: } else {
924: PetscCall(PetscSFReduceBegin(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MIN));
925: PetscCall(PetscSFReduceEnd(cellSF, MPIU_INT, cellFlags, cellFlags, MPI_MIN));
926: }
927: }
928: for (t = flt, count = cLocalStart; t <= llt; t++) {
929: p4est_tree_t *tree = &(((p4est_tree_t *)p4est->trees->array)[t]);
930: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count, i;
931: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
933: for (i = 0; i < numQuads; i++) {
934: p4est_quadrant_t *q = &quads[i];
935: q->p.user_int = cellFlags[count++];
936: }
937: }
938: PetscCall(PetscFree(cellFlags));
940: pforest->forest->user_pointer = (void *)&ctx;
941: if (adaptAny) PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_flag_any, pforest_init_determine));
942: else PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_flag_all, pforest_init_determine));
943: PetscCallP4est(p4est_refine, (pforest->forest, 0, pforest_refine_flag, NULL));
944: pforest->forest->user_pointer = (void *)dm;
945: PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
946: if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), apforest->forest, apforest->cLocalStart, pforest->forest, 0, &preCoarseToFine, &coarseToPreFine));
947: }
948: for (t = flt, count = cLocalStart; t <= llt; t++) {
949: p4est_tree_t *atree = &(((p4est_tree_t *)apforest->forest->trees->array)[t]);
950: p4est_tree_t *tree = &(((p4est_tree_t *)pforest->forest->trees->array)[t]);
951: PetscInt anumQuads = (PetscInt)atree->quadrants.elem_count, i;
952: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
953: p4est_quadrant_t *aquads = (p4est_quadrant_t *)atree->quadrants.array;
954: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
956: if (anumQuads != numQuads) {
957: ctx.anyChange = PETSC_TRUE;
958: } else {
959: for (i = 0; i < numQuads; i++) {
960: p4est_quadrant_t *aq = &aquads[i];
961: p4est_quadrant_t *q = &quads[i];
963: if (aq->level != q->level) {
964: ctx.anyChange = PETSC_TRUE;
965: break;
966: }
967: }
968: }
969: if (ctx.anyChange) break;
970: }
971: }
972: {
973: PetscInt numLabels, l;
975: PetscCall(DMGetNumLabels(adaptFrom, &numLabels));
976: for (l = 0; l < numLabels; l++) {
977: PetscBool isDepth, isCellType, isGhost, isVTK;
978: DMLabel label, labelNew;
979: PetscInt defVal;
980: const char *name;
982: PetscCall(DMGetLabelName(adaptFrom, l, &name));
983: PetscCall(DMGetLabelByNum(adaptFrom, l, &label));
984: PetscCall(PetscStrcmp(name, "depth", &isDepth));
985: if (isDepth) continue;
986: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
987: if (isCellType) continue;
988: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
989: if (isGhost) continue;
990: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
991: if (isVTK) continue;
992: PetscCall(DMCreateLabel(dm, name));
993: PetscCall(DMGetLabel(dm, name, &labelNew));
994: PetscCall(DMLabelGetDefaultValue(label, &defVal));
995: PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
996: }
997: }
998: } else { /* initial */
999: PetscInt initLevel, minLevel;
1000: #if defined(PETSC_HAVE_MPIUNI)
1001: sc_MPI_Comm comm = sc_MPI_COMM_WORLD;
1002: #else
1003: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1004: #endif
1006: PetscCall(DMForestGetInitialRefinement(dm, &initLevel));
1007: PetscCall(DMForestGetMinimumRefinement(dm, &minLevel));
1008: PetscCallP4estReturn(pforest->forest, p4est_new_ext,
1009: (comm, pforest->topo->conn, 0, /* minimum number of quadrants per processor */
1010: initLevel, /* level of refinement */
1011: 1, /* uniform refinement */
1012: 0, /* we don't allocate any per quadrant data */
1013: NULL, /* there is no special quadrant initialization */
1014: (void *)dm)); /* this dm is the user context */
1016: if (initLevel > minLevel) pforest->coarsen_hierarchy = PETSC_TRUE;
1017: if (dm->setfromoptionscalled) {
1018: PetscBool flgPattern, flgFractal;
1019: PetscInt corner = 0;
1020: PetscInt corners[P4EST_CHILDREN], ncorner = P4EST_CHILDREN;
1021: PetscReal likelihood = 1. / P4EST_DIM;
1022: PetscInt pattern;
1023: const char *prefix;
1025: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1026: PetscCall(PetscOptionsGetEList(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_pattern", DMRefinePatternName, PATTERN_COUNT, &pattern, &flgPattern));
1027: PetscCall(PetscOptionsGetInt(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_corner", &corner, NULL));
1028: PetscCall(PetscOptionsGetIntArray(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_fractal_corners", corners, &ncorner, &flgFractal));
1029: PetscCall(PetscOptionsGetReal(((PetscObject)dm)->options, prefix, "-dm_p4est_refine_hash_likelihood", &likelihood, NULL));
1031: if (flgPattern) {
1032: DMRefinePatternCtx *ctx;
1033: PetscInt maxLevel;
1035: PetscCall(DMForestGetMaximumRefinement(dm, &maxLevel));
1036: PetscCall(PetscNew(&ctx));
1037: ctx->maxLevel = PetscMin(maxLevel, P4EST_QMAXLEVEL);
1038: if (initLevel + ctx->maxLevel > minLevel) pforest->coarsen_hierarchy = PETSC_TRUE;
1039: switch (pattern) {
1040: case PATTERN_HASH:
1041: ctx->refine_fn = DMRefinePattern_Hash;
1042: ctx->hashLikelihood = likelihood;
1043: break;
1044: case PATTERN_CORNER:
1045: ctx->corner = corner;
1046: ctx->refine_fn = DMRefinePattern_Corner;
1047: break;
1048: case PATTERN_CENTER:
1049: ctx->refine_fn = DMRefinePattern_Center;
1050: break;
1051: case PATTERN_FRACTAL:
1052: if (flgFractal) {
1053: PetscInt i;
1055: for (i = 0; i < ncorner; i++) ctx->fractal[corners[i]] = PETSC_TRUE;
1056: } else {
1057: #if !defined(P4_TO_P8)
1058: ctx->fractal[0] = ctx->fractal[1] = ctx->fractal[2] = PETSC_TRUE;
1059: #else
1060: ctx->fractal[0] = ctx->fractal[3] = ctx->fractal[5] = ctx->fractal[6] = PETSC_TRUE;
1061: #endif
1062: }
1063: ctx->refine_fn = DMRefinePattern_Fractal;
1064: break;
1065: default:
1066: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Not a valid refinement pattern");
1067: }
1069: pforest->forest->user_pointer = (void *)ctx;
1070: PetscCallP4est(p4est_refine, (pforest->forest, 1, ctx->refine_fn, NULL));
1071: PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL));
1072: PetscCall(PetscFree(ctx));
1073: pforest->forest->user_pointer = (void *)dm;
1074: }
1075: }
1076: }
1077: if (pforest->coarsen_hierarchy) {
1078: PetscInt initLevel, currLevel, minLevel;
1080: PetscCall(DMPforestGetRefinementLevel(dm, &currLevel));
1081: PetscCall(DMForestGetInitialRefinement(dm, &initLevel));
1082: PetscCall(DMForestGetMinimumRefinement(dm, &minLevel));
1083: if (currLevel > minLevel) {
1084: DM_Forest_pforest *coarse_pforest;
1085: DMLabel coarsen;
1086: DM coarseDM;
1088: PetscCall(DMForestTemplate(dm, MPI_COMM_NULL, &coarseDM));
1089: PetscCall(DMForestSetAdaptivityPurpose(coarseDM, DM_ADAPT_COARSEN));
1090: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "coarsen", &coarsen));
1091: PetscCall(DMLabelSetDefaultValue(coarsen, DM_ADAPT_COARSEN));
1092: PetscCall(DMForestSetAdaptivityLabel(coarseDM, coarsen));
1093: PetscCall(DMLabelDestroy(&coarsen));
1094: PetscCall(DMSetCoarseDM(dm, coarseDM));
1095: PetscCall(PetscObjectDereference((PetscObject)coarseDM));
1096: initLevel = currLevel == initLevel ? initLevel - 1 : initLevel;
1097: PetscCall(DMForestSetInitialRefinement(coarseDM, initLevel));
1098: PetscCall(DMForestSetMinimumRefinement(coarseDM, minLevel));
1099: coarse_pforest = (DM_Forest_pforest *)((DM_Forest *)coarseDM->data)->data;
1100: coarse_pforest->coarsen_hierarchy = PETSC_TRUE;
1101: }
1102: }
1104: { /* repartitioning and overlap */
1105: PetscMPIInt size, rank;
1107: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
1108: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1109: if ((size > 1) && (pforest->partition_for_coarsening || forest->cellWeights || forest->weightCapacity != 1. || forest->weightsFactor != 1.)) {
1110: PetscBool copyForest = PETSC_FALSE;
1111: p4est_t *forest_copy = NULL;
1112: p4est_gloidx_t shipped = 0;
1114: if (preCoarseToFine || coarseToPreFine) copyForest = PETSC_TRUE;
1115: if (copyForest) PetscCallP4estReturn(forest_copy, p4est_copy, (pforest->forest, 0));
1117: if (!forest->cellWeights && forest->weightCapacity == 1. && forest->weightsFactor == 1.) {
1118: PetscCallP4estReturn(shipped, p4est_partition_ext, (pforest->forest, (int)pforest->partition_for_coarsening, NULL));
1119: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Non-uniform partition cases not implemented yet");
1120: if (shipped) ctx.anyChange = PETSC_TRUE;
1121: if (forest_copy) {
1122: if (preCoarseToFine || coarseToPreFine) {
1123: PetscSF repartSF; /* repartSF has roots in the old partition */
1124: PetscInt pStart = -1, pEnd = -1, p;
1125: PetscInt numRoots, numLeaves;
1126: PetscSFNode *repartRoots;
1127: p4est_gloidx_t postStart = pforest->forest->global_first_quadrant[rank];
1128: p4est_gloidx_t postEnd = pforest->forest->global_first_quadrant[rank + 1];
1129: p4est_gloidx_t partOffset = postStart;
1131: numRoots = (PetscInt)(forest_copy->global_first_quadrant[rank + 1] - forest_copy->global_first_quadrant[rank]);
1132: numLeaves = (PetscInt)(postEnd - postStart);
1133: PetscCall(DMPforestComputeOverlappingRanks(size, rank, pforest->forest, forest_copy, &pStart, &pEnd));
1134: PetscCall(PetscMalloc1((PetscInt)pforest->forest->local_num_quadrants, &repartRoots));
1135: for (p = pStart; p < pEnd; p++) {
1136: p4est_gloidx_t preStart = forest_copy->global_first_quadrant[p];
1137: p4est_gloidx_t preEnd = forest_copy->global_first_quadrant[p + 1];
1138: PetscInt q;
1140: if (preEnd == preStart) continue;
1141: PetscCheck(preStart <= postStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad partition overlap computation");
1142: preEnd = preEnd > postEnd ? postEnd : preEnd;
1143: for (q = partOffset; q < preEnd; q++) {
1144: repartRoots[q - postStart].rank = p;
1145: repartRoots[q - postStart].index = partOffset - preStart;
1146: }
1147: partOffset = preEnd;
1148: }
1149: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &repartSF));
1150: PetscCall(PetscSFSetGraph(repartSF, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, repartRoots, PETSC_OWN_POINTER));
1151: PetscCall(PetscSFSetUp(repartSF));
1152: if (preCoarseToFine) {
1153: PetscSF repartSFembed, preCoarseToFineNew;
1154: PetscInt nleaves;
1155: const PetscInt *leaves;
1157: PetscCall(PetscSFSetUp(preCoarseToFine));
1158: PetscCall(PetscSFGetGraph(preCoarseToFine, NULL, &nleaves, &leaves, NULL));
1159: if (leaves) {
1160: PetscCall(PetscSFCreateEmbeddedRootSF(repartSF, nleaves, leaves, &repartSFembed));
1161: } else {
1162: repartSFembed = repartSF;
1163: PetscCall(PetscObjectReference((PetscObject)repartSFembed));
1164: }
1165: PetscCall(PetscSFCompose(preCoarseToFine, repartSFembed, &preCoarseToFineNew));
1166: PetscCall(PetscSFDestroy(&preCoarseToFine));
1167: PetscCall(PetscSFDestroy(&repartSFembed));
1168: preCoarseToFine = preCoarseToFineNew;
1169: }
1170: if (coarseToPreFine) {
1171: PetscSF repartSFinv, coarseToPreFineNew;
1173: PetscCall(PetscSFCreateInverseSF(repartSF, &repartSFinv));
1174: PetscCall(PetscSFCompose(repartSFinv, coarseToPreFine, &coarseToPreFineNew));
1175: PetscCall(PetscSFDestroy(&coarseToPreFine));
1176: PetscCall(PetscSFDestroy(&repartSFinv));
1177: coarseToPreFine = coarseToPreFineNew;
1178: }
1179: PetscCall(PetscSFDestroy(&repartSF));
1180: }
1181: PetscCallP4est(p4est_destroy, (forest_copy));
1182: }
1183: }
1184: if (size > 1) {
1185: PetscInt overlap;
1187: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
1189: if (adaptFrom) {
1190: PetscInt aoverlap;
1192: PetscCall(DMForestGetPartitionOverlap(adaptFrom, &aoverlap));
1193: if (aoverlap != overlap) ctx.anyChange = PETSC_TRUE;
1194: }
1196: if (overlap > 0) {
1197: PetscInt i, cLocalStart;
1198: PetscInt cEnd;
1199: PetscSF preCellSF = NULL, cellSF = NULL;
1201: PetscCallP4estReturn(pforest->ghost, p4est_ghost_new, (pforest->forest, P4EST_CONNECT_FULL));
1202: PetscCallP4estReturn(pforest->lnodes, p4est_lnodes_new, (pforest->forest, pforest->ghost, -P4EST_DIM));
1203: PetscCallP4est(p4est_ghost_support_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));
1204: for (i = 1; i < overlap; i++) PetscCallP4est(p4est_ghost_expand_by_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));
1206: cLocalStart = pforest->cLocalStart = pforest->ghost->proc_offsets[rank];
1207: cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[size];
1209: /* shift sfs by cLocalStart, expand by cell SFs */
1210: if (preCoarseToFine || coarseToPreFine) {
1211: if (adaptFrom) PetscCall(DMForestGetCellSF(adaptFrom, &preCellSF));
1212: dm->setupcalled = PETSC_TRUE;
1213: PetscCall(DMForestGetCellSF(dm, &cellSF));
1214: }
1215: if (preCoarseToFine) {
1216: PetscSF preCoarseToFineNew;
1217: PetscInt nleaves, nroots, *leavesNew, i, nleavesNew;
1218: const PetscInt *leaves;
1219: const PetscSFNode *remotes;
1220: PetscSFNode *remotesAll;
1222: PetscCall(PetscSFSetUp(preCoarseToFine));
1223: PetscCall(PetscSFGetGraph(preCoarseToFine, &nroots, &nleaves, &leaves, &remotes));
1224: PetscCall(PetscMalloc1(cEnd, &remotesAll));
1225: for (i = 0; i < cEnd; i++) {
1226: remotesAll[i].rank = -1;
1227: remotesAll[i].index = -1;
1228: }
1229: for (i = 0; i < nleaves; i++) remotesAll[(leaves ? leaves[i] : i) + cLocalStart] = remotes[i];
1230: PetscCall(PetscSFSetUp(cellSF));
1231: PetscCall(PetscSFBcastBegin(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE));
1232: PetscCall(PetscSFBcastEnd(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE));
1233: nleavesNew = 0;
1234: for (i = 0; i < nleaves; i++) {
1235: if (remotesAll[i].rank >= 0) nleavesNew++;
1236: }
1237: PetscCall(PetscMalloc1(nleavesNew, &leavesNew));
1238: nleavesNew = 0;
1239: for (i = 0; i < nleaves; i++) {
1240: if (remotesAll[i].rank >= 0) {
1241: leavesNew[nleavesNew] = i;
1242: if (i > nleavesNew) remotesAll[nleavesNew] = remotesAll[i];
1243: nleavesNew++;
1244: }
1245: }
1246: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &preCoarseToFineNew));
1247: if (nleavesNew < cEnd) {
1248: PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, leavesNew, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES));
1249: } else { /* all cells are leaves */
1250: PetscCall(PetscFree(leavesNew));
1251: PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, NULL, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES));
1252: }
1253: PetscCall(PetscFree(remotesAll));
1254: PetscCall(PetscSFDestroy(&preCoarseToFine));
1255: preCoarseToFine = preCoarseToFineNew;
1256: preCoarseToFine = preCoarseToFineNew;
1257: }
1258: if (coarseToPreFine) {
1259: PetscSF coarseToPreFineNew;
1260: PetscInt nleaves, nroots, i, nleavesCellSF, nleavesExpanded, *leavesNew;
1261: const PetscInt *leaves;
1262: const PetscSFNode *remotes;
1263: PetscSFNode *remotesNew, *remotesNewRoot, *remotesExpanded;
1265: PetscCall(PetscSFSetUp(coarseToPreFine));
1266: PetscCall(PetscSFGetGraph(coarseToPreFine, &nroots, &nleaves, &leaves, &remotes));
1267: PetscCall(PetscSFGetGraph(preCellSF, NULL, &nleavesCellSF, NULL, NULL));
1268: PetscCall(PetscMalloc1(nroots, &remotesNewRoot));
1269: PetscCall(PetscMalloc1(nleaves, &remotesNew));
1270: for (i = 0; i < nroots; i++) {
1271: remotesNewRoot[i].rank = rank;
1272: remotesNewRoot[i].index = i + cLocalStart;
1273: }
1274: PetscCall(PetscSFBcastBegin(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE));
1275: PetscCall(PetscSFBcastEnd(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE));
1276: PetscCall(PetscFree(remotesNewRoot));
1277: PetscCall(PetscMalloc1(nleavesCellSF, &remotesExpanded));
1278: for (i = 0; i < nleavesCellSF; i++) {
1279: remotesExpanded[i].rank = -1;
1280: remotesExpanded[i].index = -1;
1281: }
1282: for (i = 0; i < nleaves; i++) remotesExpanded[leaves ? leaves[i] : i] = remotesNew[i];
1283: PetscCall(PetscFree(remotesNew));
1284: PetscCall(PetscSFBcastBegin(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE));
1285: PetscCall(PetscSFBcastEnd(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE));
1287: nleavesExpanded = 0;
1288: for (i = 0; i < nleavesCellSF; i++) {
1289: if (remotesExpanded[i].rank >= 0) nleavesExpanded++;
1290: }
1291: PetscCall(PetscMalloc1(nleavesExpanded, &leavesNew));
1292: nleavesExpanded = 0;
1293: for (i = 0; i < nleavesCellSF; i++) {
1294: if (remotesExpanded[i].rank >= 0) {
1295: leavesNew[nleavesExpanded] = i;
1296: if (i > nleavesExpanded) remotesExpanded[nleavesExpanded] = remotes[i];
1297: nleavesExpanded++;
1298: }
1299: }
1300: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &coarseToPreFineNew));
1301: if (nleavesExpanded < nleavesCellSF) {
1302: PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, leavesNew, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES));
1303: } else {
1304: PetscCall(PetscFree(leavesNew));
1305: PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, NULL, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES));
1306: }
1307: PetscCall(PetscFree(remotesExpanded));
1308: PetscCall(PetscSFDestroy(&coarseToPreFine));
1309: coarseToPreFine = coarseToPreFineNew;
1310: }
1311: }
1312: }
1313: }
1314: forest->preCoarseToFine = preCoarseToFine;
1315: forest->coarseToPreFine = coarseToPreFine;
1316: dm->setupcalled = PETSC_TRUE;
1317: PetscCall(MPIU_Allreduce(&ctx.anyChange, &pforest->adaptivitySuccess, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
1318: PetscCall(DMPforestGetPlex(dm, NULL));
1319: PetscFunctionReturn(PETSC_SUCCESS);
1320: }
1322: #define DMForestGetAdaptivitySuccess_pforest _append_pforest(DMForestGetAdaptivitySuccess)
1323: static PetscErrorCode DMForestGetAdaptivitySuccess_pforest(DM dm, PetscBool *success)
1324: {
1325: DM_Forest *forest;
1326: DM_Forest_pforest *pforest;
1328: PetscFunctionBegin;
1329: forest = (DM_Forest *)dm->data;
1330: pforest = (DM_Forest_pforest *)forest->data;
1331: *success = pforest->adaptivitySuccess;
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: #define DMView_ASCII_pforest _append_pforest(DMView_ASCII)
1336: static PetscErrorCode DMView_ASCII_pforest(PetscObject odm, PetscViewer viewer)
1337: {
1338: DM dm = (DM)odm;
1340: PetscFunctionBegin;
1343: PetscCall(DMSetUp(dm));
1344: switch (viewer->format) {
1345: case PETSC_VIEWER_DEFAULT:
1346: case PETSC_VIEWER_ASCII_INFO: {
1347: PetscInt dim;
1348: const char *name;
1350: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1351: PetscCall(DMGetDimension(dm, &dim));
1352: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "Forest %s in %" PetscInt_FMT " dimensions:\n", name, dim));
1353: else PetscCall(PetscViewerASCIIPrintf(viewer, "Forest in %" PetscInt_FMT " dimensions:\n", dim));
1354: } /* fall through */
1355: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1356: case PETSC_VIEWER_LOAD_BALANCE: {
1357: DM plex;
1359: PetscCall(DMPforestGetPlex(dm, &plex));
1360: PetscCall(DMView(plex, viewer));
1361: } break;
1362: default:
1363: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1364: }
1365: PetscFunctionReturn(PETSC_SUCCESS);
1366: }
1368: #define DMView_VTK_pforest _append_pforest(DMView_VTK)
1369: static PetscErrorCode DMView_VTK_pforest(PetscObject odm, PetscViewer viewer)
1370: {
1371: DM dm = (DM)odm;
1372: DM_Forest *forest = (DM_Forest *)dm->data;
1373: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
1374: PetscBool isvtk;
1375: PetscReal vtkScale = 1. - PETSC_MACHINE_EPSILON;
1376: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
1377: const char *name;
1378: char *filenameStrip = NULL;
1379: PetscBool hasExt;
1380: size_t len;
1381: p4est_geometry_t *geom;
1383: PetscFunctionBegin;
1386: PetscCall(DMSetUp(dm));
1387: geom = pforest->topo->geom;
1388: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1389: PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
1390: switch (viewer->format) {
1391: case PETSC_VIEWER_VTK_VTU:
1392: PetscCheck(pforest->forest, PetscObjectComm(odm), PETSC_ERR_ARG_WRONG, "DM has not been setup with a valid forest");
1393: name = vtk->filename;
1394: PetscCall(PetscStrlen(name, &len));
1395: PetscCall(PetscStrcasecmp(name + len - 4, ".vtu", &hasExt));
1396: if (hasExt) {
1397: PetscCall(PetscStrallocpy(name, &filenameStrip));
1398: filenameStrip[len - 4] = '\0';
1399: name = filenameStrip;
1400: }
1401: if (!pforest->topo->geom) PetscCallP4estReturn(geom, p4est_geometry_new_connectivity, (pforest->topo->conn));
1402: {
1403: p4est_vtk_context_t *pvtk;
1404: int footerr;
1406: PetscCallP4estReturn(pvtk, p4est_vtk_context_new, (pforest->forest, name));
1407: PetscCallP4est(p4est_vtk_context_set_geom, (pvtk, geom));
1408: PetscCallP4est(p4est_vtk_context_set_scale, (pvtk, (double)vtkScale));
1409: PetscCallP4estReturn(pvtk, p4est_vtk_write_header, (pvtk));
1410: PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_header() failed");
1411: PetscCallP4estReturn(pvtk, p4est_vtk_write_cell_dataf,
1412: (pvtk, 1, /* write tree */
1413: 1, /* write level */
1414: 1, /* write rank */
1415: 0, /* do not wrap rank */
1416: 0, /* no scalar fields */
1417: 0, /* no vector fields */
1418: pvtk));
1419: PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_cell_dataf() failed");
1420: PetscCallP4estReturn(footerr, p4est_vtk_write_footer, (pvtk));
1421: PetscCheck(!footerr, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_footer() failed");
1422: }
1423: if (!pforest->topo->geom) PetscCallP4est(p4est_geometry_destroy, (geom));
1424: PetscCall(PetscFree(filenameStrip));
1425: break;
1426: default:
1427: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1428: }
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: #define DMView_HDF5_pforest _append_pforest(DMView_HDF5)
1433: static PetscErrorCode DMView_HDF5_pforest(DM dm, PetscViewer viewer)
1434: {
1435: DM plex;
1437: PetscFunctionBegin;
1438: PetscCall(DMSetUp(dm));
1439: PetscCall(DMPforestGetPlex(dm, &plex));
1440: PetscCall(DMView(plex, viewer));
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: #define DMView_GLVis_pforest _append_pforest(DMView_GLVis)
1445: static PetscErrorCode DMView_GLVis_pforest(DM dm, PetscViewer viewer)
1446: {
1447: DM plex;
1449: PetscFunctionBegin;
1450: PetscCall(DMSetUp(dm));
1451: PetscCall(DMPforestGetPlex(dm, &plex));
1452: PetscCall(DMView(plex, viewer));
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: #define DMView_pforest _append_pforest(DMView)
1457: static PetscErrorCode DMView_pforest(DM dm, PetscViewer viewer)
1458: {
1459: PetscBool isascii, isvtk, ishdf5, isglvis;
1461: PetscFunctionBegin;
1464: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1465: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1466: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1467: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
1468: if (isascii) {
1469: PetscCall(DMView_ASCII_pforest((PetscObject)dm, viewer));
1470: } else if (isvtk) {
1471: PetscCall(DMView_VTK_pforest((PetscObject)dm, viewer));
1472: } else if (ishdf5) {
1473: PetscCall(DMView_HDF5_pforest(dm, viewer));
1474: } else if (isglvis) {
1475: PetscCall(DMView_GLVis_pforest(dm, viewer));
1476: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Viewer not supported (not VTK, HDF5, or GLVis)");
1477: PetscFunctionReturn(PETSC_SUCCESS);
1478: }
1480: static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *conn, PetscInt **tree_face_to_uniq)
1481: {
1482: PetscInt *ttf, f, t, g, count;
1483: PetscInt numFacets;
1485: PetscFunctionBegin;
1486: numFacets = conn->num_trees * P4EST_FACES;
1487: PetscCall(PetscMalloc1(numFacets, &ttf));
1488: for (f = 0; f < numFacets; f++) ttf[f] = -1;
1489: for (g = 0, count = 0, t = 0; t < conn->num_trees; t++) {
1490: for (f = 0; f < P4EST_FACES; f++, g++) {
1491: if (ttf[g] == -1) {
1492: PetscInt ng;
1494: ttf[g] = count++;
1495: ng = conn->tree_to_tree[g] * P4EST_FACES + (conn->tree_to_face[g] % P4EST_FACES);
1496: ttf[ng] = ttf[g];
1497: }
1498: }
1499: }
1500: *tree_face_to_uniq = ttf;
1501: PetscFunctionReturn(PETSC_SUCCESS);
1502: }
1504: static PetscErrorCode DMPlexCreateConnectivity_pforest(DM dm, p4est_connectivity_t **connOut, PetscInt **tree_face_to_uniq)
1505: {
1506: p4est_topidx_t numTrees, numVerts, numCorns, numCtt;
1507: PetscSection ctt;
1508: #if defined(P4_TO_P8)
1509: p4est_topidx_t numEdges, numEtt;
1510: PetscSection ett;
1511: PetscInt eStart, eEnd, e, ettSize;
1512: PetscInt vertOff = 1 + P4EST_FACES + P8EST_EDGES;
1513: PetscInt edgeOff = 1 + P4EST_FACES;
1514: #else
1515: PetscInt vertOff = 1 + P4EST_FACES;
1516: #endif
1517: p4est_connectivity_t *conn;
1518: PetscInt cStart, cEnd, c, vStart, vEnd, v, fStart, fEnd, f;
1519: PetscInt *star = NULL, *closure = NULL, closureSize, starSize, cttSize;
1520: PetscInt *ttf;
1522: PetscFunctionBegin;
1523: /* 1: count objects, allocate */
1524: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1525: PetscCall(P4estTopidxCast(cEnd - cStart, &numTrees));
1526: numVerts = P4EST_CHILDREN * numTrees;
1527: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1528: PetscCall(P4estTopidxCast(vEnd - vStart, &numCorns));
1529: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ctt));
1530: PetscCall(PetscSectionSetChart(ctt, vStart, vEnd));
1531: for (v = vStart; v < vEnd; v++) {
1532: PetscInt s;
1534: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1535: for (s = 0; s < starSize; s++) {
1536: PetscInt p = star[2 * s];
1538: if (p >= cStart && p < cEnd) {
1539: /* we want to count every time cell p references v, so we see how many times it comes up in the closure. This
1540: * only protects against periodicity problems */
1541: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1542: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " with wrong closure size %" PetscInt_FMT " != %d", p, closureSize, P4EST_INSUL);
1543: for (c = 0; c < P4EST_CHILDREN; c++) {
1544: PetscInt cellVert = closure[2 * (c + vertOff)];
1546: PetscCheck(cellVert >= vStart && cellVert < vEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: vertices");
1547: if (cellVert == v) PetscCall(PetscSectionAddDof(ctt, v, 1));
1548: }
1549: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1550: }
1551: }
1552: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1553: }
1554: PetscCall(PetscSectionSetUp(ctt));
1555: PetscCall(PetscSectionGetStorageSize(ctt, &cttSize));
1556: PetscCall(P4estTopidxCast(cttSize, &numCtt));
1557: #if defined(P4_TO_P8)
1558: PetscCall(DMPlexGetSimplexOrBoxCells(dm, P4EST_DIM - 1, &eStart, &eEnd));
1559: PetscCall(P4estTopidxCast(eEnd - eStart, &numEdges));
1560: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ett));
1561: PetscCall(PetscSectionSetChart(ett, eStart, eEnd));
1562: for (e = eStart; e < eEnd; e++) {
1563: PetscInt s;
1565: PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1566: for (s = 0; s < starSize; s++) {
1567: PetscInt p = star[2 * s];
1569: if (p >= cStart && p < cEnd) {
1570: /* we want to count every time cell p references e, so we see how many times it comes up in the closure. This
1571: * only protects against periodicity problems */
1572: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1573: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell with wrong closure size");
1574: for (c = 0; c < P8EST_EDGES; c++) {
1575: PetscInt cellEdge = closure[2 * (c + edgeOff)];
1577: PetscCheck(cellEdge >= eStart && cellEdge < eEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: edges");
1578: if (cellEdge == e) PetscCall(PetscSectionAddDof(ett, e, 1));
1579: }
1580: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1581: }
1582: }
1583: PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1584: }
1585: PetscCall(PetscSectionSetUp(ett));
1586: PetscCall(PetscSectionGetStorageSize(ett, &ettSize));
1587: PetscCall(P4estTopidxCast(ettSize, &numEtt));
1589: /* This routine allocates space for the arrays, which we fill below */
1590: PetscCallP4estReturn(conn, p8est_connectivity_new, (numVerts, numTrees, numEdges, numEtt, numCorns, numCtt));
1591: #else
1592: PetscCallP4estReturn(conn, p4est_connectivity_new, (numVerts, numTrees, numCorns, numCtt));
1593: #endif
1595: /* 2: visit every face, determine neighboring cells(trees) */
1596: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 1, &fStart, &fEnd));
1597: PetscCall(PetscMalloc1((cEnd - cStart) * P4EST_FACES, &ttf));
1598: for (f = fStart; f < fEnd; f++) {
1599: PetscInt numSupp, s;
1600: PetscInt myFace[2] = {-1, -1};
1601: PetscInt myOrnt[2] = {PETSC_MIN_INT, PETSC_MIN_INT};
1602: const PetscInt *supp;
1604: PetscCall(DMPlexGetSupportSize(dm, f, &numSupp));
1605: PetscCheck(numSupp == 1 || numSupp == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "point %" PetscInt_FMT " has facet with %" PetscInt_FMT " sides: must be 1 or 2 (boundary or conformal)", f, numSupp);
1606: PetscCall(DMPlexGetSupport(dm, f, &supp));
1608: for (s = 0; s < numSupp; s++) {
1609: PetscInt p = supp[s];
1611: if (p >= cEnd) {
1612: numSupp--;
1613: if (s) supp = &supp[1 - s];
1614: break;
1615: }
1616: }
1617: for (s = 0; s < numSupp; s++) {
1618: PetscInt p = supp[s], i;
1619: PetscInt numCone;
1620: DMPolytopeType ct;
1621: const PetscInt *cone;
1622: const PetscInt *ornt;
1623: PetscInt orient = PETSC_MIN_INT;
1625: PetscCall(DMPlexGetConeSize(dm, p, &numCone));
1626: PetscCheck(numCone == P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " has %" PetscInt_FMT " facets, expect %d", p, numCone, P4EST_FACES);
1627: PetscCall(DMPlexGetCone(dm, p, &cone));
1628: PetscCall(DMPlexGetCellType(dm, cone[0], &ct));
1629: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
1630: for (i = 0; i < P4EST_FACES; i++) {
1631: if (cone[i] == f) {
1632: orient = DMPolytopeConvertNewOrientation_Internal(ct, ornt[i]);
1633: break;
1634: }
1635: }
1636: PetscCheck(i < P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " faced %" PetscInt_FMT " mismatch", p, f);
1637: if (p < cStart || p >= cEnd) {
1638: DMPolytopeType ct;
1639: PetscCall(DMPlexGetCellType(dm, p, &ct));
1640: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " (%s) should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], cStart, cEnd);
1641: }
1642: ttf[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = f - fStart;
1643: if (numSupp == 1) {
1644: /* boundary faces indicated by self reference */
1645: conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = p - cStart;
1646: conn->tree_to_face[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = (int8_t)PetscFaceToP4estFace[i];
1647: } else {
1648: const PetscInt N = P4EST_CHILDREN / 2;
1650: conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = supp[1 - s] - cStart;
1651: myFace[s] = PetscFaceToP4estFace[i];
1652: /* get the orientation of cell p in p4est-type closure to facet f, by composing the p4est-closure to
1653: * petsc-closure permutation and the petsc-closure to facet orientation */
1654: myOrnt[s] = DihedralCompose(N, orient, DMPolytopeConvertNewOrientation_Internal(ct, P4estFaceToPetscOrnt[myFace[s]]));
1655: }
1656: }
1657: if (numSupp == 2) {
1658: for (s = 0; s < numSupp; s++) {
1659: PetscInt p = supp[s];
1660: PetscInt orntAtoB;
1661: PetscInt p4estOrient;
1662: const PetscInt N = P4EST_CHILDREN / 2;
1664: /* composing the forward permutation with the other cell's inverse permutation gives the self-to-neighbor
1665: * permutation of this cell-facet's cone */
1666: orntAtoB = DihedralCompose(N, DihedralInvert(N, myOrnt[1 - s]), myOrnt[s]);
1668: /* convert cone-description permutation (i.e., edges around facet) to cap-description permutation (i.e.,
1669: * vertices around facet) */
1670: #if !defined(P4_TO_P8)
1671: p4estOrient = orntAtoB < 0 ? -(orntAtoB + 1) : orntAtoB;
1672: #else
1673: {
1674: PetscInt firstVert = orntAtoB < 0 ? ((-orntAtoB) % N) : orntAtoB;
1675: PetscInt p4estFirstVert = firstVert < 2 ? firstVert : (firstVert ^ 1);
1677: /* swap bits */
1678: p4estOrient = ((myFace[s] <= myFace[1 - s]) || (orntAtoB < 0)) ? p4estFirstVert : ((p4estFirstVert >> 1) | ((p4estFirstVert & 1) << 1));
1679: }
1680: #endif
1681: /* encode neighbor face and orientation in tree_to_face per p4est_connectivity standard (see
1682: * p4est_connectivity.h, p8est_connectivity.h) */
1683: conn->tree_to_face[P4EST_FACES * (p - cStart) + myFace[s]] = (int8_t)myFace[1 - s] + p4estOrient * P4EST_FACES;
1684: }
1685: }
1686: }
1688: #if defined(P4_TO_P8)
1689: /* 3: visit every edge */
1690: conn->ett_offset[0] = 0;
1691: for (e = eStart; e < eEnd; e++) {
1692: PetscInt off, s;
1694: PetscCall(PetscSectionGetOffset(ett, e, &off));
1695: conn->ett_offset[e - eStart] = (p4est_topidx_t)off;
1696: PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1697: for (s = 0; s < starSize; s++) {
1698: PetscInt p = star[2 * s];
1700: if (p >= cStart && p < cEnd) {
1701: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1702: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure");
1703: for (c = 0; c < P8EST_EDGES; c++) {
1704: PetscInt cellEdge = closure[2 * (c + edgeOff)];
1705: PetscInt cellOrnt = closure[2 * (c + edgeOff) + 1];
1706: DMPolytopeType ct;
1708: PetscCall(DMPlexGetCellType(dm, cellEdge, &ct));
1709: cellOrnt = DMPolytopeConvertNewOrientation_Internal(ct, cellOrnt);
1710: if (cellEdge == e) {
1711: PetscInt p4estEdge = PetscEdgeToP4estEdge[c];
1712: PetscInt totalOrient;
1714: /* compose p4est-closure to petsc-closure permutation and petsc-closure to edge orientation */
1715: totalOrient = DihedralCompose(2, cellOrnt, DMPolytopeConvertNewOrientation_Internal(DM_POLYTOPE_SEGMENT, P4estEdgeToPetscOrnt[p4estEdge]));
1716: /* p4est orientations are positive: -2 => 1, -1 => 0 */
1717: totalOrient = (totalOrient < 0) ? -(totalOrient + 1) : totalOrient;
1718: conn->edge_to_tree[off] = (p4est_locidx_t)(p - cStart);
1719: /* encode cell-edge and orientation in edge_to_edge per p8est_connectivity standard (see
1720: * p8est_connectivity.h) */
1721: conn->edge_to_edge[off++] = (int8_t)p4estEdge + P8EST_EDGES * totalOrient;
1722: conn->tree_to_edge[P8EST_EDGES * (p - cStart) + p4estEdge] = e - eStart;
1723: }
1724: }
1725: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1726: }
1727: }
1728: PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1729: }
1730: PetscCall(PetscSectionDestroy(&ett));
1731: #endif
1733: /* 4: visit every vertex */
1734: conn->ctt_offset[0] = 0;
1735: for (v = vStart; v < vEnd; v++) {
1736: PetscInt off, s;
1738: PetscCall(PetscSectionGetOffset(ctt, v, &off));
1739: conn->ctt_offset[v - vStart] = (p4est_topidx_t)off;
1740: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1741: for (s = 0; s < starSize; s++) {
1742: PetscInt p = star[2 * s];
1744: if (p >= cStart && p < cEnd) {
1745: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1746: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure");
1747: for (c = 0; c < P4EST_CHILDREN; c++) {
1748: PetscInt cellVert = closure[2 * (c + vertOff)];
1750: if (cellVert == v) {
1751: PetscInt p4estVert = PetscVertToP4estVert[c];
1753: conn->corner_to_tree[off] = (p4est_locidx_t)(p - cStart);
1754: conn->corner_to_corner[off++] = (int8_t)p4estVert;
1755: conn->tree_to_corner[P4EST_CHILDREN * (p - cStart) + p4estVert] = v - vStart;
1756: }
1757: }
1758: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1759: }
1760: }
1761: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1762: }
1763: PetscCall(PetscSectionDestroy(&ctt));
1765: /* 5: Compute the coordinates */
1766: {
1767: PetscInt coordDim;
1769: PetscCall(DMGetCoordinateDim(dm, &coordDim));
1770: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1771: for (c = cStart; c < cEnd; c++) {
1772: PetscInt dof;
1773: PetscBool isDG;
1774: PetscScalar *cellCoords = NULL;
1775: const PetscScalar *array;
1777: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords));
1778: PetscCheck(dof == P4EST_CHILDREN * coordDim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need coordinates at the corners: (dof) %" PetscInt_FMT " != %d * %" PetscInt_FMT " (sdim)", dof, P4EST_CHILDREN, coordDim);
1779: for (v = 0; v < P4EST_CHILDREN; v++) {
1780: PetscInt i, lim = PetscMin(3, coordDim);
1781: PetscInt p4estVert = PetscVertToP4estVert[v];
1783: conn->tree_to_vertex[P4EST_CHILDREN * (c - cStart) + v] = P4EST_CHILDREN * (c - cStart) + v;
1784: /* p4est vertices are always embedded in R^3 */
1785: for (i = 0; i < 3; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = 0.;
1786: for (i = 0; i < lim; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = PetscRealPart(cellCoords[v * coordDim + i]);
1787: }
1788: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords));
1789: }
1790: }
1792: #if defined(P4EST_ENABLE_DEBUG)
1793: PetscCheck(p4est_connectivity_is_valid(conn), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Plex to p4est conversion failed");
1794: #endif
1796: *connOut = conn;
1798: *tree_face_to_uniq = ttf;
1800: PetscFunctionReturn(PETSC_SUCCESS);
1801: }
1803: static PetscErrorCode locidx_to_PetscInt(sc_array_t *array)
1804: {
1805: sc_array_t *newarray;
1806: size_t zz, count = array->elem_count;
1808: PetscFunctionBegin;
1809: PetscCheck(array->elem_size == sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size");
1811: if (sizeof(p4est_locidx_t) == sizeof(PetscInt)) PetscFunctionReturn(PETSC_SUCCESS);
1813: newarray = sc_array_new_size(sizeof(PetscInt), array->elem_count);
1814: for (zz = 0; zz < count; zz++) {
1815: p4est_locidx_t il = *((p4est_locidx_t *)sc_array_index(array, zz));
1816: PetscInt *ip = (PetscInt *)sc_array_index(newarray, zz);
1818: *ip = (PetscInt)il;
1819: }
1821: sc_array_reset(array);
1822: sc_array_init_size(array, sizeof(PetscInt), count);
1823: sc_array_copy(array, newarray);
1824: sc_array_destroy(newarray);
1825: PetscFunctionReturn(PETSC_SUCCESS);
1826: }
1828: static PetscErrorCode coords_double_to_PetscScalar(sc_array_t *array, PetscInt dim)
1829: {
1830: sc_array_t *newarray;
1831: size_t zz, count = array->elem_count;
1833: PetscFunctionBegin;
1834: PetscCheck(array->elem_size == 3 * sizeof(double), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong coordinate size");
1835: #if !defined(PETSC_USE_COMPLEX)
1836: if (sizeof(double) == sizeof(PetscScalar) && dim == 3) PetscFunctionReturn(PETSC_SUCCESS);
1837: #endif
1839: newarray = sc_array_new_size(dim * sizeof(PetscScalar), array->elem_count);
1840: for (zz = 0; zz < count; zz++) {
1841: int i;
1842: double *id = (double *)sc_array_index(array, zz);
1843: PetscScalar *ip = (PetscScalar *)sc_array_index(newarray, zz);
1845: for (i = 0; i < dim; i++) ip[i] = 0.;
1846: for (i = 0; i < PetscMin(dim, 3); i++) ip[i] = (PetscScalar)id[i];
1847: }
1849: sc_array_reset(array);
1850: sc_array_init_size(array, dim * sizeof(PetscScalar), count);
1851: sc_array_copy(array, newarray);
1852: sc_array_destroy(newarray);
1853: PetscFunctionReturn(PETSC_SUCCESS);
1854: }
1856: static PetscErrorCode locidx_pair_to_PetscSFNode(sc_array_t *array)
1857: {
1858: sc_array_t *newarray;
1859: size_t zz, count = array->elem_count;
1861: PetscFunctionBegin;
1862: PetscCheck(array->elem_size == 2 * sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size");
1864: newarray = sc_array_new_size(sizeof(PetscSFNode), array->elem_count);
1865: for (zz = 0; zz < count; zz++) {
1866: p4est_locidx_t *il = (p4est_locidx_t *)sc_array_index(array, zz);
1867: PetscSFNode *ip = (PetscSFNode *)sc_array_index(newarray, zz);
1869: ip->rank = (PetscInt)il[0];
1870: ip->index = (PetscInt)il[1];
1871: }
1873: sc_array_reset(array);
1874: sc_array_init_size(array, sizeof(PetscSFNode), count);
1875: sc_array_copy(array, newarray);
1876: sc_array_destroy(newarray);
1877: PetscFunctionReturn(PETSC_SUCCESS);
1878: }
1880: static PetscErrorCode P4estToPlex_Local(p4est_t *p4est, DM *plex)
1881: {
1882: PetscFunctionBegin;
1883: {
1884: sc_array_t *points_per_dim = sc_array_new(sizeof(p4est_locidx_t));
1885: sc_array_t *cone_sizes = sc_array_new(sizeof(p4est_locidx_t));
1886: sc_array_t *cones = sc_array_new(sizeof(p4est_locidx_t));
1887: sc_array_t *cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
1888: sc_array_t *coords = sc_array_new(3 * sizeof(double));
1889: sc_array_t *children = sc_array_new(sizeof(p4est_locidx_t));
1890: sc_array_t *parents = sc_array_new(sizeof(p4est_locidx_t));
1891: sc_array_t *childids = sc_array_new(sizeof(p4est_locidx_t));
1892: sc_array_t *leaves = sc_array_new(sizeof(p4est_locidx_t));
1893: sc_array_t *remotes = sc_array_new(2 * sizeof(p4est_locidx_t));
1894: p4est_locidx_t first_local_quad;
1896: 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));
1898: PetscCall(locidx_to_PetscInt(points_per_dim));
1899: PetscCall(locidx_to_PetscInt(cone_sizes));
1900: PetscCall(locidx_to_PetscInt(cones));
1901: PetscCall(locidx_to_PetscInt(cone_orientations));
1902: PetscCall(coords_double_to_PetscScalar(coords, P4EST_DIM));
1904: PetscCall(DMPlexCreate(PETSC_COMM_SELF, plex));
1905: PetscCall(DMSetDimension(*plex, P4EST_DIM));
1906: PetscCall(DMPlexCreateFromDAG(*plex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array));
1907: PetscCall(DMPlexConvertOldOrientations_Internal(*plex));
1908: sc_array_destroy(points_per_dim);
1909: sc_array_destroy(cone_sizes);
1910: sc_array_destroy(cones);
1911: sc_array_destroy(cone_orientations);
1912: sc_array_destroy(coords);
1913: sc_array_destroy(children);
1914: sc_array_destroy(parents);
1915: sc_array_destroy(childids);
1916: sc_array_destroy(leaves);
1917: sc_array_destroy(remotes);
1918: }
1919: PetscFunctionReturn(PETSC_SUCCESS);
1920: }
1922: #define DMReferenceTreeGetChildSymmetry_pforest _append_pforest(DMReferenceTreeGetChildSymmetry)
1923: static PetscErrorCode DMReferenceTreeGetChildSymmetry_pforest(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
1924: {
1925: PetscInt coneSize, dStart, dEnd, vStart, vEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
1927: PetscFunctionBegin;
1928: if (parentOrientA == parentOrientB) {
1929: if (childOrientB) *childOrientB = childOrientA;
1930: if (childB) *childB = childA;
1931: PetscFunctionReturn(PETSC_SUCCESS);
1932: }
1933: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1934: if (childA >= vStart && childA < vEnd) { /* vertices (always in the middle) are invariant under rotation */
1935: if (childOrientB) *childOrientB = 0;
1936: if (childB) *childB = childA;
1937: PetscFunctionReturn(PETSC_SUCCESS);
1938: }
1939: for (dim = 0; dim < 3; dim++) {
1940: PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd));
1941: if (parent >= dStart && parent <= dEnd) break;
1942: }
1943: PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim);
1944: PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children");
1945: if (childA < dStart || childA >= dEnd) { /* a 1-cell in a 2-cell */
1946: /* this is a lower-dimensional child: bootstrap */
1947: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
1948: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
1950: PetscCall(DMPlexGetSupportSize(dm, childA, &size));
1951: PetscCall(DMPlexGetSupport(dm, childA, &supp));
1953: /* find a point sA in supp(childA) that has the same parent */
1954: for (i = 0; i < size; i++) {
1955: PetscInt sParent;
1957: sA = supp[i];
1958: if (sA == parent) continue;
1959: PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL));
1960: if (sParent == parent) break;
1961: }
1962: PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children");
1963: /* find out which point sB is in an equivalent position to sA under
1964: * parentOrientB */
1965: PetscCall(DMReferenceTreeGetChildSymmetry_pforest(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB));
1966: PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize));
1967: PetscCall(DMPlexGetCone(dm, sA, &coneA));
1968: PetscCall(DMPlexGetCone(dm, sB, &coneB));
1969: PetscCall(DMPlexGetConeOrientation(dm, sA, &oA));
1970: PetscCall(DMPlexGetConeOrientation(dm, sB, &oB));
1971: /* step through the cone of sA in natural order */
1972: for (i = 0; i < sConeSize; i++) {
1973: if (coneA[i] == childA) {
1974: /* if childA is at position i in coneA,
1975: * then we want the point that is at sOrientB*i in coneB */
1976: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
1977: if (childB) *childB = coneB[j];
1978: if (childOrientB) {
1979: DMPolytopeType ct;
1980: PetscInt oBtrue;
1982: PetscCall(DMPlexGetConeSize(dm, childA, &coneSize));
1983: /* compose sOrientB and oB[j] */
1984: PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge");
1985: ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1986: /* we may have to flip an edge */
1987: oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientationInv(ct, -1, oB[j]);
1988: oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
1989: ABswap = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
1990: *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
1991: }
1992: break;
1993: }
1994: }
1995: PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch");
1996: PetscFunctionReturn(PETSC_SUCCESS);
1997: }
1998: /* get the cone size and symmetry swap */
1999: PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
2000: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
2001: if (dim == 2) {
2002: /* orientations refer to cones: we want them to refer to vertices:
2003: * if it's a rotation, they are the same, but if the order is reversed, a
2004: * permutation that puts side i first does *not* put vertex i first */
2005: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
2006: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
2007: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
2008: } else {
2009: oAvert = parentOrientA;
2010: oBvert = parentOrientB;
2011: ABswapVert = ABswap;
2012: }
2013: if (childB) {
2014: /* assume that each child corresponds to a vertex, in the same order */
2015: PetscInt p, posA = -1, numChildren, i;
2016: const PetscInt *children;
2018: /* count which position the child is in */
2019: PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children));
2020: for (i = 0; i < numChildren; i++) {
2021: p = children[i];
2022: if (p == childA) {
2023: if (dim == 1) {
2024: posA = i;
2025: } else { /* 2D Morton to rotation */
2026: posA = (i & 2) ? (i ^ 1) : i;
2027: }
2028: break;
2029: }
2030: }
2031: if (posA >= coneSize) {
2032: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find childA in children of parent");
2033: } else {
2034: /* figure out position B by applying ABswapVert */
2035: PetscInt posB, childIdB;
2037: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
2038: if (dim == 1) {
2039: childIdB = posB;
2040: } else { /* 2D rotation to Morton */
2041: childIdB = (posB & 2) ? (posB ^ 1) : posB;
2042: }
2043: if (childB) *childB = children[childIdB];
2044: }
2045: }
2046: if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
2047: PetscFunctionReturn(PETSC_SUCCESS);
2048: }
2050: #define DMCreateReferenceTree_pforest _append_pforest(DMCreateReferenceTree)
2051: static PetscErrorCode DMCreateReferenceTree_pforest(MPI_Comm comm, DM *dm)
2052: {
2053: p4est_connectivity_t *refcube;
2054: p4est_t *root, *refined;
2055: DM dmRoot, dmRefined;
2056: DM_Plex *mesh;
2057: PetscMPIInt rank;
2058: #if defined(PETSC_HAVE_MPIUNI)
2059: sc_MPI_Comm comm_self = sc_MPI_COMM_SELF;
2060: #else
2061: MPI_Comm comm_self = PETSC_COMM_SELF;
2062: #endif
2064: PetscFunctionBegin;
2065: PetscCallP4estReturn(refcube, p4est_connectivity_new_byname, ("unit"));
2066: { /* [-1,1]^d geometry */
2067: PetscInt i, j;
2069: for (i = 0; i < P4EST_CHILDREN; i++) {
2070: for (j = 0; j < 3; j++) {
2071: refcube->vertices[3 * i + j] *= 2.;
2072: refcube->vertices[3 * i + j] -= 1.;
2073: }
2074: }
2075: }
2076: PetscCallP4estReturn(root, p4est_new, (comm_self, refcube, 0, NULL, NULL));
2077: PetscCallP4estReturn(refined, p4est_new_ext, (comm_self, refcube, 0, 1, 1, 0, NULL, NULL));
2078: PetscCall(P4estToPlex_Local(root, &dmRoot));
2079: PetscCall(P4estToPlex_Local(refined, &dmRefined));
2080: {
2081: #if !defined(P4_TO_P8)
2082: PetscInt nPoints = 25;
2083: 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};
2084: 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};
2085: #else
2086: PetscInt nPoints = 125;
2087: 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,
2088: 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,
2089: 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};
2090: 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,
2091: 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};
2093: #endif
2094: IS permIS;
2095: DM dmPerm;
2097: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nPoints, perm, PETSC_USE_POINTER, &permIS));
2098: PetscCall(DMPlexPermute(dmRefined, permIS, &dmPerm));
2099: if (dmPerm) {
2100: PetscCall(DMDestroy(&dmRefined));
2101: dmRefined = dmPerm;
2102: }
2103: PetscCall(ISDestroy(&permIS));
2104: {
2105: PetscInt p;
2106: PetscCall(DMCreateLabel(dmRoot, "identity"));
2107: PetscCall(DMCreateLabel(dmRefined, "identity"));
2108: for (p = 0; p < P4EST_INSUL; p++) PetscCall(DMSetLabelValue(dmRoot, "identity", p, p));
2109: for (p = 0; p < nPoints; p++) PetscCall(DMSetLabelValue(dmRefined, "identity", p, ident[p]));
2110: }
2111: }
2112: PetscCall(DMPlexCreateReferenceTree_Union(dmRoot, dmRefined, "identity", dm));
2113: mesh = (DM_Plex *)(*dm)->data;
2114: mesh->getchildsymmetry = DMReferenceTreeGetChildSymmetry_pforest;
2115: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2116: if (rank == 0) {
2117: PetscCall(DMViewFromOptions(dmRoot, NULL, "-dm_p4est_ref_root_view"));
2118: PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_refined_view"));
2119: PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_tree_view"));
2120: }
2121: PetscCall(DMDestroy(&dmRefined));
2122: PetscCall(DMDestroy(&dmRoot));
2123: PetscCallP4est(p4est_destroy, (refined));
2124: PetscCallP4est(p4est_destroy, (root));
2125: PetscCallP4est(p4est_connectivity_destroy, (refcube));
2126: PetscFunctionReturn(PETSC_SUCCESS);
2127: }
2129: static PetscErrorCode DMShareDiscretization(DM dmA, DM dmB)
2130: {
2131: void *ctx;
2132: PetscInt num;
2133: PetscReal val;
2135: PetscFunctionBegin;
2136: PetscCall(DMGetApplicationContext(dmA, &ctx));
2137: PetscCall(DMSetApplicationContext(dmB, ctx));
2138: PetscCall(DMCopyDisc(dmA, dmB));
2139: PetscCall(DMGetOutputSequenceNumber(dmA, &num, &val));
2140: PetscCall(DMSetOutputSequenceNumber(dmB, num, val));
2141: if (dmB->localSection != dmA->localSection || dmB->globalSection != dmA->globalSection) {
2142: PetscCall(DMClearLocalVectors(dmB));
2143: PetscCall(PetscObjectReference((PetscObject)dmA->localSection));
2144: PetscCall(PetscSectionDestroy(&(dmB->localSection)));
2145: dmB->localSection = dmA->localSection;
2146: PetscCall(DMClearGlobalVectors(dmB));
2147: PetscCall(PetscObjectReference((PetscObject)dmA->globalSection));
2148: PetscCall(PetscSectionDestroy(&(dmB->globalSection)));
2149: dmB->globalSection = dmA->globalSection;
2150: PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.section));
2151: PetscCall(PetscSectionDestroy(&(dmB->defaultConstraint.section)));
2152: dmB->defaultConstraint.section = dmA->defaultConstraint.section;
2153: PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.mat));
2154: PetscCall(MatDestroy(&(dmB->defaultConstraint.mat)));
2155: dmB->defaultConstraint.mat = dmA->defaultConstraint.mat;
2156: if (dmA->map) PetscCall(PetscLayoutReference(dmA->map, &dmB->map));
2157: }
2158: if (dmB->sectionSF != dmA->sectionSF) {
2159: PetscCall(PetscObjectReference((PetscObject)dmA->sectionSF));
2160: PetscCall(PetscSFDestroy(&dmB->sectionSF));
2161: dmB->sectionSF = dmA->sectionSF;
2162: }
2163: PetscFunctionReturn(PETSC_SUCCESS);
2164: }
2166: /* Get an SF that broadcasts a coarse-cell covering of the local fine cells */
2167: static PetscErrorCode DMPforestGetCellCoveringSF(MPI_Comm comm, p4est_t *p4estC, p4est_t *p4estF, PetscInt cStart, PetscInt cEnd, PetscSF *coveringSF)
2168: {
2169: PetscInt startF, endF, startC, endC, p, nLeaves;
2170: PetscSFNode *leaves;
2171: PetscSF sf;
2172: PetscInt *recv, *send;
2173: PetscMPIInt tag;
2174: MPI_Request *recvReqs, *sendReqs;
2175: PetscSection section;
2177: PetscFunctionBegin;
2178: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estF, p4estC, &startC, &endC));
2179: PetscCall(PetscMalloc2(2 * (endC - startC), &recv, endC - startC, &recvReqs));
2180: PetscCall(PetscCommGetNewTag(comm, &tag));
2181: for (p = startC; p < endC; p++) {
2182: recvReqs[p - startC] = MPI_REQUEST_NULL; /* just in case we don't initiate a receive */
2183: if (p4estC->global_first_quadrant[p] == p4estC->global_first_quadrant[p + 1]) { /* empty coarse partition */
2184: recv[2 * (p - startC)] = 0;
2185: recv[2 * (p - startC) + 1] = 0;
2186: continue;
2187: }
2189: PetscCallMPI(MPI_Irecv(&recv[2 * (p - startC)], 2, MPIU_INT, p, tag, comm, &recvReqs[p - startC]));
2190: }
2191: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estC, p4estF, &startF, &endF));
2192: PetscCall(PetscMalloc2(2 * (endF - startF), &send, endF - startF, &sendReqs));
2193: /* count the quadrants rank will send to each of [startF,endF) */
2194: for (p = startF; p < endF; p++) {
2195: p4est_quadrant_t *myFineStart = &p4estF->global_first_position[p];
2196: p4est_quadrant_t *myFineEnd = &p4estF->global_first_position[p + 1];
2197: PetscInt tStart = (PetscInt)myFineStart->p.which_tree;
2198: PetscInt tEnd = (PetscInt)myFineEnd->p.which_tree;
2199: PetscInt firstCell = -1, lastCell = -1;
2200: p4est_tree_t *treeStart = &(((p4est_tree_t *)p4estC->trees->array)[tStart]);
2201: p4est_tree_t *treeEnd = (size_t)tEnd < p4estC->trees->elem_count ? &(((p4est_tree_t *)p4estC->trees->array)[tEnd]) : NULL;
2202: ssize_t overlapIndex;
2204: sendReqs[p - startF] = MPI_REQUEST_NULL; /* just in case we don't initiate a send */
2205: if (p4estF->global_first_quadrant[p] == p4estF->global_first_quadrant[p + 1]) continue;
2207: /* locate myFineStart in (or before) a cell */
2208: if (treeStart->quadrants.elem_count) {
2209: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&(treeStart->quadrants), myFineStart, p4est_quadrant_disjoint));
2210: if (overlapIndex < 0) {
2211: firstCell = 0;
2212: } else {
2213: firstCell = treeStart->quadrants_offset + overlapIndex;
2214: }
2215: } else {
2216: firstCell = 0;
2217: }
2218: if (treeEnd && treeEnd->quadrants.elem_count) {
2219: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&(treeEnd->quadrants), myFineEnd, p4est_quadrant_disjoint));
2220: if (overlapIndex < 0) { /* all of this local section is overlapped */
2221: lastCell = p4estC->local_num_quadrants;
2222: } else {
2223: p4est_quadrant_t *container = &(((p4est_quadrant_t *)treeEnd->quadrants.array)[overlapIndex]);
2224: p4est_quadrant_t first_desc;
2225: int equal;
2227: PetscCallP4est(p4est_quadrant_first_descendant, (container, &first_desc, P4EST_QMAXLEVEL));
2228: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (myFineEnd, &first_desc));
2229: if (equal) {
2230: lastCell = treeEnd->quadrants_offset + overlapIndex;
2231: } else {
2232: lastCell = treeEnd->quadrants_offset + overlapIndex + 1;
2233: }
2234: }
2235: } else {
2236: lastCell = p4estC->local_num_quadrants;
2237: }
2238: send[2 * (p - startF)] = firstCell;
2239: send[2 * (p - startF) + 1] = lastCell - firstCell;
2240: PetscCallMPI(MPI_Isend(&send[2 * (p - startF)], 2, MPIU_INT, p, tag, comm, &sendReqs[p - startF]));
2241: }
2242: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endC - startC), recvReqs, MPI_STATUSES_IGNORE));
2243: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
2244: PetscCall(PetscSectionSetChart(section, startC, endC));
2245: for (p = startC; p < endC; p++) {
2246: PetscInt numCells = recv[2 * (p - startC) + 1];
2247: PetscCall(PetscSectionSetDof(section, p, numCells));
2248: }
2249: PetscCall(PetscSectionSetUp(section));
2250: PetscCall(PetscSectionGetStorageSize(section, &nLeaves));
2251: PetscCall(PetscMalloc1(nLeaves, &leaves));
2252: for (p = startC; p < endC; p++) {
2253: PetscInt firstCell = recv[2 * (p - startC)];
2254: PetscInt numCells = recv[2 * (p - startC) + 1];
2255: PetscInt off, i;
2257: PetscCall(PetscSectionGetOffset(section, p, &off));
2258: for (i = 0; i < numCells; i++) {
2259: leaves[off + i].rank = p;
2260: leaves[off + i].index = firstCell + i;
2261: }
2262: }
2263: PetscCall(PetscSFCreate(comm, &sf));
2264: PetscCall(PetscSFSetGraph(sf, cEnd - cStart, nLeaves, NULL, PETSC_OWN_POINTER, leaves, PETSC_OWN_POINTER));
2265: PetscCall(PetscSectionDestroy(§ion));
2266: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endF - startF), sendReqs, MPI_STATUSES_IGNORE));
2267: PetscCall(PetscFree2(send, sendReqs));
2268: PetscCall(PetscFree2(recv, recvReqs));
2269: *coveringSF = sf;
2270: PetscFunctionReturn(PETSC_SUCCESS);
2271: }
2273: /* closure points for locally-owned cells */
2274: static PetscErrorCode DMPforestGetCellSFNodes(DM dm, PetscInt numClosureIndices, PetscInt *numClosurePoints, PetscSFNode **closurePoints, PetscBool redirect)
2275: {
2276: PetscInt cStart, cEnd;
2277: PetscInt count, c;
2278: PetscMPIInt rank;
2279: PetscInt closureSize = -1;
2280: PetscInt *closure = NULL;
2281: PetscSF pointSF;
2282: PetscInt nleaves, nroots;
2283: const PetscInt *ilocal;
2284: const PetscSFNode *iremote;
2285: DM plex;
2286: DM_Forest *forest;
2287: DM_Forest_pforest *pforest;
2289: PetscFunctionBegin;
2290: forest = (DM_Forest *)dm->data;
2291: pforest = (DM_Forest_pforest *)forest->data;
2292: cStart = pforest->cLocalStart;
2293: cEnd = pforest->cLocalEnd;
2294: PetscCall(DMPforestGetPlex(dm, &plex));
2295: PetscCall(DMGetPointSF(dm, &pointSF));
2296: PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
2297: nleaves = PetscMax(0, nleaves);
2298: nroots = PetscMax(0, nroots);
2299: *numClosurePoints = numClosureIndices * (cEnd - cStart);
2300: PetscCall(PetscMalloc1(*numClosurePoints, closurePoints));
2301: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2302: for (c = cStart, count = 0; c < cEnd; c++) {
2303: PetscInt i;
2304: PetscCall(DMPlexGetTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2306: for (i = 0; i < numClosureIndices; i++, count++) {
2307: PetscInt p = closure[2 * i];
2308: PetscInt loc = -1;
2310: PetscCall(PetscFindInt(p, nleaves, ilocal, &loc));
2311: if (redirect && loc >= 0) {
2312: (*closurePoints)[count].rank = iremote[loc].rank;
2313: (*closurePoints)[count].index = iremote[loc].index;
2314: } else {
2315: (*closurePoints)[count].rank = rank;
2316: (*closurePoints)[count].index = p;
2317: }
2318: }
2319: PetscCall(DMPlexRestoreTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2320: }
2321: PetscFunctionReturn(PETSC_SUCCESS);
2322: }
2324: static void MPIAPI DMPforestMaxSFNode(void *a, void *b, PetscMPIInt *len, MPI_Datatype *type)
2325: {
2326: PetscMPIInt i;
2328: for (i = 0; i < *len; i++) {
2329: PetscSFNode *A = (PetscSFNode *)a;
2330: PetscSFNode *B = (PetscSFNode *)b;
2332: if (B->rank < 0) *B = *A;
2333: }
2334: }
2336: static PetscErrorCode DMPforestGetTransferSF_Point(DM coarse, DM fine, PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2337: {
2338: MPI_Comm comm;
2339: PetscMPIInt rank, size;
2340: DM_Forest_pforest *pforestC, *pforestF;
2341: p4est_t *p4estC, *p4estF;
2342: PetscInt numClosureIndices;
2343: PetscInt numClosurePointsC, numClosurePointsF;
2344: PetscSFNode *closurePointsC, *closurePointsF;
2345: p4est_quadrant_t *coverQuads = NULL;
2346: p4est_quadrant_t **treeQuads;
2347: PetscInt *treeQuadCounts;
2348: MPI_Datatype nodeType;
2349: MPI_Datatype nodeClosureType;
2350: MPI_Op sfNodeReduce;
2351: p4est_topidx_t fltF, lltF, t;
2352: DM plexC, plexF;
2353: PetscInt pStartF, pEndF, pStartC, pEndC;
2354: PetscBool saveInCoarse = PETSC_FALSE;
2355: PetscBool saveInFine = PETSC_FALSE;
2356: PetscBool formCids = (childIds != NULL) ? PETSC_TRUE : PETSC_FALSE;
2357: PetscInt *cids = NULL;
2359: PetscFunctionBegin;
2360: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2361: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2362: p4estC = pforestC->forest;
2363: p4estF = pforestF->forest;
2364: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2365: comm = PetscObjectComm((PetscObject)coarse);
2366: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2367: PetscCallMPI(MPI_Comm_size(comm, &size));
2368: PetscCall(DMPforestGetPlex(fine, &plexF));
2369: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2370: PetscCall(DMPforestGetPlex(coarse, &plexC));
2371: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2372: { /* check if the results have been cached */
2373: DM adaptCoarse, adaptFine;
2375: PetscCall(DMForestGetAdaptivityForest(coarse, &adaptCoarse));
2376: PetscCall(DMForestGetAdaptivityForest(fine, &adaptFine));
2377: if (adaptCoarse && adaptCoarse->data == fine->data) { /* coarse is adapted from fine */
2378: if (pforestC->pointSelfToAdaptSF) {
2379: PetscCall(PetscObjectReference((PetscObject)(pforestC->pointSelfToAdaptSF)));
2380: *sf = pforestC->pointSelfToAdaptSF;
2381: if (childIds) {
2382: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2383: PetscCall(PetscArraycpy(cids, pforestC->pointSelfToAdaptCids, pEndF - pStartF));
2384: *childIds = cids;
2385: }
2386: PetscFunctionReturn(PETSC_SUCCESS);
2387: } else {
2388: saveInCoarse = PETSC_TRUE;
2389: formCids = PETSC_TRUE;
2390: }
2391: } else if (adaptFine && adaptFine->data == coarse->data) { /* fine is adapted from coarse */
2392: if (pforestF->pointAdaptToSelfSF) {
2393: PetscCall(PetscObjectReference((PetscObject)(pforestF->pointAdaptToSelfSF)));
2394: *sf = pforestF->pointAdaptToSelfSF;
2395: if (childIds) {
2396: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2397: PetscCall(PetscArraycpy(cids, pforestF->pointAdaptToSelfCids, pEndF - pStartF));
2398: *childIds = cids;
2399: }
2400: PetscFunctionReturn(PETSC_SUCCESS);
2401: } else {
2402: saveInFine = PETSC_TRUE;
2403: formCids = PETSC_TRUE;
2404: }
2405: }
2406: }
2408: /* count the number of closure points that have dofs and create a list */
2409: numClosureIndices = P4EST_INSUL;
2410: /* create the datatype */
2411: PetscCallMPI(MPI_Type_contiguous(2, MPIU_INT, &nodeType));
2412: PetscCallMPI(MPI_Type_commit(&nodeType));
2413: PetscCallMPI(MPI_Op_create(DMPforestMaxSFNode, PETSC_FALSE, &sfNodeReduce));
2414: PetscCallMPI(MPI_Type_contiguous(numClosureIndices * 2, MPIU_INT, &nodeClosureType));
2415: PetscCallMPI(MPI_Type_commit(&nodeClosureType));
2416: /* everything has to go through cells: for each cell, create a list of the sfnodes in its closure */
2417: /* get lists of closure point SF nodes for every cell */
2418: PetscCall(DMPforestGetCellSFNodes(coarse, numClosureIndices, &numClosurePointsC, &closurePointsC, PETSC_TRUE));
2419: PetscCall(DMPforestGetCellSFNodes(fine, numClosureIndices, &numClosurePointsF, &closurePointsF, PETSC_FALSE));
2420: /* create pointers for tree lists */
2421: fltF = p4estF->first_local_tree;
2422: lltF = p4estF->last_local_tree;
2423: PetscCall(PetscCalloc2(lltF + 1 - fltF, &treeQuads, lltF + 1 - fltF, &treeQuadCounts));
2424: /* if the partitions don't match, ship the coarse to cover the fine */
2425: if (size > 1) {
2426: PetscInt p;
2428: for (p = 0; p < size; p++) {
2429: int equal;
2431: PetscCallP4estReturn(equal, p4est_quadrant_is_equal_piggy, (&p4estC->global_first_position[p], &p4estF->global_first_position[p]));
2432: if (!equal) break;
2433: }
2434: if (p < size) { /* non-matching distribution: send the coarse to cover the fine */
2435: PetscInt cStartC, cEndC;
2436: PetscSF coveringSF;
2437: PetscInt nleaves;
2438: PetscInt count;
2439: PetscSFNode *newClosurePointsC;
2440: p4est_quadrant_t *coverQuadsSend;
2441: p4est_topidx_t fltC = p4estC->first_local_tree;
2442: p4est_topidx_t lltC = p4estC->last_local_tree;
2443: p4est_topidx_t t;
2444: PetscMPIInt blockSizes[4] = {P4EST_DIM, 2, 1, 1};
2445: MPI_Aint blockOffsets[4] = {offsetof(p4est_quadrant_t, x), offsetof(p4est_quadrant_t, level), offsetof(p4est_quadrant_t, pad16), offsetof(p4est_quadrant_t, p)};
2446: MPI_Datatype blockTypes[4] = {MPI_INT32_T, MPI_INT8_T, MPI_INT16_T, MPI_INT32_T /* p.which_tree */};
2447: MPI_Datatype quadStruct, quadType;
2449: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, 0, &cStartC, &cEndC));
2450: PetscCall(DMPforestGetCellCoveringSF(comm, p4estC, p4estF, pforestC->cLocalStart, pforestC->cLocalEnd, &coveringSF));
2451: PetscCall(PetscSFGetGraph(coveringSF, NULL, &nleaves, NULL, NULL));
2452: PetscCall(PetscMalloc1(numClosureIndices * nleaves, &newClosurePointsC));
2453: PetscCall(PetscMalloc1(nleaves, &coverQuads));
2454: PetscCall(PetscMalloc1(cEndC - cStartC, &coverQuadsSend));
2455: count = 0;
2456: for (t = fltC; t <= lltC; t++) { /* unfortunately, we need to pack a send array, since quads are not stored packed in p4est */
2457: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2458: PetscInt q;
2460: PetscCall(PetscMemcpy(&coverQuadsSend[count], tree->quadrants.array, tree->quadrants.elem_count * sizeof(p4est_quadrant_t)));
2461: for (q = 0; (size_t)q < tree->quadrants.elem_count; q++) coverQuadsSend[count + q].p.which_tree = t;
2462: count += tree->quadrants.elem_count;
2463: }
2464: /* p is of a union type p4est_quadrant_data, but only the p.which_tree field is active at this time. So, we
2465: have a simple blockTypes[] to use. Note that quadStruct does not count potential padding in array of
2466: p4est_quadrant_t. We have to call MPI_Type_create_resized() to change upper-bound of quadStruct.
2467: */
2468: PetscCallMPI(MPI_Type_create_struct(4, blockSizes, blockOffsets, blockTypes, &quadStruct));
2469: PetscCallMPI(MPI_Type_create_resized(quadStruct, 0, sizeof(p4est_quadrant_t), &quadType));
2470: PetscCallMPI(MPI_Type_commit(&quadType));
2471: PetscCall(PetscSFBcastBegin(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2472: PetscCall(PetscSFBcastBegin(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2473: PetscCall(PetscSFBcastEnd(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2474: PetscCall(PetscSFBcastEnd(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2475: PetscCallMPI(MPI_Type_free(&quadStruct));
2476: PetscCallMPI(MPI_Type_free(&quadType));
2477: PetscCall(PetscFree(coverQuadsSend));
2478: PetscCall(PetscFree(closurePointsC));
2479: PetscCall(PetscSFDestroy(&coveringSF));
2480: closurePointsC = newClosurePointsC;
2482: /* assign tree quads based on locations in coverQuads */
2483: {
2484: PetscInt q;
2485: for (q = 0; q < nleaves; q++) {
2486: p4est_locidx_t t = coverQuads[q].p.which_tree;
2487: if (!treeQuadCounts[t - fltF]++) treeQuads[t - fltF] = &coverQuads[q];
2488: }
2489: }
2490: }
2491: }
2492: if (!coverQuads) { /* matching partitions: assign tree quads based on locations in p4est native arrays */
2493: for (t = fltF; t <= lltF; t++) {
2494: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2496: treeQuadCounts[t - fltF] = tree->quadrants.elem_count;
2497: treeQuads[t - fltF] = (p4est_quadrant_t *)tree->quadrants.array;
2498: }
2499: }
2501: {
2502: PetscInt p;
2503: PetscInt cLocalStartF;
2504: PetscSF pointSF;
2505: PetscSFNode *roots;
2506: PetscInt *rootType;
2507: DM refTree = NULL;
2508: DMLabel canonical;
2509: PetscInt *childClosures[P4EST_CHILDREN] = {NULL};
2510: PetscInt *rootClosure = NULL;
2511: PetscInt coarseOffset;
2512: PetscInt numCoarseQuads;
2514: PetscCall(PetscMalloc1(pEndF - pStartF, &roots));
2515: PetscCall(PetscMalloc1(pEndF - pStartF, &rootType));
2516: PetscCall(DMGetPointSF(fine, &pointSF));
2517: for (p = pStartF; p < pEndF; p++) {
2518: roots[p - pStartF].rank = -1;
2519: roots[p - pStartF].index = -1;
2520: rootType[p - pStartF] = -1;
2521: }
2522: if (formCids) {
2523: PetscInt child;
2525: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2526: for (p = pStartF; p < pEndF; p++) cids[p - pStartF] = -2;
2527: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2528: PetscCall(DMPlexGetTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2529: for (child = 0; child < P4EST_CHILDREN; child++) { /* get the closures of the child cells in the reference tree */
2530: PetscCall(DMPlexGetTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2531: }
2532: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2533: }
2534: cLocalStartF = pforestF->cLocalStart;
2535: for (t = fltF, coarseOffset = 0, numCoarseQuads = 0; t <= lltF; t++, coarseOffset += numCoarseQuads) {
2536: p4est_tree_t *tree = &(((p4est_tree_t *)p4estF->trees->array)[t]);
2537: PetscInt numFineQuads = tree->quadrants.elem_count;
2538: p4est_quadrant_t *coarseQuads = treeQuads[t - fltF];
2539: p4est_quadrant_t *fineQuads = (p4est_quadrant_t *)tree->quadrants.array;
2540: PetscInt i, coarseCount = 0;
2541: PetscInt offset = tree->quadrants_offset;
2542: sc_array_t coarseQuadsArray;
2544: numCoarseQuads = treeQuadCounts[t - fltF];
2545: PetscCallP4est(sc_array_init_data, (&coarseQuadsArray, coarseQuads, sizeof(p4est_quadrant_t), (size_t)numCoarseQuads));
2546: for (i = 0; i < numFineQuads; i++) {
2547: PetscInt c = i + offset;
2548: p4est_quadrant_t *quad = &fineQuads[i];
2549: p4est_quadrant_t *quadCoarse = NULL;
2550: ssize_t disjoint = -1;
2552: while (disjoint < 0 && coarseCount < numCoarseQuads) {
2553: quadCoarse = &coarseQuads[coarseCount];
2554: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2555: if (disjoint < 0) coarseCount++;
2556: }
2557: PetscCheck(disjoint == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "did not find overlapping coarse quad");
2558: if (quadCoarse->level > quad->level || (quadCoarse->level == quad->level && !transferIdent)) { /* the "coarse" mesh is finer than the fine mesh at the point: continue */
2559: if (transferIdent) { /* find corners */
2560: PetscInt j = 0;
2562: do {
2563: if (j < P4EST_CHILDREN) {
2564: p4est_quadrant_t cornerQuad;
2565: int equal;
2567: PetscCallP4est(p4est_quadrant_corner_descendant, (quad, &cornerQuad, j, quadCoarse->level));
2568: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (&cornerQuad, quadCoarse));
2569: if (equal) {
2570: PetscInt petscJ = P4estVertToPetscVert[j];
2571: PetscInt p = closurePointsF[numClosureIndices * c + (P4EST_INSUL - P4EST_CHILDREN) + petscJ].index;
2572: PetscSFNode q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + (P4EST_INSUL - P4EST_CHILDREN) + petscJ];
2574: roots[p - pStartF] = q;
2575: rootType[p - pStartF] = PETSC_MAX_INT;
2576: cids[p - pStartF] = -1;
2577: j++;
2578: }
2579: }
2580: coarseCount++;
2581: disjoint = 1;
2582: if (coarseCount < numCoarseQuads) {
2583: quadCoarse = &coarseQuads[coarseCount];
2584: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2585: }
2586: } while (!disjoint);
2587: }
2588: continue;
2589: }
2590: if (quadCoarse->level == quad->level) { /* same quad present in coarse and fine mesh */
2591: PetscInt j;
2592: for (j = 0; j < numClosureIndices; j++) {
2593: PetscInt p = closurePointsF[numClosureIndices * c + j].index;
2595: roots[p - pStartF] = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + j];
2596: rootType[p - pStartF] = PETSC_MAX_INT; /* unconditionally accept */
2597: cids[p - pStartF] = -1;
2598: }
2599: } else {
2600: PetscInt levelDiff = quad->level - quadCoarse->level;
2601: PetscInt proposedCids[P4EST_INSUL] = {0};
2603: if (formCids) {
2604: PetscInt cl;
2605: PetscInt *pointClosure = NULL;
2606: int cid;
2608: PetscCheck(levelDiff <= 1, PETSC_COMM_SELF, PETSC_ERR_USER, "Recursive child ids not implemented");
2609: PetscCallP4estReturn(cid, p4est_quadrant_child_id, (quad));
2610: PetscCall(DMPlexGetTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2611: for (cl = 0; cl < P4EST_INSUL; cl++) {
2612: PetscInt p = pointClosure[2 * cl];
2613: PetscInt point = childClosures[cid][2 * cl];
2614: PetscInt ornt = childClosures[cid][2 * cl + 1];
2615: PetscInt newcid = -1;
2616: DMPolytopeType ct;
2618: if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2619: PetscCall(DMPlexGetCellType(refTree, point, &ct));
2620: ornt = DMPolytopeConvertNewOrientation_Internal(ct, ornt);
2621: if (!cl) {
2622: newcid = cid + 1;
2623: } else {
2624: PetscInt rcl, parent, parentOrnt = 0;
2626: PetscCall(DMPlexGetTreeParent(refTree, point, &parent, NULL));
2627: if (parent == point) {
2628: newcid = -1;
2629: } else if (!parent) { /* in the root */
2630: newcid = point;
2631: } else {
2632: DMPolytopeType rct = DM_POLYTOPE_UNKNOWN;
2634: for (rcl = 1; rcl < P4EST_INSUL; rcl++) {
2635: if (rootClosure[2 * rcl] == parent) {
2636: PetscCall(DMPlexGetCellType(refTree, parent, &rct));
2637: parentOrnt = DMPolytopeConvertNewOrientation_Internal(rct, rootClosure[2 * rcl + 1]);
2638: break;
2639: }
2640: }
2641: PetscCheck(rcl < P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Couldn't find parent in root closure");
2642: PetscCall(DMPlexReferenceTreeGetChildSymmetry(refTree, parent, parentOrnt, ornt, point, DMPolytopeConvertNewOrientation_Internal(rct, pointClosure[2 * rcl + 1]), NULL, &newcid));
2643: }
2644: }
2645: if (newcid >= 0) {
2646: if (canonical) PetscCall(DMLabelGetValue(canonical, newcid, &newcid));
2647: proposedCids[cl] = newcid;
2648: }
2649: }
2650: PetscCall(DMPlexRestoreTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2651: }
2652: p4est_qcoord_t coarseBound[2][P4EST_DIM] = {
2653: {quadCoarse->x, quadCoarse->y,
2654: #if defined(P4_TO_P8)
2655: quadCoarse->z
2656: #endif
2657: },
2658: {0}
2659: };
2660: p4est_qcoord_t fineBound[2][P4EST_DIM] = {
2661: {quad->x, quad->y,
2662: #if defined(P4_TO_P8)
2663: quad->z
2664: #endif
2665: },
2666: {0}
2667: };
2668: PetscInt j;
2669: for (j = 0; j < P4EST_DIM; j++) { /* get the coordinates of cell boundaries in each direction */
2670: coarseBound[1][j] = coarseBound[0][j] + P4EST_QUADRANT_LEN(quadCoarse->level);
2671: fineBound[1][j] = fineBound[0][j] + P4EST_QUADRANT_LEN(quad->level);
2672: }
2673: for (j = 0; j < numClosureIndices; j++) {
2674: PetscInt l, p;
2675: PetscSFNode q;
2677: p = closurePointsF[numClosureIndices * c + j].index;
2678: if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2679: if (j == 0) { /* volume: ancestor is volume */
2680: l = 0;
2681: } else if (j < 1 + P4EST_FACES) { /* facet */
2682: PetscInt face = PetscFaceToP4estFace[j - 1];
2683: PetscInt direction = face / 2;
2684: PetscInt coarseFace = -1;
2686: if (coarseBound[face % 2][direction] == fineBound[face % 2][direction]) {
2687: coarseFace = face;
2688: l = 1 + P4estFaceToPetscFace[coarseFace];
2689: } else {
2690: l = 0;
2691: }
2692: #if defined(P4_TO_P8)
2693: } else if (j < 1 + P4EST_FACES + P8EST_EDGES) {
2694: PetscInt edge = PetscEdgeToP4estEdge[j - (1 + P4EST_FACES)];
2695: PetscInt direction = edge / 4;
2696: PetscInt mod = edge % 4;
2697: PetscInt coarseEdge = -1, coarseFace = -1;
2698: PetscInt minDir = PetscMin((direction + 1) % 3, (direction + 2) % 3);
2699: PetscInt maxDir = PetscMax((direction + 1) % 3, (direction + 2) % 3);
2700: PetscBool dirTest[2];
2702: dirTest[0] = (PetscBool)(coarseBound[mod % 2][minDir] == fineBound[mod % 2][minDir]);
2703: dirTest[1] = (PetscBool)(coarseBound[mod / 2][maxDir] == fineBound[mod / 2][maxDir]);
2705: if (dirTest[0] && dirTest[1]) { /* fine edge falls on coarse edge */
2706: coarseEdge = edge;
2707: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2708: } else if (dirTest[0]) { /* fine edge falls on a coarse face in the minDir direction */
2709: coarseFace = 2 * minDir + (mod % 2);
2710: l = 1 + P4estFaceToPetscFace[coarseFace];
2711: } else if (dirTest[1]) { /* fine edge falls on a coarse face in the maxDir direction */
2712: coarseFace = 2 * maxDir + (mod / 2);
2713: l = 1 + P4estFaceToPetscFace[coarseFace];
2714: } else {
2715: l = 0;
2716: }
2717: #endif
2718: } else {
2719: PetscInt vertex = PetscVertToP4estVert[P4EST_CHILDREN - (P4EST_INSUL - j)];
2720: PetscBool dirTest[P4EST_DIM];
2721: PetscInt m;
2722: PetscInt numMatch = 0;
2723: PetscInt coarseVertex = -1, coarseFace = -1;
2724: #if defined(P4_TO_P8)
2725: PetscInt coarseEdge = -1;
2726: #endif
2728: for (m = 0; m < P4EST_DIM; m++) {
2729: dirTest[m] = (PetscBool)(coarseBound[(vertex >> m) & 1][m] == fineBound[(vertex >> m) & 1][m]);
2730: if (dirTest[m]) numMatch++;
2731: }
2732: if (numMatch == P4EST_DIM) { /* vertex on vertex */
2733: coarseVertex = vertex;
2734: l = P4EST_INSUL - (P4EST_CHILDREN - P4estVertToPetscVert[coarseVertex]);
2735: } else if (numMatch == 1) { /* vertex on face */
2736: for (m = 0; m < P4EST_DIM; m++) {
2737: if (dirTest[m]) {
2738: coarseFace = 2 * m + ((vertex >> m) & 1);
2739: break;
2740: }
2741: }
2742: l = 1 + P4estFaceToPetscFace[coarseFace];
2743: #if defined(P4_TO_P8)
2744: } else if (numMatch == 2) { /* vertex on edge */
2745: for (m = 0; m < P4EST_DIM; m++) {
2746: if (!dirTest[m]) {
2747: PetscInt otherDir1 = (m + 1) % 3;
2748: PetscInt otherDir2 = (m + 2) % 3;
2749: PetscInt minDir = PetscMin(otherDir1, otherDir2);
2750: PetscInt maxDir = PetscMax(otherDir1, otherDir2);
2752: coarseEdge = m * 4 + 2 * ((vertex >> maxDir) & 1) + ((vertex >> minDir) & 1);
2753: break;
2754: }
2755: }
2756: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2757: #endif
2758: } else { /* volume */
2759: l = 0;
2760: }
2761: }
2762: q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + l];
2763: if (l > rootType[p - pStartF]) {
2764: if (l >= P4EST_INSUL - P4EST_CHILDREN) { /* vertex on vertex: unconditional acceptance */
2765: if (transferIdent) {
2766: roots[p - pStartF] = q;
2767: rootType[p - pStartF] = PETSC_MAX_INT;
2768: if (formCids) cids[p - pStartF] = -1;
2769: }
2770: } else {
2771: PetscInt k, thisp = p, limit;
2773: roots[p - pStartF] = q;
2774: rootType[p - pStartF] = l;
2775: if (formCids) cids[p - pStartF] = proposedCids[j];
2776: limit = transferIdent ? levelDiff : (levelDiff - 1);
2777: for (k = 0; k < limit; k++) {
2778: PetscInt parent;
2780: PetscCall(DMPlexGetTreeParent(plexF, thisp, &parent, NULL));
2781: if (parent == thisp) break;
2783: roots[parent - pStartF] = q;
2784: rootType[parent - pStartF] = PETSC_MAX_INT;
2785: if (formCids) cids[parent - pStartF] = -1;
2786: thisp = parent;
2787: }
2788: }
2789: }
2790: }
2791: }
2792: }
2793: }
2795: /* 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 */
2796: if (size > 1) {
2797: PetscInt *rootTypeCopy, p;
2799: PetscCall(PetscMalloc1(pEndF - pStartF, &rootTypeCopy));
2800: PetscCall(PetscArraycpy(rootTypeCopy, rootType, pEndF - pStartF));
2801: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2802: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2803: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2804: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2805: for (p = pStartF; p < pEndF; p++) {
2806: 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 */
2807: roots[p - pStartF].rank = -1;
2808: roots[p - pStartF].index = -1;
2809: }
2810: if (formCids && rootTypeCopy[p - pStartF] == PETSC_MAX_INT) { cids[p - pStartF] = -1; /* we have found an antecedent that is the same: no child id */ }
2811: }
2812: PetscCall(PetscFree(rootTypeCopy));
2813: PetscCall(PetscSFReduceBegin(pointSF, nodeType, roots, roots, sfNodeReduce));
2814: PetscCall(PetscSFReduceEnd(pointSF, nodeType, roots, roots, sfNodeReduce));
2815: PetscCall(PetscSFBcastBegin(pointSF, nodeType, roots, roots, MPI_REPLACE));
2816: PetscCall(PetscSFBcastEnd(pointSF, nodeType, roots, roots, MPI_REPLACE));
2817: }
2818: PetscCall(PetscFree(rootType));
2820: {
2821: PetscInt numRoots;
2822: PetscInt numLeaves;
2823: PetscInt *leaves;
2824: PetscSFNode *iremote;
2825: /* count leaves */
2827: numRoots = pEndC - pStartC;
2829: numLeaves = 0;
2830: for (p = pStartF; p < pEndF; p++) {
2831: if (roots[p - pStartF].index >= 0) numLeaves++;
2832: }
2833: PetscCall(PetscMalloc1(numLeaves, &leaves));
2834: PetscCall(PetscMalloc1(numLeaves, &iremote));
2835: numLeaves = 0;
2836: for (p = pStartF; p < pEndF; p++) {
2837: if (roots[p - pStartF].index >= 0) {
2838: leaves[numLeaves] = p - pStartF;
2839: iremote[numLeaves] = roots[p - pStartF];
2840: numLeaves++;
2841: }
2842: }
2843: PetscCall(PetscFree(roots));
2844: PetscCall(PetscSFCreate(comm, sf));
2845: if (numLeaves == (pEndF - pStartF)) {
2846: PetscCall(PetscFree(leaves));
2847: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2848: } else {
2849: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, leaves, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2850: }
2851: }
2852: if (formCids) {
2853: PetscSF pointSF;
2854: PetscInt child;
2856: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2857: PetscCall(DMGetPointSF(plexF, &pointSF));
2858: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2859: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2860: if (childIds) *childIds = cids;
2861: for (child = 0; child < P4EST_CHILDREN; child++) PetscCall(DMPlexRestoreTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2862: PetscCall(DMPlexRestoreTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2863: }
2864: }
2865: if (saveInCoarse) { /* cache results */
2866: PetscCall(PetscObjectReference((PetscObject)*sf));
2867: pforestC->pointSelfToAdaptSF = *sf;
2868: if (!childIds) {
2869: pforestC->pointSelfToAdaptCids = cids;
2870: } else {
2871: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestC->pointSelfToAdaptCids));
2872: PetscCall(PetscArraycpy(pforestC->pointSelfToAdaptCids, cids, pEndF - pStartF));
2873: }
2874: } else if (saveInFine) {
2875: PetscCall(PetscObjectReference((PetscObject)*sf));
2876: pforestF->pointAdaptToSelfSF = *sf;
2877: if (!childIds) {
2878: pforestF->pointAdaptToSelfCids = cids;
2879: } else {
2880: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestF->pointAdaptToSelfCids));
2881: PetscCall(PetscArraycpy(pforestF->pointAdaptToSelfCids, cids, pEndF - pStartF));
2882: }
2883: }
2884: PetscCall(PetscFree2(treeQuads, treeQuadCounts));
2885: PetscCall(PetscFree(coverQuads));
2886: PetscCall(PetscFree(closurePointsC));
2887: PetscCall(PetscFree(closurePointsF));
2888: PetscCallMPI(MPI_Type_free(&nodeClosureType));
2889: PetscCallMPI(MPI_Op_free(&sfNodeReduce));
2890: PetscCallMPI(MPI_Type_free(&nodeType));
2891: PetscFunctionReturn(PETSC_SUCCESS);
2892: }
2894: /* children are sf leaves of parents */
2895: static PetscErrorCode DMPforestGetTransferSF_Internal(DM coarse, DM fine, const PetscInt dofPerDim[], PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2896: {
2897: MPI_Comm comm;
2898: PetscMPIInt rank;
2899: DM_Forest_pforest *pforestC, *pforestF;
2900: DM plexC, plexF;
2901: PetscInt pStartC, pEndC, pStartF, pEndF;
2902: PetscSF pointTransferSF;
2903: PetscBool allOnes = PETSC_TRUE;
2905: PetscFunctionBegin;
2906: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2907: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2908: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2909: comm = PetscObjectComm((PetscObject)coarse);
2910: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2912: {
2913: PetscInt i;
2914: for (i = 0; i <= P4EST_DIM; i++) {
2915: if (dofPerDim[i] != 1) {
2916: allOnes = PETSC_FALSE;
2917: break;
2918: }
2919: }
2920: }
2921: PetscCall(DMPforestGetTransferSF_Point(coarse, fine, &pointTransferSF, transferIdent, childIds));
2922: if (allOnes) {
2923: *sf = pointTransferSF;
2924: PetscFunctionReturn(PETSC_SUCCESS);
2925: }
2927: PetscCall(DMPforestGetPlex(fine, &plexF));
2928: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2929: PetscCall(DMPforestGetPlex(coarse, &plexC));
2930: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2931: {
2932: PetscInt numRoots;
2933: PetscInt numLeaves;
2934: const PetscInt *leaves;
2935: const PetscSFNode *iremote;
2936: PetscInt d;
2937: PetscSection leafSection, rootSection;
2938: /* count leaves */
2940: PetscCall(PetscSFGetGraph(pointTransferSF, &numRoots, &numLeaves, &leaves, &iremote));
2941: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &rootSection));
2942: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &leafSection));
2943: PetscCall(PetscSectionSetChart(rootSection, pStartC, pEndC));
2944: PetscCall(PetscSectionSetChart(leafSection, pStartF, pEndF));
2946: for (d = 0; d <= P4EST_DIM; d++) {
2947: PetscInt startC, endC, e;
2949: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, P4EST_DIM - d, &startC, &endC));
2950: for (e = startC; e < endC; e++) PetscCall(PetscSectionSetDof(rootSection, e, dofPerDim[d]));
2951: }
2953: for (d = 0; d <= P4EST_DIM; d++) {
2954: PetscInt startF, endF, e;
2956: PetscCall(DMPlexGetSimplexOrBoxCells(plexF, P4EST_DIM - d, &startF, &endF));
2957: for (e = startF; e < endF; e++) PetscCall(PetscSectionSetDof(leafSection, e, dofPerDim[d]));
2958: }
2960: PetscCall(PetscSectionSetUp(rootSection));
2961: PetscCall(PetscSectionSetUp(leafSection));
2962: {
2963: PetscInt nroots, nleaves;
2964: PetscInt *mine, i, p;
2965: PetscInt *offsets, *offsetsRoot;
2966: PetscSFNode *remote;
2968: PetscCall(PetscMalloc1(pEndF - pStartF, &offsets));
2969: PetscCall(PetscMalloc1(pEndC - pStartC, &offsetsRoot));
2970: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionGetOffset(rootSection, p, &offsetsRoot[p - pStartC]));
2971: PetscCall(PetscSFBcastBegin(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2972: PetscCall(PetscSFBcastEnd(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2973: PetscCall(PetscSectionGetStorageSize(rootSection, &nroots));
2974: nleaves = 0;
2975: for (i = 0; i < numLeaves; i++) {
2976: PetscInt leaf = leaves ? leaves[i] : i;
2977: PetscInt dof;
2979: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2980: nleaves += dof;
2981: }
2982: PetscCall(PetscMalloc1(nleaves, &mine));
2983: PetscCall(PetscMalloc1(nleaves, &remote));
2984: nleaves = 0;
2985: for (i = 0; i < numLeaves; i++) {
2986: PetscInt leaf = leaves ? leaves[i] : i;
2987: PetscInt dof;
2988: PetscInt off, j;
2990: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2991: PetscCall(PetscSectionGetOffset(leafSection, leaf, &off));
2992: for (j = 0; j < dof; j++) {
2993: remote[nleaves].rank = iremote[i].rank;
2994: remote[nleaves].index = offsets[leaf] + j;
2995: mine[nleaves++] = off + j;
2996: }
2997: }
2998: PetscCall(PetscFree(offsetsRoot));
2999: PetscCall(PetscFree(offsets));
3000: PetscCall(PetscSFCreate(comm, sf));
3001: PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
3002: }
3003: PetscCall(PetscSectionDestroy(&leafSection));
3004: PetscCall(PetscSectionDestroy(&rootSection));
3005: PetscCall(PetscSFDestroy(&pointTransferSF));
3006: }
3007: PetscFunctionReturn(PETSC_SUCCESS);
3008: }
3010: static PetscErrorCode DMPforestGetTransferSF(DM dmA, DM dmB, const PetscInt dofPerDim[], PetscSF *sfAtoB, PetscSF *sfBtoA)
3011: {
3012: DM adaptA, adaptB;
3013: DMAdaptFlag purpose;
3015: PetscFunctionBegin;
3016: PetscCall(DMForestGetAdaptivityForest(dmA, &adaptA));
3017: PetscCall(DMForestGetAdaptivityForest(dmB, &adaptB));
3018: /* it is more efficient when the coarser mesh is the first argument: reorder if we know one is coarser than the other */
3019: if (adaptA && adaptA->data == dmB->data) { /* dmA was adapted from dmB */
3020: PetscCall(DMForestGetAdaptivityPurpose(dmA, &purpose));
3021: if (purpose == DM_ADAPT_REFINE) {
3022: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3023: PetscFunctionReturn(PETSC_SUCCESS);
3024: }
3025: } else if (adaptB && adaptB->data == dmA->data) { /* dmB was adapted from dmA */
3026: PetscCall(DMForestGetAdaptivityPurpose(dmB, &purpose));
3027: if (purpose == DM_ADAPT_COARSEN) {
3028: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3029: PetscFunctionReturn(PETSC_SUCCESS);
3030: }
3031: }
3032: if (sfAtoB) PetscCall(DMPforestGetTransferSF_Internal(dmA, dmB, dofPerDim, sfAtoB, PETSC_TRUE, NULL));
3033: if (sfBtoA) PetscCall(DMPforestGetTransferSF_Internal(dmB, dmA, dofPerDim, sfBtoA, (PetscBool)(sfAtoB == NULL), NULL));
3034: PetscFunctionReturn(PETSC_SUCCESS);
3035: }
3037: static PetscErrorCode DMPforestLabelsInitialize(DM dm, DM plex)
3038: {
3039: DM_Forest *forest = (DM_Forest *)dm->data;
3040: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
3041: PetscInt cLocalStart, cLocalEnd, cStart, cEnd, fStart, fEnd, eStart, eEnd, vStart, vEnd;
3042: PetscInt cStartBase, cEndBase, fStartBase, fEndBase, vStartBase, vEndBase, eStartBase, eEndBase;
3043: PetscInt pStart, pEnd, pStartBase, pEndBase, p;
3044: DM base;
3045: PetscInt *star = NULL, starSize;
3046: DMLabelLink next = dm->labels;
3047: PetscInt guess = 0;
3048: p4est_topidx_t num_trees = pforest->topo->conn->num_trees;
3050: PetscFunctionBegin;
3051: pforest->labelsFinalized = PETSC_TRUE;
3052: cLocalStart = pforest->cLocalStart;
3053: cLocalEnd = pforest->cLocalEnd;
3054: PetscCall(DMForestGetBaseDM(dm, &base));
3055: if (!base) {
3056: if (pforest->ghostName) { /* insert a label to make the boundaries, with stratum values denoting which face of the element touches the boundary */
3057: p4est_connectivity_t *conn = pforest->topo->conn;
3058: p4est_t *p4est = pforest->forest;
3059: p4est_tree_t *trees = (p4est_tree_t *)p4est->trees->array;
3060: p4est_topidx_t t, flt = p4est->first_local_tree;
3061: p4est_topidx_t llt = pforest->forest->last_local_tree;
3062: DMLabel ghostLabel;
3063: PetscInt c;
3065: PetscCall(DMCreateLabel(plex, pforest->ghostName));
3066: PetscCall(DMGetLabel(plex, pforest->ghostName, &ghostLabel));
3067: for (c = cLocalStart, t = flt; t <= llt; t++) {
3068: p4est_tree_t *tree = &trees[t];
3069: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3070: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3071: PetscInt q;
3073: for (q = 0; q < numQuads; q++, c++) {
3074: p4est_quadrant_t *quad = &quads[q];
3075: PetscInt f;
3077: for (f = 0; f < P4EST_FACES; f++) {
3078: p4est_quadrant_t neigh;
3079: int isOutside;
3081: PetscCallP4est(p4est_quadrant_face_neighbor, (quad, f, &neigh));
3082: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&neigh));
3083: if (isOutside) {
3084: p4est_topidx_t nt;
3085: PetscInt nf;
3087: nt = conn->tree_to_tree[t * P4EST_FACES + f];
3088: nf = (PetscInt)conn->tree_to_face[t * P4EST_FACES + f];
3089: nf = nf % P4EST_FACES;
3090: if (nt == t && nf == f) {
3091: PetscInt plexF = P4estFaceToPetscFace[f];
3092: const PetscInt *cone;
3094: PetscCall(DMPlexGetCone(plex, c, &cone));
3095: PetscCall(DMLabelSetValue(ghostLabel, cone[plexF], plexF + 1));
3096: }
3097: }
3098: }
3099: }
3100: }
3101: }
3102: PetscFunctionReturn(PETSC_SUCCESS);
3103: }
3104: PetscCall(DMPlexGetSimplexOrBoxCells(base, 0, &cStartBase, &cEndBase));
3105: PetscCall(DMPlexGetSimplexOrBoxCells(base, 1, &fStartBase, &fEndBase));
3106: PetscCall(DMPlexGetSimplexOrBoxCells(base, P4EST_DIM - 1, &eStartBase, &eEndBase));
3107: PetscCall(DMPlexGetDepthStratum(base, 0, &vStartBase, &vEndBase));
3109: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
3110: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 1, &fStart, &fEnd));
3111: PetscCall(DMPlexGetSimplexOrBoxCells(plex, P4EST_DIM - 1, &eStart, &eEnd));
3112: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3114: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3115: PetscCall(DMPlexGetChart(base, &pStartBase, &pEndBase));
3116: /* go through the mesh: use star to find a quadrant that borders a point. Use the closure to determine the
3117: * orientation of the quadrant relative to that point. Use that to relate the point to the numbering in the base
3118: * mesh, and extract a label value (since the base mesh is redundantly distributed, can be found locally). */
3119: while (next) {
3120: DMLabel baseLabel;
3121: DMLabel label = next->label;
3122: PetscBool isDepth, isCellType, isGhost, isVTK, isSpmap;
3123: const char *name;
3125: PetscCall(PetscObjectGetName((PetscObject)label, &name));
3126: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3127: if (isDepth) {
3128: next = next->next;
3129: continue;
3130: }
3131: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3132: if (isCellType) {
3133: next = next->next;
3134: continue;
3135: }
3136: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3137: if (isGhost) {
3138: next = next->next;
3139: continue;
3140: }
3141: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3142: if (isVTK) {
3143: next = next->next;
3144: continue;
3145: }
3146: PetscCall(PetscStrcmp(name, "_forest_base_subpoint_map", &isSpmap));
3147: if (!isSpmap) {
3148: PetscCall(DMGetLabel(base, name, &baseLabel));
3149: if (!baseLabel) {
3150: next = next->next;
3151: continue;
3152: }
3153: PetscCall(DMLabelCreateIndex(baseLabel, pStartBase, pEndBase));
3154: } else baseLabel = NULL;
3156: for (p = pStart; p < pEnd; p++) {
3157: PetscInt s, c = -1, l;
3158: PetscInt *closure = NULL, closureSize;
3159: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3160: p4est_tree_t *trees = (p4est_tree_t *)pforest->forest->trees->array;
3161: p4est_quadrant_t *q;
3162: PetscInt t, val;
3163: PetscBool zerosupportpoint = PETSC_FALSE;
3165: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3166: for (s = 0; s < starSize; s++) {
3167: PetscInt point = star[2 * s];
3169: if (cStart <= point && point < cEnd) {
3170: PetscCall(DMPlexGetTransitiveClosure(plex, point, PETSC_TRUE, &closureSize, &closure));
3171: for (l = 0; l < closureSize; l++) {
3172: PetscInt qParent = closure[2 * l], q, pp = p, pParent = p;
3173: do { /* check parents of q */
3174: q = qParent;
3175: if (q == p) {
3176: c = point;
3177: break;
3178: }
3179: PetscCall(DMPlexGetTreeParent(plex, q, &qParent, NULL));
3180: } while (qParent != q);
3181: if (c != -1) break;
3182: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3183: q = closure[2 * l];
3184: while (pParent != pp) { /* check parents of p */
3185: pp = pParent;
3186: if (pp == q) {
3187: c = point;
3188: break;
3189: }
3190: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3191: }
3192: if (c != -1) break;
3193: }
3194: PetscCall(DMPlexRestoreTransitiveClosure(plex, point, PETSC_TRUE, NULL, &closure));
3195: if (l < closureSize) break;
3196: } else {
3197: PetscInt supportSize;
3199: PetscCall(DMPlexGetSupportSize(plex, point, &supportSize));
3200: zerosupportpoint = (PetscBool)(zerosupportpoint || !supportSize);
3201: }
3202: }
3203: if (c < 0) {
3204: const char *prefix;
3205: PetscBool print = PETSC_FALSE;
3207: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
3208: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_forest_print_label_error", &print, NULL));
3209: if (print) {
3210: PetscInt i;
3212: PetscCall(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));
3213: for (i = 0; i < starSize; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, " star[%" PetscInt_FMT "] = %" PetscInt_FMT ",%" PetscInt_FMT "\n", i, star[2 * i], star[2 * i + 1]));
3214: }
3215: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3216: if (zerosupportpoint) continue;
3217: else
3218: 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");
3219: }
3220: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3222: if (c < cLocalStart) {
3223: /* get from the beginning of the ghost layer */
3224: q = &(ghosts[c]);
3225: t = (PetscInt)q->p.which_tree;
3226: } else if (c < cLocalEnd) {
3227: PetscInt lo = 0, hi = num_trees;
3228: /* get from local quadrants: have to find the right tree */
3230: c -= cLocalStart;
3232: do {
3233: p4est_tree_t *tree;
3235: PetscCheck(guess >= lo && guess < num_trees && lo < hi, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed binary search");
3236: tree = &trees[guess];
3237: if (c < tree->quadrants_offset) {
3238: hi = guess;
3239: } else if (c < tree->quadrants_offset + (PetscInt)tree->quadrants.elem_count) {
3240: q = &((p4est_quadrant_t *)tree->quadrants.array)[c - (PetscInt)tree->quadrants_offset];
3241: t = guess;
3242: break;
3243: } else {
3244: lo = guess + 1;
3245: }
3246: guess = lo + (hi - lo) / 2;
3247: } while (1);
3248: } else {
3249: /* get from the end of the ghost layer */
3250: c -= (cLocalEnd - cLocalStart);
3252: q = &(ghosts[c]);
3253: t = (PetscInt)q->p.which_tree;
3254: }
3256: if (l == 0) { /* cell */
3257: if (baseLabel) {
3258: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3259: } else {
3260: val = t + cStartBase;
3261: }
3262: PetscCall(DMLabelSetValue(label, p, val));
3263: } else if (l >= 1 && l < 1 + P4EST_FACES) { /* facet */
3264: p4est_quadrant_t nq;
3265: int isInside;
3267: l = PetscFaceToP4estFace[l - 1];
3268: PetscCallP4est(p4est_quadrant_face_neighbor, (q, l, &nq));
3269: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3270: if (isInside) {
3271: /* this facet is in the interior of a tree, so it inherits the label of the tree */
3272: if (baseLabel) {
3273: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3274: } else {
3275: val = t + cStartBase;
3276: }
3277: PetscCall(DMLabelSetValue(label, p, val));
3278: } else {
3279: PetscInt f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + l];
3281: if (baseLabel) {
3282: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3283: } else {
3284: val = f + fStartBase;
3285: }
3286: PetscCall(DMLabelSetValue(label, p, val));
3287: }
3288: #if defined(P4_TO_P8)
3289: } else if (l >= 1 + P4EST_FACES && l < 1 + P4EST_FACES + P8EST_EDGES) { /* edge */
3290: p4est_quadrant_t nq;
3291: int isInside;
3293: l = PetscEdgeToP4estEdge[l - (1 + P4EST_FACES)];
3294: PetscCallP4est(p8est_quadrant_edge_neighbor, (q, l, &nq));
3295: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3296: if (isInside) {
3297: /* this edge is in the interior of a tree, so it inherits the label of the tree */
3298: if (baseLabel) {
3299: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3300: } else {
3301: val = t + cStartBase;
3302: }
3303: PetscCall(DMLabelSetValue(label, p, val));
3304: } else {
3305: int isOutsideFace;
3307: PetscCallP4estReturn(isOutsideFace, p4est_quadrant_is_outside_face, (&nq));
3308: if (isOutsideFace) {
3309: PetscInt f;
3311: if (nq.x < 0) {
3312: f = 0;
3313: } else if (nq.x >= P4EST_ROOT_LEN) {
3314: f = 1;
3315: } else if (nq.y < 0) {
3316: f = 2;
3317: } else if (nq.y >= P4EST_ROOT_LEN) {
3318: f = 3;
3319: } else if (nq.z < 0) {
3320: f = 4;
3321: } else {
3322: f = 5;
3323: }
3324: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3325: if (baseLabel) {
3326: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3327: } else {
3328: val = f + fStartBase;
3329: }
3330: PetscCall(DMLabelSetValue(label, p, val));
3331: } else { /* the quadrant edge corresponds to the tree edge */
3332: PetscInt e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + l];
3334: if (baseLabel) {
3335: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3336: } else {
3337: val = e + eStartBase;
3338: }
3339: PetscCall(DMLabelSetValue(label, p, val));
3340: }
3341: }
3342: #endif
3343: } else { /* vertex */
3344: p4est_quadrant_t nq;
3345: int isInside;
3347: #if defined(P4_TO_P8)
3348: l = PetscVertToP4estVert[l - (1 + P4EST_FACES + P8EST_EDGES)];
3349: #else
3350: l = PetscVertToP4estVert[l - (1 + P4EST_FACES)];
3351: #endif
3352: PetscCallP4est(p4est_quadrant_corner_neighbor, (q, l, &nq));
3353: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3354: if (isInside) {
3355: if (baseLabel) {
3356: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3357: } else {
3358: val = t + cStartBase;
3359: }
3360: PetscCall(DMLabelSetValue(label, p, val));
3361: } else {
3362: int isOutside;
3364: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&nq));
3365: if (isOutside) {
3366: PetscInt f = -1;
3368: if (nq.x < 0) {
3369: f = 0;
3370: } else if (nq.x >= P4EST_ROOT_LEN) {
3371: f = 1;
3372: } else if (nq.y < 0) {
3373: f = 2;
3374: } else if (nq.y >= P4EST_ROOT_LEN) {
3375: f = 3;
3376: #if defined(P4_TO_P8)
3377: } else if (nq.z < 0) {
3378: f = 4;
3379: } else {
3380: f = 5;
3381: #endif
3382: }
3383: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3384: if (baseLabel) {
3385: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3386: } else {
3387: val = f + fStartBase;
3388: }
3389: PetscCall(DMLabelSetValue(label, p, val));
3390: continue;
3391: }
3392: #if defined(P4_TO_P8)
3393: PetscCallP4estReturn(isOutside, p8est_quadrant_is_outside_edge, (&nq));
3394: if (isOutside) {
3395: /* outside edge */
3396: PetscInt e = -1;
3398: if (nq.x >= 0 && nq.x < P4EST_ROOT_LEN) {
3399: if (nq.z < 0) {
3400: if (nq.y < 0) {
3401: e = 0;
3402: } else {
3403: e = 1;
3404: }
3405: } else {
3406: if (nq.y < 0) {
3407: e = 2;
3408: } else {
3409: e = 3;
3410: }
3411: }
3412: } else if (nq.y >= 0 && nq.y < P4EST_ROOT_LEN) {
3413: if (nq.z < 0) {
3414: if (nq.x < 0) {
3415: e = 4;
3416: } else {
3417: e = 5;
3418: }
3419: } else {
3420: if (nq.x < 0) {
3421: e = 6;
3422: } else {
3423: e = 7;
3424: }
3425: }
3426: } else {
3427: if (nq.y < 0) {
3428: if (nq.x < 0) {
3429: e = 8;
3430: } else {
3431: e = 9;
3432: }
3433: } else {
3434: if (nq.x < 0) {
3435: e = 10;
3436: } else {
3437: e = 11;
3438: }
3439: }
3440: }
3442: e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + e];
3443: if (baseLabel) {
3444: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3445: } else {
3446: val = e + eStartBase;
3447: }
3448: PetscCall(DMLabelSetValue(label, p, val));
3449: continue;
3450: }
3451: #endif
3452: {
3453: /* outside vertex: same corner as quadrant corner */
3454: PetscInt v = pforest->topo->conn->tree_to_corner[P4EST_CHILDREN * t + l];
3456: if (baseLabel) {
3457: PetscCall(DMLabelGetValue(baseLabel, v + vStartBase, &val));
3458: } else {
3459: val = v + vStartBase;
3460: }
3461: PetscCall(DMLabelSetValue(label, p, val));
3462: }
3463: }
3464: }
3465: }
3466: next = next->next;
3467: }
3468: PetscFunctionReturn(PETSC_SUCCESS);
3469: }
3471: static PetscErrorCode DMPforestLabelsFinalize(DM dm, DM plex)
3472: {
3473: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
3474: DM adapt;
3476: PetscFunctionBegin;
3477: if (pforest->labelsFinalized) PetscFunctionReturn(PETSC_SUCCESS);
3478: pforest->labelsFinalized = PETSC_TRUE;
3479: PetscCall(DMForestGetAdaptivityForest(dm, &adapt));
3480: if (!adapt) {
3481: /* Initialize labels from the base dm */
3482: PetscCall(DMPforestLabelsInitialize(dm, plex));
3483: } else {
3484: PetscInt dofPerDim[4] = {1, 1, 1, 1};
3485: PetscSF transferForward, transferBackward, pointSF;
3486: PetscInt pStart, pEnd, pStartA, pEndA;
3487: PetscInt *values, *adaptValues;
3488: DMLabelLink next = adapt->labels;
3489: DMLabel adaptLabel;
3490: DM adaptPlex;
3492: PetscCall(DMForestGetAdaptivityLabel(dm, &adaptLabel));
3493: PetscCall(DMPforestGetPlex(adapt, &adaptPlex));
3494: PetscCall(DMPforestGetTransferSF(adapt, dm, dofPerDim, &transferForward, &transferBackward));
3495: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3496: PetscCall(DMPlexGetChart(adaptPlex, &pStartA, &pEndA));
3497: PetscCall(PetscMalloc2(pEnd - pStart, &values, pEndA - pStartA, &adaptValues));
3498: PetscCall(DMGetPointSF(plex, &pointSF));
3499: if (PetscDefined(USE_DEBUG)) {
3500: PetscInt p;
3501: for (p = pStartA; p < pEndA; p++) adaptValues[p - pStartA] = -1;
3502: for (p = pStart; p < pEnd; p++) values[p - pStart] = -2;
3503: if (transferForward) {
3504: PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3505: PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3506: }
3507: if (transferBackward) {
3508: PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3509: PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3510: }
3511: for (p = pStart; p < pEnd; p++) {
3512: PetscInt q = p, parent;
3514: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3515: while (parent != q) {
3516: if (values[parent] == -2) values[parent] = values[q];
3517: q = parent;
3518: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3519: }
3520: }
3521: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3522: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3523: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3524: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3525: for (p = pStart; p < pEnd; p++) PetscCheck(values[p - pStart] != -2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "uncovered point %" PetscInt_FMT, p);
3526: }
3527: while (next) {
3528: DMLabel nextLabel = next->label;
3529: const char *name;
3530: PetscBool isDepth, isCellType, isGhost, isVTK;
3531: DMLabel label;
3532: PetscInt p;
3534: PetscCall(PetscObjectGetName((PetscObject)nextLabel, &name));
3535: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3536: if (isDepth) {
3537: next = next->next;
3538: continue;
3539: }
3540: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3541: if (isCellType) {
3542: next = next->next;
3543: continue;
3544: }
3545: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3546: if (isGhost) {
3547: next = next->next;
3548: continue;
3549: }
3550: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3551: if (isVTK) {
3552: next = next->next;
3553: continue;
3554: }
3555: if (nextLabel == adaptLabel) {
3556: next = next->next;
3557: continue;
3558: }
3559: /* label was created earlier */
3560: PetscCall(DMGetLabel(dm, name, &label));
3561: for (p = pStartA; p < pEndA; p++) PetscCall(DMLabelGetValue(nextLabel, p, &adaptValues[p]));
3562: for (p = pStart; p < pEnd; p++) values[p] = PETSC_MIN_INT;
3564: if (transferForward) PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3565: if (transferBackward) PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3566: if (transferForward) PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3567: if (transferBackward) PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3568: for (p = pStart; p < pEnd; p++) {
3569: PetscInt q = p, parent;
3571: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3572: while (parent != q) {
3573: if (values[parent] == PETSC_MIN_INT) values[parent] = values[q];
3574: q = parent;
3575: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3576: }
3577: }
3578: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3579: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3580: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3581: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3583: for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(label, p, values[p]));
3584: next = next->next;
3585: }
3586: PetscCall(PetscFree2(values, adaptValues));
3587: PetscCall(PetscSFDestroy(&transferForward));
3588: PetscCall(PetscSFDestroy(&transferBackward));
3589: pforest->labelsFinalized = PETSC_TRUE;
3590: }
3591: PetscFunctionReturn(PETSC_SUCCESS);
3592: }
3594: 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)
3595: {
3596: PetscInt closureSize, c, coordStart, coordEnd, coordDim;
3597: PetscInt *closure = NULL;
3598: PetscSection coordSec;
3600: PetscFunctionBegin;
3601: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3602: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3603: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3604: PetscCall(DMPlexGetTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3605: for (c = 0; c < closureSize; c++) {
3606: PetscInt point = closure[2 * c];
3608: if (point >= coordStart && point < coordEnd) {
3609: PetscInt dof, off;
3610: PetscInt nCoords, i;
3611: PetscCall(PetscSectionGetDof(coordSec, point, &dof));
3612: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3613: nCoords = dof / coordDim;
3614: PetscCall(PetscSectionGetOffset(coordSec, point, &off));
3615: for (i = 0; i < nCoords; i++) {
3616: PetscScalar *coord = &coords[off + i * coordDim];
3617: double coordP4est[3] = {0.};
3618: double coordP4estMapped[3] = {0.};
3619: PetscInt j;
3620: PetscReal treeCoords[P4EST_CHILDREN][3] = {{0.}};
3621: PetscReal eta[3] = {0.};
3622: PetscInt numRounds = 10;
3623: PetscReal coordGuess[3] = {0.};
3625: eta[0] = (PetscReal)q->x / (PetscReal)P4EST_ROOT_LEN;
3626: eta[1] = (PetscReal)q->y / (PetscReal)P4EST_ROOT_LEN;
3627: #if defined(P4_TO_P8)
3628: eta[2] = (PetscReal)q->z / (PetscReal)P4EST_ROOT_LEN;
3629: #endif
3631: for (j = 0; j < P4EST_CHILDREN; j++) {
3632: PetscInt k;
3634: for (k = 0; k < 3; k++) treeCoords[j][k] = conn->vertices[3 * conn->tree_to_vertex[P4EST_CHILDREN * t + j] + k];
3635: }
3637: for (j = 0; j < P4EST_CHILDREN; j++) {
3638: PetscInt k;
3639: PetscReal prod = 1.;
3641: for (k = 0; k < P4EST_DIM; k++) prod *= (j & (1 << k)) ? eta[k] : (1. - eta[k]);
3642: for (k = 0; k < 3; k++) coordGuess[k] += prod * treeCoords[j][k];
3643: }
3645: for (j = 0; j < numRounds; j++) {
3646: PetscInt dir;
3648: for (dir = 0; dir < P4EST_DIM; dir++) {
3649: PetscInt k;
3650: PetscReal diff[3];
3651: PetscReal dXdeta[3] = {0.};
3652: PetscReal rhs, scale, update;
3654: for (k = 0; k < 3; k++) diff[k] = coordP4est[k] - coordGuess[k];
3655: for (k = 0; k < P4EST_CHILDREN; k++) {
3656: PetscInt l;
3657: PetscReal prod = 1.;
3659: for (l = 0; l < P4EST_DIM; l++) {
3660: if (l == dir) {
3661: prod *= (k & (1 << l)) ? 1. : -1.;
3662: } else {
3663: prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3664: }
3665: }
3666: for (l = 0; l < 3; l++) dXdeta[l] += prod * treeCoords[k][l];
3667: }
3668: rhs = 0.;
3669: scale = 0;
3670: for (k = 0; k < 3; k++) {
3671: rhs += diff[k] * dXdeta[k];
3672: scale += dXdeta[k] * dXdeta[k];
3673: }
3674: update = rhs / scale;
3675: eta[dir] += update;
3676: eta[dir] = PetscMin(eta[dir], 1.);
3677: eta[dir] = PetscMax(eta[dir], 0.);
3679: coordGuess[0] = coordGuess[1] = coordGuess[2] = 0.;
3680: for (k = 0; k < P4EST_CHILDREN; k++) {
3681: PetscInt l;
3682: PetscReal prod = 1.;
3684: for (l = 0; l < P4EST_DIM; l++) prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3685: for (l = 0; l < 3; l++) coordGuess[l] += prod * treeCoords[k][l];
3686: }
3687: }
3688: }
3689: for (j = 0; j < 3; j++) coordP4est[j] = (double)eta[j];
3691: if (geom) {
3692: (geom->X)(geom, t, coordP4est, coordP4estMapped);
3693: for (j = 0; j < coordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3694: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded");
3695: }
3696: }
3697: }
3698: PetscCall(DMPlexRestoreTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3699: PetscFunctionReturn(PETSC_SUCCESS);
3700: }
3702: static PetscErrorCode DMPforestMapCoordinates(DM dm, DM plex)
3703: {
3704: DM_Forest *forest;
3705: DM_Forest_pforest *pforest;
3706: p4est_geometry_t *geom;
3707: PetscInt cLocalStart, cLocalEnd;
3708: Vec coordLocalVec;
3709: PetscScalar *coords;
3710: p4est_topidx_t flt, llt, t;
3711: p4est_tree_t *trees;
3712: PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
3713: void *mapCtx;
3715: PetscFunctionBegin;
3716: forest = (DM_Forest *)dm->data;
3717: pforest = (DM_Forest_pforest *)forest->data;
3718: geom = pforest->topo->geom;
3719: PetscCall(DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx));
3720: if (!geom && !map) PetscFunctionReturn(PETSC_SUCCESS);
3721: PetscCall(DMGetCoordinatesLocal(plex, &coordLocalVec));
3722: PetscCall(VecGetArray(coordLocalVec, &coords));
3723: cLocalStart = pforest->cLocalStart;
3724: cLocalEnd = pforest->cLocalEnd;
3725: flt = pforest->forest->first_local_tree;
3726: llt = pforest->forest->last_local_tree;
3727: trees = (p4est_tree_t *)pforest->forest->trees->array;
3728: if (map) { /* apply the map directly to the existing coordinates */
3729: PetscSection coordSec;
3730: PetscInt coordStart, coordEnd, p, coordDim, p4estCoordDim, cStart, cEnd, cEndInterior;
3731: DM base;
3733: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3734: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3735: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3736: PetscCall(DMForestGetBaseDM(dm, &base));
3737: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3738: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3739: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3740: p4estCoordDim = PetscMin(coordDim, 3);
3741: for (p = coordStart; p < coordEnd; p++) {
3742: PetscInt *star = NULL, starSize;
3743: PetscInt dof, off, cell = -1, coarsePoint = -1;
3744: PetscInt nCoords, i;
3745: PetscCall(PetscSectionGetDof(coordSec, p, &dof));
3746: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3747: nCoords = dof / coordDim;
3748: PetscCall(PetscSectionGetOffset(coordSec, p, &off));
3749: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3750: for (i = 0; i < starSize; i++) {
3751: PetscInt point = star[2 * i];
3753: if (cStart <= point && point < cEnd) {
3754: cell = point;
3755: break;
3756: }
3757: }
3758: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3759: if (cell >= 0) {
3760: if (cell < cLocalStart) {
3761: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3763: coarsePoint = ghosts[cell].p.which_tree;
3764: } else if (cell < cLocalEnd) {
3765: cell -= cLocalStart;
3766: for (t = flt; t <= llt; t++) {
3767: p4est_tree_t *tree = &(trees[t]);
3769: if (cell >= tree->quadrants_offset && (size_t)cell < tree->quadrants_offset + tree->quadrants.elem_count) {
3770: coarsePoint = t;
3771: break;
3772: }
3773: }
3774: } else {
3775: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3777: coarsePoint = ghosts[cell - cLocalEnd].p.which_tree;
3778: }
3779: }
3780: for (i = 0; i < nCoords; i++) {
3781: PetscScalar *coord = &coords[off + i * coordDim];
3782: PetscReal coordP4est[3] = {0.};
3783: PetscReal coordP4estMapped[3] = {0.};
3784: PetscInt j;
3786: for (j = 0; j < p4estCoordDim; j++) coordP4est[j] = PetscRealPart(coord[j]);
3787: PetscCall((map)(base, coarsePoint, p4estCoordDim, coordP4est, coordP4estMapped, mapCtx));
3788: for (j = 0; j < p4estCoordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3789: }
3790: }
3791: } else { /* we have to transform coordinates back to the unit cube (where geom is defined), and then apply geom */
3792: PetscInt cStart, cEnd, cEndInterior;
3794: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3795: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3796: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3797: if (cLocalStart > 0) {
3798: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3799: PetscInt count;
3801: for (count = 0; count < cLocalStart; count++) {
3802: p4est_quadrant_t *quad = &ghosts[count];
3803: p4est_topidx_t t = quad->p.which_tree;
3805: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, quad, t, pforest->topo->conn, coords));
3806: }
3807: }
3808: for (t = flt; t <= llt; t++) {
3809: p4est_tree_t *tree = &(trees[t]);
3810: PetscInt offset = cLocalStart + tree->quadrants_offset, i;
3811: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3812: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3814: for (i = 0; i < numQuads; i++) {
3815: PetscInt count = i + offset;
3817: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, &quads[i], t, pforest->topo->conn, coords));
3818: }
3819: }
3820: if (cLocalEnd - cLocalStart < cEnd - cStart) {
3821: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3822: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
3823: PetscInt count;
3825: for (count = 0; count < numGhosts - cLocalStart; count++) {
3826: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
3827: p4est_topidx_t t = quad->p.which_tree;
3829: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count + cLocalEnd, quad, t, pforest->topo->conn, coords));
3830: }
3831: }
3832: }
3833: PetscCall(VecRestoreArray(coordLocalVec, &coords));
3834: PetscFunctionReturn(PETSC_SUCCESS);
3835: }
3837: static PetscErrorCode PforestQuadrantIsInterior(p4est_quadrant_t *quad, PetscBool *is_interior)
3838: {
3839: PetscFunctionBegin;
3840: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3841: if ((quad->x > 0) && (quad->x + h < P4EST_ROOT_LEN)
3842: #ifdef P4_TO_P8
3843: && (quad->z > 0) && (quad->z + h < P4EST_ROOT_LEN)
3844: #endif
3845: && (quad->y > 0) && (quad->y + h < P4EST_ROOT_LEN)) {
3846: *is_interior = PETSC_TRUE;
3847: } else {
3848: *is_interior = PETSC_FALSE;
3849: }
3850: PetscFunctionReturn(PETSC_SUCCESS);
3851: }
3853: /* We always use DG coordinates with p4est: if they do not match the vertex
3854: coordinates, add space for them in the section */
3855: 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)
3856: {
3857: PetscBool is_interior;
3859: PetscFunctionBegin;
3860: PetscCall(PforestQuadrantIsInterior(quad, &is_interior));
3861: if (is_interior) { // quads in the interior of a coarse cell can't touch periodic interfaces
3862: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3863: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3864: } else {
3865: PetscInt cSize;
3866: PetscScalar *values = NULL;
3867: PetscBool same_coords = PETSC_TRUE;
3869: PetscCall(DMPlexVecGetClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3870: PetscAssert(cSize == cDim * P4EST_CHILDREN, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected closure size");
3871: for (int c = 0; c < P4EST_CHILDREN; c++) {
3872: p4est_qcoord_t quad_coords[3];
3873: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3874: double corner_coords[3];
3875: double vert_coords[3];
3876: PetscInt corner = PetscVertToP4estVert[c];
3878: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) vert_coords[d] = PetscRealPart(values[c * cDim + d]);
3880: quad_coords[0] = quad->x;
3881: quad_coords[1] = quad->y;
3882: #ifdef P4_TO_P8
3883: quad_coords[2] = quad->z;
3884: #endif
3885: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3886: #ifndef P4_TO_P8
3887: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3888: #else
3889: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3890: #endif
3891: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) {
3892: if (fabs(vert_coords[d] - corner_coords[d]) > PETSC_SMALL) {
3893: same_coords = PETSC_FALSE;
3894: break;
3895: }
3896: }
3897: }
3898: if (same_coords) {
3899: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3900: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3901: } else {
3902: PetscCall(PetscSectionSetDof(newSection, cell, cSize));
3903: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, cSize));
3904: }
3905: PetscCall(DMPlexVecRestoreClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3906: }
3907: PetscFunctionReturn(PETSC_SUCCESS);
3908: }
3910: static PetscErrorCode PforestLocalizeCell(DM plex, PetscInt cDim, DM_Forest_pforest *pforest, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad, PetscScalar coords[])
3911: {
3912: PetscInt cdof, off;
3914: PetscFunctionBegin;
3915: PetscCall(PetscSectionGetDof(newSection, cell, &cdof));
3916: if (!cdof) PetscFunctionReturn(PETSC_SUCCESS);
3918: PetscCall(PetscSectionGetOffset(newSection, cell, &off));
3919: for (PetscInt c = 0, pos = off; c < P4EST_CHILDREN; c++) {
3920: p4est_qcoord_t quad_coords[3];
3921: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3922: double corner_coords[3];
3923: PetscInt corner = PetscVertToP4estVert[c];
3925: quad_coords[0] = quad->x;
3926: quad_coords[1] = quad->y;
3927: #ifdef P4_TO_P8
3928: quad_coords[2] = quad->z;
3929: #endif
3930: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3931: #ifndef P4_TO_P8
3932: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3933: #else
3934: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3935: #endif
3936: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) coords[pos++] = corner_coords[d];
3937: for (PetscInt d = PetscMin(cDim, 3); d < cDim; d++) coords[pos++] = 0.;
3938: }
3939: PetscFunctionReturn(PETSC_SUCCESS);
3940: }
3942: static PetscErrorCode DMPforestLocalizeCoordinates(DM dm, DM plex)
3943: {
3944: DM_Forest *forest;
3945: DM_Forest_pforest *pforest;
3946: DM base, cdm, cdmCell;
3947: Vec cVec, cVecOld;
3948: PetscSection oldSection, newSection;
3949: PetscScalar *coords2;
3950: const PetscReal *L;
3951: PetscInt cLocalStart, cLocalEnd, coarsePoint;
3952: PetscInt cDim, newStart, newEnd;
3953: PetscInt v, vStart, vEnd, cp, cStart, cEnd, cEndInterior;
3954: p4est_topidx_t flt, llt, t;
3955: p4est_tree_t *trees;
3956: PetscBool baseLocalized = PETSC_FALSE;
3958: PetscFunctionBegin;
3959: PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
3960: /* we localize on all cells if we don't have a base DM or the base DM coordinates have not been localized */
3961: PetscCall(DMGetCoordinateDim(dm, &cDim));
3962: PetscCall(DMForestGetBaseDM(dm, &base));
3963: if (base) PetscCall(DMGetCoordinatesLocalized(base, &baseLocalized));
3964: if (!baseLocalized) base = NULL;
3965: if (!baseLocalized && !L) PetscFunctionReturn(PETSC_SUCCESS);
3966: PetscCall(DMPlexGetChart(plex, &newStart, &newEnd));
3968: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newSection));
3969: PetscCall(PetscSectionSetNumFields(newSection, 1));
3970: PetscCall(PetscSectionSetFieldComponents(newSection, 0, cDim));
3971: PetscCall(PetscSectionSetChart(newSection, newStart, newEnd));
3973: PetscCall(DMGetCoordinateSection(plex, &oldSection));
3974: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3975: PetscCall(DMGetCoordinatesLocal(plex, &cVecOld));
3977: forest = (DM_Forest *)dm->data;
3978: pforest = (DM_Forest_pforest *)forest->data;
3979: cLocalStart = pforest->cLocalStart;
3980: cLocalEnd = pforest->cLocalEnd;
3981: flt = pforest->forest->first_local_tree;
3982: llt = pforest->forest->last_local_tree;
3983: trees = (p4est_tree_t *)pforest->forest->trees->array;
3985: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3986: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3987: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3988: cp = 0;
3989: if (cLocalStart > 0) {
3990: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3991: PetscInt cell;
3993: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
3994: p4est_quadrant_t *quad = &ghosts[cell];
3996: coarsePoint = quad->p.which_tree;
3997: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
3998: }
3999: }
4000: for (t = flt; t <= llt; t++) {
4001: p4est_tree_t *tree = &(trees[t]);
4002: PetscInt offset = cLocalStart + tree->quadrants_offset;
4003: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
4004: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4005: PetscInt i;
4007: if (!numQuads) continue;
4008: coarsePoint = t;
4009: for (i = 0; i < numQuads; i++, cp++) {
4010: PetscInt cell = i + offset;
4011: p4est_quadrant_t *quad = &quads[i];
4013: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4014: }
4015: }
4016: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4017: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4018: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4019: PetscInt count;
4021: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4022: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4023: coarsePoint = quad->p.which_tree;
4024: PetscInt cell = count + cLocalEnd;
4026: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4027: }
4028: }
4029: PetscAssert(cp == cEnd - cStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of fine cells %" PetscInt_FMT " != %" PetscInt_FMT, cp, cEnd - cStart);
4031: PetscCall(PetscSectionSetUp(newSection));
4032: PetscCall(DMGetCoordinateDM(plex, &cdm));
4033: PetscCall(DMClone(cdm, &cdmCell));
4034: PetscCall(DMSetCellCoordinateDM(plex, cdmCell));
4035: PetscCall(DMDestroy(&cdmCell));
4036: PetscCall(DMSetCellCoordinateSection(plex, cDim, newSection));
4037: PetscCall(PetscSectionGetStorageSize(newSection, &v));
4038: PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
4039: PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
4040: PetscCall(VecSetBlockSize(cVec, cDim));
4041: PetscCall(VecSetSizes(cVec, v, PETSC_DETERMINE));
4042: PetscCall(VecSetType(cVec, VECSTANDARD));
4043: PetscCall(VecSet(cVec, PETSC_MIN_REAL));
4045: /* Localize coordinates on cells if needed */
4046: PetscCall(VecGetArray(cVec, &coords2));
4047: cp = 0;
4048: if (cLocalStart > 0) {
4049: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4050: PetscInt cell;
4052: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
4053: p4est_quadrant_t *quad = &ghosts[cell];
4055: coarsePoint = quad->p.which_tree;
4056: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4057: }
4058: }
4059: for (t = flt; t <= llt; t++) {
4060: p4est_tree_t *tree = &(trees[t]);
4061: PetscInt offset = cLocalStart + tree->quadrants_offset;
4062: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
4063: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4064: PetscInt i;
4066: if (!numQuads) continue;
4067: coarsePoint = t;
4068: for (i = 0; i < numQuads; i++, cp++) {
4069: PetscInt cell = i + offset;
4070: p4est_quadrant_t *quad = &quads[i];
4072: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4073: }
4074: }
4075: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4076: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4077: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4078: PetscInt count;
4080: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4081: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4082: coarsePoint = quad->p.which_tree;
4083: PetscInt cell = count + cLocalEnd;
4085: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4086: }
4087: }
4088: PetscCall(VecRestoreArray(cVec, &coords2));
4089: PetscCall(DMSetCellCoordinatesLocal(plex, cVec));
4090: PetscCall(VecDestroy(&cVec));
4091: PetscCall(PetscSectionDestroy(&newSection));
4092: PetscFunctionReturn(PETSC_SUCCESS);
4093: }
4095: #define DMForestClearAdaptivityForest_pforest _append_pforest(DMForestClearAdaptivityForest)
4096: static PetscErrorCode DMForestClearAdaptivityForest_pforest(DM dm)
4097: {
4098: DM_Forest *forest;
4099: DM_Forest_pforest *pforest;
4101: PetscFunctionBegin;
4102: forest = (DM_Forest *)dm->data;
4103: pforest = (DM_Forest_pforest *)forest->data;
4104: PetscCall(PetscSFDestroy(&(pforest->pointAdaptToSelfSF)));
4105: PetscCall(PetscSFDestroy(&(pforest->pointSelfToAdaptSF)));
4106: PetscCall(PetscFree(pforest->pointAdaptToSelfCids));
4107: PetscCall(PetscFree(pforest->pointSelfToAdaptCids));
4108: PetscFunctionReturn(PETSC_SUCCESS);
4109: }
4111: static PetscErrorCode DMConvert_pforest_plex(DM dm, DMType newtype, DM *plex)
4112: {
4113: DM_Forest *forest;
4114: DM_Forest_pforest *pforest;
4115: DM refTree, newPlex, base;
4116: PetscInt adjDim, adjCodim, coordDim;
4117: MPI_Comm comm;
4118: PetscBool isPforest;
4119: PetscInt dim;
4120: PetscInt overlap;
4121: p4est_connect_type_t ctype;
4122: p4est_locidx_t first_local_quad = -1;
4123: sc_array_t *points_per_dim, *cone_sizes, *cones, *cone_orientations, *coords, *children, *parents, *childids, *leaves, *remotes;
4124: PetscSection parentSection;
4125: PetscSF pointSF;
4126: size_t zz, count;
4127: PetscInt pStart, pEnd;
4128: DMLabel ghostLabelBase = NULL;
4130: PetscFunctionBegin;
4133: comm = PetscObjectComm((PetscObject)dm);
4134: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPFOREST, &isPforest));
4135: PetscCheck(isPforest, comm, PETSC_ERR_ARG_WRONG, "Expected DM type %s, got %s", DMPFOREST, ((PetscObject)dm)->type_name);
4136: PetscCall(DMGetDimension(dm, &dim));
4137: PetscCheck(dim == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Expected DM dimension %d, got %" PetscInt_FMT, P4EST_DIM, dim);
4138: forest = (DM_Forest *)dm->data;
4139: pforest = (DM_Forest_pforest *)forest->data;
4140: PetscCall(DMForestGetBaseDM(dm, &base));
4141: if (base) PetscCall(DMGetLabel(base, "ghost", &ghostLabelBase));
4142: if (!pforest->plex) {
4143: PetscMPIInt size;
4145: PetscCallMPI(MPI_Comm_size(comm, &size));
4146: PetscCall(DMCreate(comm, &newPlex));
4147: PetscCall(DMSetType(newPlex, DMPLEX));
4148: PetscCall(DMSetMatType(newPlex, dm->mattype));
4149: /* share labels */
4150: PetscCall(DMCopyLabels(dm, newPlex, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4151: PetscCall(DMForestGetAdjacencyDimension(dm, &adjDim));
4152: PetscCall(DMForestGetAdjacencyCodimension(dm, &adjCodim));
4153: PetscCall(DMGetCoordinateDim(dm, &coordDim));
4154: if (adjDim == 0) {
4155: ctype = P4EST_CONNECT_FULL;
4156: } else if (adjCodim == 1) {
4157: ctype = P4EST_CONNECT_FACE;
4158: #if defined(P4_TO_P8)
4159: } else if (adjDim == 1) {
4160: ctype = P8EST_CONNECT_EDGE;
4161: #endif
4162: } else {
4163: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid adjacency dimension %" PetscInt_FMT, adjDim);
4164: }
4165: PetscCheck(ctype == P4EST_CONNECT_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Adjacency dimension %" PetscInt_FMT " / codimension %" PetscInt_FMT " not supported yet", adjDim, adjCodim);
4166: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
4167: PetscCall(DMPlexSetOverlap_Plex(newPlex, NULL, overlap));
4169: points_per_dim = sc_array_new(sizeof(p4est_locidx_t));
4170: cone_sizes = sc_array_new(sizeof(p4est_locidx_t));
4171: cones = sc_array_new(sizeof(p4est_locidx_t));
4172: cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
4173: coords = sc_array_new(3 * sizeof(double));
4174: children = sc_array_new(sizeof(p4est_locidx_t));
4175: parents = sc_array_new(sizeof(p4est_locidx_t));
4176: childids = sc_array_new(sizeof(p4est_locidx_t));
4177: leaves = sc_array_new(sizeof(p4est_locidx_t));
4178: remotes = sc_array_new(2 * sizeof(p4est_locidx_t));
4180: 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));
4182: pforest->cLocalStart = (PetscInt)first_local_quad;
4183: pforest->cLocalEnd = pforest->cLocalStart + (PetscInt)pforest->forest->local_num_quadrants;
4184: PetscCall(locidx_to_PetscInt(points_per_dim));
4185: PetscCall(locidx_to_PetscInt(cone_sizes));
4186: PetscCall(locidx_to_PetscInt(cones));
4187: PetscCall(locidx_to_PetscInt(cone_orientations));
4188: PetscCall(coords_double_to_PetscScalar(coords, coordDim));
4189: PetscCall(locidx_to_PetscInt(children));
4190: PetscCall(locidx_to_PetscInt(parents));
4191: PetscCall(locidx_to_PetscInt(childids));
4192: PetscCall(locidx_to_PetscInt(leaves));
4193: PetscCall(locidx_pair_to_PetscSFNode(remotes));
4195: PetscCall(DMSetDimension(newPlex, P4EST_DIM));
4196: PetscCall(DMSetCoordinateDim(newPlex, coordDim));
4197: PetscCall(DMPlexSetMaxProjectionHeight(newPlex, P4EST_DIM - 1));
4198: PetscCall(DMPlexCreateFromDAG(newPlex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array));
4199: PetscCall(DMPlexConvertOldOrientations_Internal(newPlex));
4200: PetscCall(DMCreateReferenceTree_pforest(comm, &refTree));
4201: PetscCall(DMPlexSetReferenceTree(newPlex, refTree));
4202: PetscCall(PetscSectionCreate(comm, &parentSection));
4203: PetscCall(DMPlexGetChart(newPlex, &pStart, &pEnd));
4204: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
4205: count = children->elem_count;
4206: for (zz = 0; zz < count; zz++) {
4207: PetscInt child = *((PetscInt *)sc_array_index(children, zz));
4209: PetscCall(PetscSectionSetDof(parentSection, child, 1));
4210: }
4211: PetscCall(PetscSectionSetUp(parentSection));
4212: PetscCall(DMPlexSetTree(newPlex, parentSection, (PetscInt *)parents->array, (PetscInt *)childids->array));
4213: PetscCall(PetscSectionDestroy(&parentSection));
4214: PetscCall(PetscSFCreate(comm, &pointSF));
4215: /*
4216: These arrays defining the sf are from the p4est library, but the code there shows the leaves being populated in increasing order.
4217: https://gitlab.com/petsc/petsc/merge_requests/2248#note_240186391
4218: */
4219: PetscCall(PetscSFSetGraph(pointSF, pEnd - pStart, (PetscInt)leaves->elem_count, (PetscInt *)leaves->array, PETSC_COPY_VALUES, (PetscSFNode *)remotes->array, PETSC_COPY_VALUES));
4220: PetscCall(DMSetPointSF(newPlex, pointSF));
4221: PetscCall(DMSetPointSF(dm, pointSF));
4222: {
4223: DM coordDM;
4225: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4226: PetscCall(DMSetPointSF(coordDM, pointSF));
4227: }
4228: PetscCall(PetscSFDestroy(&pointSF));
4229: sc_array_destroy(points_per_dim);
4230: sc_array_destroy(cone_sizes);
4231: sc_array_destroy(cones);
4232: sc_array_destroy(cone_orientations);
4233: sc_array_destroy(coords);
4234: sc_array_destroy(children);
4235: sc_array_destroy(parents);
4236: sc_array_destroy(childids);
4237: sc_array_destroy(leaves);
4238: sc_array_destroy(remotes);
4240: {
4241: const PetscReal *maxCell, *Lstart, *L;
4243: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
4244: PetscCall(DMSetPeriodicity(newPlex, maxCell, Lstart, L));
4245: PetscCall(DMPforestLocalizeCoordinates(dm, newPlex));
4246: }
4248: if (overlap > 0) { /* the p4est routine can't set all of the coordinates in its routine if there is overlap */
4249: Vec coordsGlobal, coordsLocal;
4250: const PetscScalar *globalArray;
4251: PetscScalar *localArray;
4252: PetscSF coordSF;
4253: DM coordDM;
4255: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4256: PetscCall(DMGetSectionSF(coordDM, &coordSF));
4257: PetscCall(DMGetCoordinates(newPlex, &coordsGlobal));
4258: PetscCall(DMGetCoordinatesLocal(newPlex, &coordsLocal));
4259: PetscCall(VecGetArrayRead(coordsGlobal, &globalArray));
4260: PetscCall(VecGetArray(coordsLocal, &localArray));
4261: PetscCall(PetscSFBcastBegin(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4262: PetscCall(PetscSFBcastEnd(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4263: PetscCall(VecRestoreArray(coordsLocal, &localArray));
4264: PetscCall(VecRestoreArrayRead(coordsGlobal, &globalArray));
4265: PetscCall(DMSetCoordinatesLocal(newPlex, coordsLocal));
4266: }
4267: PetscCall(DMPforestMapCoordinates(dm, newPlex));
4269: pforest->plex = newPlex;
4271: /* copy labels */
4272: PetscCall(DMPforestLabelsFinalize(dm, newPlex));
4274: if (ghostLabelBase || pforest->ghostName) { /* we have to do this after copying labels because the labels drive the construction of ghost cells */
4275: PetscInt numAdded;
4276: DM newPlexGhosted;
4277: void *ctx;
4279: PetscCall(DMPlexConstructGhostCells(newPlex, pforest->ghostName, &numAdded, &newPlexGhosted));
4280: PetscCall(DMGetApplicationContext(newPlex, &ctx));
4281: PetscCall(DMSetApplicationContext(newPlexGhosted, ctx));
4282: /* we want the sf for the ghost dm to be the one for the p4est dm as well */
4283: PetscCall(DMGetPointSF(newPlexGhosted, &pointSF));
4284: PetscCall(DMSetPointSF(dm, pointSF));
4285: PetscCall(DMDestroy(&newPlex));
4286: PetscCall(DMPlexSetReferenceTree(newPlexGhosted, refTree));
4287: PetscCall(DMForestClearAdaptivityForest_pforest(dm));
4288: newPlex = newPlexGhosted;
4290: /* share the labels back */
4291: PetscCall(DMDestroyLabelLinkList_Internal(dm));
4292: PetscCall(DMCopyLabels(newPlex, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4293: pforest->plex = newPlex;
4294: }
4295: PetscCall(DMDestroy(&refTree));
4296: if (dm->setfromoptionscalled) {
4297: PetscObjectOptionsBegin((PetscObject)newPlex);
4298: PetscCall(DMSetFromOptions_NonRefinement_Plex(newPlex, PetscOptionsObject));
4299: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)newPlex, PetscOptionsObject));
4300: PetscOptionsEnd();
4301: }
4302: PetscCall(DMViewFromOptions(newPlex, NULL, "-dm_p4est_plex_view"));
4303: {
4304: DM cdm;
4305: PetscSection coordsSec;
4306: Vec coords;
4307: PetscInt cDim;
4309: PetscCall(DMGetCoordinateDim(newPlex, &cDim));
4310: PetscCall(DMGetCoordinateSection(newPlex, &coordsSec));
4311: PetscCall(DMSetCoordinateSection(dm, cDim, coordsSec));
4312: PetscCall(DMGetCoordinatesLocal(newPlex, &coords));
4313: PetscCall(DMSetCoordinatesLocal(dm, coords));
4314: PetscCall(DMGetCellCoordinateDM(newPlex, &cdm));
4315: if (cdm) PetscCall(DMSetCellCoordinateDM(dm, cdm));
4316: PetscCall(DMGetCellCoordinateSection(newPlex, &coordsSec));
4317: if (coordsSec) PetscCall(DMSetCellCoordinateSection(dm, cDim, coordsSec));
4318: PetscCall(DMGetCellCoordinatesLocal(newPlex, &coords));
4319: if (coords) PetscCall(DMSetCellCoordinatesLocal(dm, coords));
4320: }
4321: }
4322: newPlex = pforest->plex;
4323: if (plex) {
4324: PetscCall(DMClone(newPlex, plex));
4325: #if 0
4326: PetscCall(DMGetCoordinateDM(newPlex,&coordDM));
4327: PetscCall(DMSetCoordinateDM(*plex,coordDM));
4328: PetscCall(DMGetCellCoordinateDM(newPlex,&coordDM));
4329: PetscCall(DMSetCellCoordinateDM(*plex,coordDM));
4330: #endif
4331: PetscCall(DMShareDiscretization(dm, *plex));
4332: }
4333: PetscFunctionReturn(PETSC_SUCCESS);
4334: }
4336: static PetscErrorCode DMSetFromOptions_pforest(DM dm, PetscOptionItems *PetscOptionsObject)
4337: {
4338: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4339: char stringBuffer[256];
4340: PetscBool flg;
4342: PetscFunctionBegin;
4343: PetscCall(DMSetFromOptions_Forest(dm, PetscOptionsObject));
4344: PetscOptionsHeadBegin(PetscOptionsObject, "DM" P4EST_STRING " options");
4345: PetscCall(PetscOptionsBool("-dm_p4est_partition_for_coarsening", "partition forest to allow for coarsening", "DMP4estSetPartitionForCoarsening", pforest->partition_for_coarsening, &(pforest->partition_for_coarsening), NULL));
4346: PetscCall(PetscOptionsString("-dm_p4est_ghost_label_name", "the name of the ghost label when converting from a DMPlex", NULL, NULL, stringBuffer, sizeof(stringBuffer), &flg));
4347: PetscOptionsHeadEnd();
4348: if (flg) {
4349: PetscCall(PetscFree(pforest->ghostName));
4350: PetscCall(PetscStrallocpy(stringBuffer, &pforest->ghostName));
4351: }
4352: PetscFunctionReturn(PETSC_SUCCESS);
4353: }
4355: #if !defined(P4_TO_P8)
4356: #define DMPforestGetPartitionForCoarsening DMP4estGetPartitionForCoarsening
4357: #define DMPforestSetPartitionForCoarsening DMP4estSetPartitionForCoarsening
4358: #else
4359: #define DMPforestGetPartitionForCoarsening DMP8estGetPartitionForCoarsening
4360: #define DMPforestSetPartitionForCoarsening DMP8estSetPartitionForCoarsening
4361: #endif
4363: PETSC_EXTERN PetscErrorCode DMPforestGetPartitionForCoarsening(DM dm, PetscBool *flg)
4364: {
4365: DM_Forest_pforest *pforest;
4367: PetscFunctionBegin;
4369: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4370: *flg = pforest->partition_for_coarsening;
4371: PetscFunctionReturn(PETSC_SUCCESS);
4372: }
4374: PETSC_EXTERN PetscErrorCode DMPforestSetPartitionForCoarsening(DM dm, PetscBool flg)
4375: {
4376: DM_Forest_pforest *pforest;
4378: PetscFunctionBegin;
4380: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4381: pforest->partition_for_coarsening = flg;
4382: PetscFunctionReturn(PETSC_SUCCESS);
4383: }
4385: static PetscErrorCode DMPforestGetPlex(DM dm, DM *plex)
4386: {
4387: DM_Forest_pforest *pforest;
4389: PetscFunctionBegin;
4390: if (plex) *plex = NULL;
4391: PetscCall(DMSetUp(dm));
4392: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4393: if (!pforest->plex) PetscCall(DMConvert_pforest_plex(dm, DMPLEX, NULL));
4394: PetscCall(DMShareDiscretization(dm, pforest->plex));
4395: if (plex) *plex = pforest->plex;
4396: PetscFunctionReturn(PETSC_SUCCESS);
4397: }
4399: #define DMCreateInterpolation_pforest _append_pforest(DMCreateInterpolation)
4400: static PetscErrorCode DMCreateInterpolation_pforest(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
4401: {
4402: PetscSection gsc, gsf;
4403: PetscInt m, n;
4404: DM cdm;
4406: PetscFunctionBegin;
4407: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4408: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
4409: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4410: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &n));
4412: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), interpolation));
4413: PetscCall(MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4414: PetscCall(MatSetType(*interpolation, MATAIJ));
4416: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4417: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only interpolation from coarse DM for now");
4419: {
4420: DM plexF, plexC;
4421: PetscSF sf;
4422: PetscInt *cids;
4423: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4425: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4426: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4427: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4428: PetscCall(PetscSFSetUp(sf));
4429: PetscCall(DMPlexComputeInterpolatorTree(plexC, plexF, sf, cids, *interpolation));
4430: PetscCall(PetscSFDestroy(&sf));
4431: PetscCall(PetscFree(cids));
4432: }
4433: PetscCall(MatViewFromOptions(*interpolation, NULL, "-interp_mat_view"));
4434: /* Use naive scaling */
4435: PetscCall(DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling));
4436: PetscFunctionReturn(PETSC_SUCCESS);
4437: }
4439: #define DMCreateInjection_pforest _append_pforest(DMCreateInjection)
4440: static PetscErrorCode DMCreateInjection_pforest(DM dmCoarse, DM dmFine, Mat *injection)
4441: {
4442: PetscSection gsc, gsf;
4443: PetscInt m, n;
4444: DM cdm;
4446: PetscFunctionBegin;
4447: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4448: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &n));
4449: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4450: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &m));
4452: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), injection));
4453: PetscCall(MatSetSizes(*injection, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4454: PetscCall(MatSetType(*injection, MATAIJ));
4456: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4457: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only injection to coarse DM for now");
4459: {
4460: DM plexF, plexC;
4461: PetscSF sf;
4462: PetscInt *cids;
4463: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4465: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4466: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4467: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4468: PetscCall(PetscSFSetUp(sf));
4469: PetscCall(DMPlexComputeInjectorTree(plexC, plexF, sf, cids, *injection));
4470: PetscCall(PetscSFDestroy(&sf));
4471: PetscCall(PetscFree(cids));
4472: }
4473: PetscCall(MatViewFromOptions(*injection, NULL, "-inject_mat_view"));
4474: /* Use naive scaling */
4475: PetscFunctionReturn(PETSC_SUCCESS);
4476: }
4478: #define DMForestTransferVecFromBase_pforest _append_pforest(DMForestTransferVecFromBase)
4479: static PetscErrorCode DMForestTransferVecFromBase_pforest(DM dm, Vec vecIn, Vec vecOut)
4480: {
4481: DM dmIn, dmVecIn, base, basec, plex, coarseDM;
4482: DM *hierarchy;
4483: PetscSF sfRed = NULL;
4484: PetscDS ds;
4485: Vec vecInLocal, vecOutLocal;
4486: DMLabel subpointMap;
4487: PetscInt minLevel, mh, n_hi, i;
4488: PetscBool hiforest, *hierarchy_forest;
4490: PetscFunctionBegin;
4491: PetscCall(VecGetDM(vecIn, &dmVecIn));
4492: PetscCall(DMGetDS(dmVecIn, &ds));
4493: PetscCheck(ds, PetscObjectComm((PetscObject)dmVecIn), PETSC_ERR_SUP, "Cannot transfer without a PetscDS object");
4494: { /* we cannot stick user contexts into function callbacks for DMProjectFieldLocal! */
4495: PetscSection section;
4496: PetscInt Nf;
4498: PetscCall(DMGetLocalSection(dmVecIn, §ion));
4499: PetscCall(PetscSectionGetNumFields(section, &Nf));
4500: PetscCheck(Nf <= 3, PetscObjectComm((PetscObject)dmVecIn), PETSC_ERR_SUP, "Number of fields %" PetscInt_FMT " are currently not supported! Send an email at petsc-dev@mcs.anl.gov", Nf);
4501: }
4502: PetscCall(DMForestGetMinimumRefinement(dm, &minLevel));
4503: PetscCheck(!minLevel, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot transfer with minimum refinement set to %" PetscInt_FMT ". Rerun with DMForestSetMinimumRefinement(dm,0)", minLevel);
4504: PetscCall(DMForestGetBaseDM(dm, &base));
4505: PetscCheck(base, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing base DM");
4507: PetscCall(VecSet(vecOut, 0.0));
4508: if (dmVecIn == base) { /* sequential runs */
4509: PetscCall(PetscObjectReference((PetscObject)vecIn));
4510: } else {
4511: PetscSection secIn, secInRed;
4512: Vec vecInRed, vecInLocal;
4514: PetscCall(PetscObjectQuery((PetscObject)base, "_base_migration_sf", (PetscObject *)&sfRed));
4515: PetscCheck(sfRed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not the DM set with DMForestSetBaseDM()");
4516: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmVecIn), &secInRed));
4517: PetscCall(VecCreate(PETSC_COMM_SELF, &vecInRed));
4518: PetscCall(DMGetLocalSection(dmVecIn, &secIn));
4519: PetscCall(DMGetLocalVector(dmVecIn, &vecInLocal));
4520: PetscCall(DMGlobalToLocalBegin(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4521: PetscCall(DMGlobalToLocalEnd(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4522: PetscCall(DMPlexDistributeField(dmVecIn, sfRed, secIn, vecInLocal, secInRed, vecInRed));
4523: PetscCall(DMRestoreLocalVector(dmVecIn, &vecInLocal));
4524: PetscCall(PetscSectionDestroy(&secInRed));
4525: vecIn = vecInRed;
4526: }
4528: /* we first search through the AdaptivityForest hierarchy
4529: once we found the first disconnected forest, we upsweep the DM hierarchy */
4530: hiforest = PETSC_TRUE;
4532: /* upsweep to the coarsest DM */
4533: n_hi = 0;
4534: coarseDM = dm;
4535: do {
4536: PetscBool isforest;
4538: dmIn = coarseDM;
4539: /* need to call DMSetUp to have the hierarchy recursively setup */
4540: PetscCall(DMSetUp(dmIn));
4541: PetscCall(DMIsForest(dmIn, &isforest));
4542: PetscCheck(isforest, PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Cannot currently transfer through a mixed hierarchy! Found DM type %s", ((PetscObject)dmIn)->type_name);
4543: coarseDM = NULL;
4544: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4545: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4546: hiforest = PETSC_FALSE;
4547: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4548: }
4549: n_hi++;
4550: } while (coarseDM);
4552: PetscCall(PetscMalloc2(n_hi, &hierarchy, n_hi, &hierarchy_forest));
4554: i = 0;
4555: hiforest = PETSC_TRUE;
4556: coarseDM = dm;
4557: do {
4558: dmIn = coarseDM;
4559: coarseDM = NULL;
4560: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4561: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4562: hiforest = PETSC_FALSE;
4563: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4564: }
4565: i++;
4566: hierarchy[n_hi - i] = dmIn;
4567: } while (coarseDM);
4569: /* project base vector on the coarsest forest (minimum refinement = 0) */
4570: PetscCall(DMPforestGetPlex(dmIn, &plex));
4572: /* Check this plex is compatible with the base */
4573: {
4574: IS gnum[2];
4575: PetscInt ncells[2], gncells[2];
4577: PetscCall(DMPlexGetCellNumbering(base, &gnum[0]));
4578: PetscCall(DMPlexGetCellNumbering(plex, &gnum[1]));
4579: PetscCall(ISGetMinMax(gnum[0], NULL, &ncells[0]));
4580: PetscCall(ISGetMinMax(gnum[1], NULL, &ncells[1]));
4581: PetscCall(MPIU_Allreduce(ncells, gncells, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
4582: PetscCheck(gncells[0] == gncells[1], PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Invalid number of base cells! Expected %" PetscInt_FMT ", found %" PetscInt_FMT, gncells[0] + 1, gncells[1] + 1);
4583: }
4585: PetscCall(DMGetLabel(dmIn, "_forest_base_subpoint_map", &subpointMap));
4586: PetscCheck(subpointMap, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing _forest_base_subpoint_map label");
4588: PetscCall(DMPlexGetMaxProjectionHeight(base, &mh));
4589: PetscCall(DMPlexSetMaxProjectionHeight(plex, mh));
4591: PetscCall(DMClone(base, &basec));
4592: PetscCall(DMCopyDisc(dmVecIn, basec));
4593: if (sfRed) {
4594: PetscCall(PetscObjectReference((PetscObject)vecIn));
4595: vecInLocal = vecIn;
4596: } else {
4597: PetscCall(DMCreateLocalVector(basec, &vecInLocal));
4598: PetscCall(DMGlobalToLocalBegin(basec, vecIn, INSERT_VALUES, vecInLocal));
4599: PetscCall(DMGlobalToLocalEnd(basec, vecIn, INSERT_VALUES, vecInLocal));
4600: }
4602: PetscCall(DMGetLocalVector(dmIn, &vecOutLocal));
4603: { /* get degrees of freedom ordered onto dmIn */
4604: PetscSF basetocoarse;
4605: PetscInt bStart, bEnd, nroots;
4606: PetscInt iStart, iEnd, nleaves, leaf;
4607: PetscMPIInt rank;
4608: PetscSFNode *remotes;
4609: PetscSection secIn, secOut;
4610: PetscInt *remoteOffsets;
4611: PetscSF transferSF;
4612: const PetscScalar *inArray;
4613: PetscScalar *outArray;
4615: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)basec), &rank));
4616: PetscCall(DMPlexGetChart(basec, &bStart, &bEnd));
4617: nroots = PetscMax(bEnd - bStart, 0);
4618: PetscCall(DMPlexGetChart(plex, &iStart, &iEnd));
4619: nleaves = PetscMax(iEnd - iStart, 0);
4621: PetscCall(PetscMalloc1(nleaves, &remotes));
4622: for (leaf = iStart; leaf < iEnd; leaf++) {
4623: PetscInt index;
4625: remotes[leaf - iStart].rank = rank;
4626: PetscCall(DMLabelGetValue(subpointMap, leaf, &index));
4627: remotes[leaf - iStart].index = index;
4628: }
4630: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)basec), &basetocoarse));
4631: PetscCall(PetscSFSetGraph(basetocoarse, nroots, nleaves, NULL, PETSC_OWN_POINTER, remotes, PETSC_OWN_POINTER));
4632: PetscCall(PetscSFSetUp(basetocoarse));
4633: PetscCall(DMGetLocalSection(basec, &secIn));
4634: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmIn), &secOut));
4635: PetscCall(PetscSFDistributeSection(basetocoarse, secIn, &remoteOffsets, secOut));
4636: PetscCall(PetscSFCreateSectionSF(basetocoarse, secIn, remoteOffsets, secOut, &transferSF));
4637: PetscCall(PetscFree(remoteOffsets));
4638: PetscCall(VecGetArrayWrite(vecOutLocal, &outArray));
4639: PetscCall(VecGetArrayRead(vecInLocal, &inArray));
4640: PetscCall(PetscSFBcastBegin(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4641: PetscCall(PetscSFBcastEnd(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4642: PetscCall(VecRestoreArrayRead(vecInLocal, &inArray));
4643: PetscCall(VecRestoreArrayWrite(vecOutLocal, &outArray));
4644: PetscCall(PetscSFDestroy(&transferSF));
4645: PetscCall(PetscSectionDestroy(&secOut));
4646: PetscCall(PetscSFDestroy(&basetocoarse));
4647: }
4648: PetscCall(VecDestroy(&vecInLocal));
4649: PetscCall(DMDestroy(&basec));
4650: PetscCall(VecDestroy(&vecIn));
4652: /* output */
4653: if (n_hi > 1) { /* downsweep the stored hierarchy */
4654: Vec vecOut1, vecOut2;
4655: DM fineDM;
4657: PetscCall(DMGetGlobalVector(dmIn, &vecOut1));
4658: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut1));
4659: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4660: for (i = 1; i < n_hi - 1; i++) {
4661: fineDM = hierarchy[i];
4662: PetscCall(DMGetGlobalVector(fineDM, &vecOut2));
4663: PetscCall(DMForestTransferVec(dmIn, vecOut1, fineDM, vecOut2, PETSC_TRUE, 0.0));
4664: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4665: vecOut1 = vecOut2;
4666: dmIn = fineDM;
4667: }
4668: PetscCall(DMForestTransferVec(dmIn, vecOut1, dm, vecOut, PETSC_TRUE, 0.0));
4669: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4670: } else {
4671: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut));
4672: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4673: }
4674: PetscCall(PetscFree2(hierarchy, hierarchy_forest));
4675: PetscFunctionReturn(PETSC_SUCCESS);
4676: }
4678: #define DMForestTransferVec_pforest _append_pforest(DMForestTransferVec)
4679: static PetscErrorCode DMForestTransferVec_pforest(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscBool useBCs, PetscReal time)
4680: {
4681: DM adaptIn, adaptOut, plexIn, plexOut;
4682: DM_Forest *forestIn, *forestOut, *forestAdaptIn, *forestAdaptOut;
4683: PetscInt dofPerDim[] = {1, 1, 1, 1};
4684: PetscSF inSF = NULL, outSF = NULL;
4685: PetscInt *inCids = NULL, *outCids = NULL;
4686: DMAdaptFlag purposeIn, purposeOut;
4688: PetscFunctionBegin;
4689: forestOut = (DM_Forest *)dmOut->data;
4690: forestIn = (DM_Forest *)dmIn->data;
4692: PetscCall(DMForestGetAdaptivityForest(dmOut, &adaptOut));
4693: PetscCall(DMForestGetAdaptivityPurpose(dmOut, &purposeOut));
4694: forestAdaptOut = adaptOut ? (DM_Forest *)adaptOut->data : NULL;
4696: PetscCall(DMForestGetAdaptivityForest(dmIn, &adaptIn));
4697: PetscCall(DMForestGetAdaptivityPurpose(dmIn, &purposeIn));
4698: forestAdaptIn = adaptIn ? (DM_Forest *)adaptIn->data : NULL;
4700: if (forestAdaptOut == forestIn) {
4701: switch (purposeOut) {
4702: case DM_ADAPT_REFINE:
4703: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4704: PetscCall(PetscSFSetUp(inSF));
4705: break;
4706: case DM_ADAPT_COARSEN:
4707: case DM_ADAPT_COARSEN_LAST:
4708: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &outCids));
4709: PetscCall(PetscSFSetUp(outSF));
4710: break;
4711: default:
4712: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4713: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4714: PetscCall(PetscSFSetUp(inSF));
4715: PetscCall(PetscSFSetUp(outSF));
4716: }
4717: } else if (forestAdaptIn == forestOut) {
4718: switch (purposeIn) {
4719: case DM_ADAPT_REFINE:
4720: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &inCids));
4721: PetscCall(PetscSFSetUp(outSF));
4722: break;
4723: case DM_ADAPT_COARSEN:
4724: case DM_ADAPT_COARSEN_LAST:
4725: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4726: PetscCall(PetscSFSetUp(inSF));
4727: break;
4728: default:
4729: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4730: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4731: PetscCall(PetscSFSetUp(inSF));
4732: PetscCall(PetscSFSetUp(outSF));
4733: }
4734: } else SETERRQ(PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Only support transfer from pre-adaptivity to post-adaptivity right now");
4735: PetscCall(DMPforestGetPlex(dmIn, &plexIn));
4736: PetscCall(DMPforestGetPlex(dmOut, &plexOut));
4738: PetscCall(DMPlexTransferVecTree(plexIn, vecIn, plexOut, vecOut, inSF, outSF, inCids, outCids, useBCs, time));
4739: PetscCall(PetscFree(inCids));
4740: PetscCall(PetscFree(outCids));
4741: PetscCall(PetscSFDestroy(&inSF));
4742: PetscCall(PetscSFDestroy(&outSF));
4743: PetscCall(PetscFree(inCids));
4744: PetscCall(PetscFree(outCids));
4745: PetscFunctionReturn(PETSC_SUCCESS);
4746: }
4748: #define DMCreateCoordinateDM_pforest _append_pforest(DMCreateCoordinateDM)
4749: static PetscErrorCode DMCreateCoordinateDM_pforest(DM dm, DM *cdm)
4750: {
4751: DM plex;
4753: PetscFunctionBegin;
4755: PetscCall(DMPforestGetPlex(dm, &plex));
4756: PetscCall(DMGetCoordinateDM(plex, cdm));
4757: PetscCall(PetscObjectReference((PetscObject)*cdm));
4758: PetscFunctionReturn(PETSC_SUCCESS);
4759: }
4761: #define VecViewLocal_pforest _append_pforest(VecViewLocal)
4762: static PetscErrorCode VecViewLocal_pforest(Vec vec, PetscViewer viewer)
4763: {
4764: DM dm, plex;
4766: PetscFunctionBegin;
4767: PetscCall(VecGetDM(vec, &dm));
4768: PetscCall(DMPforestGetPlex(dm, &plex));
4769: PetscCall(VecSetDM(vec, plex));
4770: PetscCall(VecView_Plex_Local(vec, viewer));
4771: PetscCall(VecSetDM(vec, dm));
4772: PetscFunctionReturn(PETSC_SUCCESS);
4773: }
4775: #define VecView_pforest _append_pforest(VecView)
4776: static PetscErrorCode VecView_pforest(Vec vec, PetscViewer viewer)
4777: {
4778: DM dm, plex;
4780: PetscFunctionBegin;
4781: PetscCall(VecGetDM(vec, &dm));
4782: PetscCall(DMPforestGetPlex(dm, &plex));
4783: PetscCall(VecSetDM(vec, plex));
4784: PetscCall(VecView_Plex(vec, viewer));
4785: PetscCall(VecSetDM(vec, dm));
4786: PetscFunctionReturn(PETSC_SUCCESS);
4787: }
4789: #define VecView_pforest_Native _infix_pforest(VecView, _Native)
4790: static PetscErrorCode VecView_pforest_Native(Vec vec, PetscViewer viewer)
4791: {
4792: DM dm, plex;
4794: PetscFunctionBegin;
4795: PetscCall(VecGetDM(vec, &dm));
4796: PetscCall(DMPforestGetPlex(dm, &plex));
4797: PetscCall(VecSetDM(vec, plex));
4798: PetscCall(VecView_Plex_Native(vec, viewer));
4799: PetscCall(VecSetDM(vec, dm));
4800: PetscFunctionReturn(PETSC_SUCCESS);
4801: }
4803: #define VecLoad_pforest _append_pforest(VecLoad)
4804: static PetscErrorCode VecLoad_pforest(Vec vec, PetscViewer viewer)
4805: {
4806: DM dm, plex;
4808: PetscFunctionBegin;
4809: PetscCall(VecGetDM(vec, &dm));
4810: PetscCall(DMPforestGetPlex(dm, &plex));
4811: PetscCall(VecSetDM(vec, plex));
4812: PetscCall(VecLoad_Plex(vec, viewer));
4813: PetscCall(VecSetDM(vec, dm));
4814: PetscFunctionReturn(PETSC_SUCCESS);
4815: }
4817: #define VecLoad_pforest_Native _infix_pforest(VecLoad, _Native)
4818: static PetscErrorCode VecLoad_pforest_Native(Vec vec, PetscViewer viewer)
4819: {
4820: DM dm, plex;
4822: PetscFunctionBegin;
4823: PetscCall(VecGetDM(vec, &dm));
4824: PetscCall(DMPforestGetPlex(dm, &plex));
4825: PetscCall(VecSetDM(vec, plex));
4826: PetscCall(VecLoad_Plex_Native(vec, viewer));
4827: PetscCall(VecSetDM(vec, dm));
4828: PetscFunctionReturn(PETSC_SUCCESS);
4829: }
4831: #define DMCreateGlobalVector_pforest _append_pforest(DMCreateGlobalVector)
4832: static PetscErrorCode DMCreateGlobalVector_pforest(DM dm, Vec *vec)
4833: {
4834: PetscFunctionBegin;
4835: PetscCall(DMCreateGlobalVector_Section_Private(dm, vec));
4836: /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */
4837: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_pforest));
4838: PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void))VecView_pforest_Native));
4839: PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_pforest));
4840: PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void))VecLoad_pforest_Native));
4841: PetscFunctionReturn(PETSC_SUCCESS);
4842: }
4844: #define DMCreateLocalVector_pforest _append_pforest(DMCreateLocalVector)
4845: static PetscErrorCode DMCreateLocalVector_pforest(DM dm, Vec *vec)
4846: {
4847: PetscFunctionBegin;
4848: PetscCall(DMCreateLocalVector_Section_Private(dm, vec));
4849: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecViewLocal_pforest));
4850: PetscFunctionReturn(PETSC_SUCCESS);
4851: }
4853: #define DMCreateMatrix_pforest _append_pforest(DMCreateMatrix)
4854: static PetscErrorCode DMCreateMatrix_pforest(DM dm, Mat *mat)
4855: {
4856: DM plex;
4858: PetscFunctionBegin;
4860: PetscCall(DMPforestGetPlex(dm, &plex));
4861: if (plex->prealloc_only != dm->prealloc_only) plex->prealloc_only = dm->prealloc_only; /* maybe this should go into forest->plex */
4862: PetscCall(DMCreateMatrix(plex, mat));
4863: PetscCall(MatSetDM(*mat, dm));
4864: PetscFunctionReturn(PETSC_SUCCESS);
4865: }
4867: #define DMProjectFunctionLocal_pforest _append_pforest(DMProjectFunctionLocal)
4868: static PetscErrorCode DMProjectFunctionLocal_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4869: {
4870: DM plex;
4872: PetscFunctionBegin;
4874: PetscCall(DMPforestGetPlex(dm, &plex));
4875: PetscCall(DMProjectFunctionLocal(plex, time, funcs, ctxs, mode, localX));
4876: PetscFunctionReturn(PETSC_SUCCESS);
4877: }
4879: #define DMProjectFunctionLabelLocal_pforest _append_pforest(DMProjectFunctionLabelLocal)
4880: 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)
4881: {
4882: DM plex;
4884: PetscFunctionBegin;
4886: PetscCall(DMPforestGetPlex(dm, &plex));
4887: PetscCall(DMProjectFunctionLabelLocal(plex, time, label, numIds, ids, Ncc, comps, funcs, ctxs, mode, localX));
4888: PetscFunctionReturn(PETSC_SUCCESS);
4889: }
4891: #define DMProjectFieldLocal_pforest _append_pforest(DMProjectFieldLocal)
4892: 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)
4893: {
4894: DM plex;
4896: PetscFunctionBegin;
4898: PetscCall(DMPforestGetPlex(dm, &plex));
4899: PetscCall(DMProjectFieldLocal(plex, time, localU, funcs, mode, localX));
4900: PetscFunctionReturn(PETSC_SUCCESS);
4901: }
4903: #define DMComputeL2Diff_pforest _append_pforest(DMComputeL2Diff)
4904: PetscErrorCode DMComputeL2Diff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
4905: {
4906: DM plex;
4908: PetscFunctionBegin;
4910: PetscCall(DMPforestGetPlex(dm, &plex));
4911: PetscCall(DMComputeL2Diff(plex, time, funcs, ctxs, X, diff));
4912: PetscFunctionReturn(PETSC_SUCCESS);
4913: }
4915: #define DMComputeL2FieldDiff_pforest _append_pforest(DMComputeL2FieldDiff)
4916: PetscErrorCode DMComputeL2FieldDiff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
4917: {
4918: DM plex;
4920: PetscFunctionBegin;
4922: PetscCall(DMPforestGetPlex(dm, &plex));
4923: PetscCall(DMComputeL2FieldDiff(plex, time, funcs, ctxs, X, diff));
4924: PetscFunctionReturn(PETSC_SUCCESS);
4925: }
4927: #define DMCreatelocalsection_pforest _append_pforest(DMCreatelocalsection)
4928: static PetscErrorCode DMCreatelocalsection_pforest(DM dm)
4929: {
4930: DM plex;
4931: PetscSection section;
4933: PetscFunctionBegin;
4935: PetscCall(DMPforestGetPlex(dm, &plex));
4936: PetscCall(DMGetLocalSection(plex, §ion));
4937: PetscCall(DMSetLocalSection(dm, section));
4938: PetscFunctionReturn(PETSC_SUCCESS);
4939: }
4941: #define DMCreateDefaultConstraints_pforest _append_pforest(DMCreateDefaultConstraints)
4942: static PetscErrorCode DMCreateDefaultConstraints_pforest(DM dm)
4943: {
4944: DM plex;
4945: Mat mat;
4946: Vec bias;
4947: PetscSection section;
4949: PetscFunctionBegin;
4951: PetscCall(DMPforestGetPlex(dm, &plex));
4952: PetscCall(DMGetDefaultConstraints(plex, §ion, &mat, &bias));
4953: PetscCall(DMSetDefaultConstraints(dm, section, mat, bias));
4954: PetscFunctionReturn(PETSC_SUCCESS);
4955: }
4957: #define DMGetDimPoints_pforest _append_pforest(DMGetDimPoints)
4958: static PetscErrorCode DMGetDimPoints_pforest(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
4959: {
4960: DM plex;
4962: PetscFunctionBegin;
4964: PetscCall(DMPforestGetPlex(dm, &plex));
4965: PetscCall(DMGetDimPoints(plex, dim, cStart, cEnd));
4966: PetscFunctionReturn(PETSC_SUCCESS);
4967: }
4969: /* Need to forward declare */
4970: #define DMInitialize_pforest _append_pforest(DMInitialize)
4971: static PetscErrorCode DMInitialize_pforest(DM dm);
4973: #define DMClone_pforest _append_pforest(DMClone)
4974: static PetscErrorCode DMClone_pforest(DM dm, DM *newdm)
4975: {
4976: PetscFunctionBegin;
4977: PetscCall(DMClone_Forest(dm, newdm));
4978: PetscCall(DMInitialize_pforest(*newdm));
4979: PetscFunctionReturn(PETSC_SUCCESS);
4980: }
4982: #define DMForestCreateCellChart_pforest _append_pforest(DMForestCreateCellChart)
4983: static PetscErrorCode DMForestCreateCellChart_pforest(DM dm, PetscInt *cStart, PetscInt *cEnd)
4984: {
4985: DM_Forest *forest;
4986: DM_Forest_pforest *pforest;
4987: PetscInt overlap;
4989: PetscFunctionBegin;
4990: PetscCall(DMSetUp(dm));
4991: forest = (DM_Forest *)dm->data;
4992: pforest = (DM_Forest_pforest *)forest->data;
4993: *cStart = 0;
4994: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
4995: if (overlap && pforest->ghost) {
4996: *cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[pforest->forest->mpisize];
4997: } else {
4998: *cEnd = pforest->forest->local_num_quadrants;
4999: }
5000: PetscFunctionReturn(PETSC_SUCCESS);
5001: }
5003: #define DMForestCreateCellSF_pforest _append_pforest(DMForestCreateCellSF)
5004: static PetscErrorCode DMForestCreateCellSF_pforest(DM dm, PetscSF *cellSF)
5005: {
5006: DM_Forest *forest;
5007: DM_Forest_pforest *pforest;
5008: PetscMPIInt rank;
5009: PetscInt overlap;
5010: PetscInt cStart, cEnd, cLocalStart, cLocalEnd;
5011: PetscInt nRoots, nLeaves, *mine = NULL;
5012: PetscSFNode *remote = NULL;
5013: PetscSF sf;
5015: PetscFunctionBegin;
5016: PetscCall(DMForestGetCellChart(dm, &cStart, &cEnd));
5017: forest = (DM_Forest *)dm->data;
5018: pforest = (DM_Forest_pforest *)forest->data;
5019: nRoots = cEnd - cStart;
5020: cLocalStart = pforest->cLocalStart;
5021: cLocalEnd = pforest->cLocalEnd;
5022: nLeaves = 0;
5023: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
5024: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
5025: if (overlap && pforest->ghost) {
5026: PetscSFNode *mirror;
5027: p4est_quadrant_t *mirror_array;
5028: PetscInt nMirror, nGhostPre, nSelf, q;
5029: void **mirrorPtrs;
5031: nMirror = (PetscInt)pforest->ghost->mirrors.elem_count;
5032: nSelf = cLocalEnd - cLocalStart;
5033: nLeaves = nRoots - nSelf;
5034: nGhostPre = (PetscInt)pforest->ghost->proc_offsets[rank];
5035: PetscCall(PetscMalloc1(nLeaves, &mine));
5036: PetscCall(PetscMalloc1(nLeaves, &remote));
5037: PetscCall(PetscMalloc2(nMirror, &mirror, nMirror, &mirrorPtrs));
5038: mirror_array = (p4est_quadrant_t *)pforest->ghost->mirrors.array;
5039: for (q = 0; q < nMirror; q++) {
5040: p4est_quadrant_t *mir = &(mirror_array[q]);
5042: mirror[q].rank = rank;
5043: mirror[q].index = (PetscInt)mir->p.piggy3.local_num + cLocalStart;
5044: mirrorPtrs[q] = (void *)&(mirror[q]);
5045: }
5046: PetscCallP4est(p4est_ghost_exchange_custom, (pforest->forest, pforest->ghost, sizeof(PetscSFNode), mirrorPtrs, remote));
5047: PetscCall(PetscFree2(mirror, mirrorPtrs));
5048: for (q = 0; q < nGhostPre; q++) mine[q] = q;
5049: for (; q < nLeaves; q++) mine[q] = (q - nGhostPre) + cLocalEnd;
5050: }
5051: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf));
5052: PetscCall(PetscSFSetGraph(sf, nRoots, nLeaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
5053: *cellSF = sf;
5054: PetscFunctionReturn(PETSC_SUCCESS);
5055: }
5057: static PetscErrorCode DMCreateNeumannOverlap_pforest(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx)
5058: {
5059: DM plex;
5061: PetscFunctionBegin;
5062: PetscCall(DMPforestGetPlex(dm, &plex));
5063: PetscCall(DMCreateNeumannOverlap_Plex(plex, ovl, J, setup, setup_ctx));
5064: if (!*setup) {
5065: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup));
5066: if (*setup) PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm));
5067: }
5068: PetscFunctionReturn(PETSC_SUCCESS);
5069: }
5071: static PetscErrorCode DMInitialize_pforest(DM dm)
5072: {
5073: PetscFunctionBegin;
5074: dm->ops->setup = DMSetUp_pforest;
5075: dm->ops->view = DMView_pforest;
5076: dm->ops->clone = DMClone_pforest;
5077: dm->ops->createinterpolation = DMCreateInterpolation_pforest;
5078: dm->ops->createinjection = DMCreateInjection_pforest;
5079: dm->ops->setfromoptions = DMSetFromOptions_pforest;
5080: dm->ops->createcoordinatedm = DMCreateCoordinateDM_pforest;
5081: dm->ops->createglobalvector = DMCreateGlobalVector_pforest;
5082: dm->ops->createlocalvector = DMCreateLocalVector_pforest;
5083: dm->ops->creatematrix = DMCreateMatrix_pforest;
5084: dm->ops->projectfunctionlocal = DMProjectFunctionLocal_pforest;
5085: dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_pforest;
5086: dm->ops->projectfieldlocal = DMProjectFieldLocal_pforest;
5087: dm->ops->createlocalsection = DMCreatelocalsection_pforest;
5088: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_pforest;
5089: dm->ops->computel2diff = DMComputeL2Diff_pforest;
5090: dm->ops->computel2fielddiff = DMComputeL2FieldDiff_pforest;
5091: dm->ops->getdimpoints = DMGetDimPoints_pforest;
5093: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", DMConvert_plex_pforest));
5094: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", DMConvert_pforest_plex));
5095: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_pforest));
5096: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMForestGetPartitionOverlap));
5097: PetscFunctionReturn(PETSC_SUCCESS);
5098: }
5100: #define DMCreate_pforest _append_pforest(DMCreate)
5101: PETSC_EXTERN PetscErrorCode DMCreate_pforest(DM dm)
5102: {
5103: DM_Forest *forest;
5104: DM_Forest_pforest *pforest;
5106: PetscFunctionBegin;
5107: PetscCall(PetscP4estInitialize());
5108: PetscCall(DMCreate_Forest(dm));
5109: PetscCall(DMInitialize_pforest(dm));
5110: PetscCall(DMSetDimension(dm, P4EST_DIM));
5112: /* set forest defaults */
5113: PetscCall(DMForestSetTopology(dm, "unit"));
5114: PetscCall(DMForestSetMinimumRefinement(dm, 0));
5115: PetscCall(DMForestSetInitialRefinement(dm, 0));
5116: PetscCall(DMForestSetMaximumRefinement(dm, P4EST_QMAXLEVEL));
5117: PetscCall(DMForestSetGradeFactor(dm, 2));
5118: PetscCall(DMForestSetAdjacencyDimension(dm, 0));
5119: PetscCall(DMForestSetPartitionOverlap(dm, 0));
5121: /* create p4est data */
5122: PetscCall(PetscNew(&pforest));
5124: forest = (DM_Forest *)dm->data;
5125: forest->data = pforest;
5126: forest->destroy = DMForestDestroy_pforest;
5127: forest->ftemplate = DMForestTemplate_pforest;
5128: forest->transfervec = DMForestTransferVec_pforest;
5129: forest->transfervecfrombase = DMForestTransferVecFromBase_pforest;
5130: forest->createcellchart = DMForestCreateCellChart_pforest;
5131: forest->createcellsf = DMForestCreateCellSF_pforest;
5132: forest->clearadaptivityforest = DMForestClearAdaptivityForest_pforest;
5133: forest->getadaptivitysuccess = DMForestGetAdaptivitySuccess_pforest;
5134: pforest->topo = NULL;
5135: pforest->forest = NULL;
5136: pforest->ghost = NULL;
5137: pforest->lnodes = NULL;
5138: pforest->partition_for_coarsening = PETSC_TRUE;
5139: pforest->coarsen_hierarchy = PETSC_FALSE;
5140: pforest->cLocalStart = -1;
5141: pforest->cLocalEnd = -1;
5142: pforest->labelsFinalized = PETSC_FALSE;
5143: pforest->ghostName = NULL;
5144: PetscFunctionReturn(PETSC_SUCCESS);
5145: }
5147: #endif /* defined(PETSC_HAVE_P4EST) */