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