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