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 defined(PETSC_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: PetscInt maxlevelloc = 0;
594: p4est_t *p4est;
596: PetscFunctionBegin;
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: maxlevelloc = PetscMax((PetscInt)tree->maxlevel, maxlevelloc);
605: }
606: PetscCallMPI(MPIU_Allreduce(&maxlevelloc, 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 numValues, 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, &numValues));
873: PetscCallMPI(MPIU_Allreduce(&numValues, &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 defined(PETSC_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, MPIU_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 !defined(PETSC_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 defined(PETSC_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: }
2154: if (dmB->sectionSF != dmA->sectionSF) {
2155: PetscCall(PetscObjectReference((PetscObject)dmA->sectionSF));
2156: PetscCall(PetscSFDestroy(&dmB->sectionSF));
2157: dmB->sectionSF = dmA->sectionSF;
2158: }
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /* Get an SF that broadcasts a coarse-cell covering of the local fine cells */
2163: static PetscErrorCode DMPforestGetCellCoveringSF(MPI_Comm comm, p4est_t *p4estC, p4est_t *p4estF, PetscInt cStart, PetscInt cEnd, PetscSF *coveringSF)
2164: {
2165: PetscInt startF, endF, startC, endC, p, nLeaves;
2166: PetscSFNode *leaves;
2167: PetscSF sf;
2168: PetscInt *recv, *send;
2169: PetscMPIInt tag;
2170: MPI_Request *recvReqs, *sendReqs;
2171: PetscSection section;
2173: PetscFunctionBegin;
2174: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estF, p4estC, &startC, &endC));
2175: PetscCall(PetscMalloc2(2 * (endC - startC), &recv, endC - startC, &recvReqs));
2176: PetscCall(PetscCommGetNewTag(comm, &tag));
2177: for (p = startC; p < endC; p++) {
2178: recvReqs[p - startC] = MPI_REQUEST_NULL; /* just in case we don't initiate a receive */
2179: if (p4estC->global_first_quadrant[p] == p4estC->global_first_quadrant[p + 1]) { /* empty coarse partition */
2180: recv[2 * (p - startC)] = 0;
2181: recv[2 * (p - startC) + 1] = 0;
2182: continue;
2183: }
2185: PetscCallMPI(MPIU_Irecv(&recv[2 * (p - startC)], 2, MPIU_INT, p, tag, comm, &recvReqs[p - startC]));
2186: }
2187: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estC, p4estF, &startF, &endF));
2188: PetscCall(PetscMalloc2(2 * (endF - startF), &send, endF - startF, &sendReqs));
2189: /* count the quadrants rank will send to each of [startF,endF) */
2190: for (p = startF; p < endF; p++) {
2191: p4est_quadrant_t *myFineStart = &p4estF->global_first_position[p];
2192: p4est_quadrant_t *myFineEnd = &p4estF->global_first_position[p + 1];
2193: PetscInt tStart = (PetscInt)myFineStart->p.which_tree;
2194: PetscInt tEnd = (PetscInt)myFineEnd->p.which_tree;
2195: PetscInt firstCell = -1, lastCell = -1;
2196: p4est_tree_t *treeStart = &(((p4est_tree_t *)p4estC->trees->array)[tStart]);
2197: p4est_tree_t *treeEnd = (size_t)tEnd < p4estC->trees->elem_count ? &(((p4est_tree_t *)p4estC->trees->array)[tEnd]) : NULL;
2198: ssize_t overlapIndex;
2200: sendReqs[p - startF] = MPI_REQUEST_NULL; /* just in case we don't initiate a send */
2201: if (p4estF->global_first_quadrant[p] == p4estF->global_first_quadrant[p + 1]) continue;
2203: /* locate myFineStart in (or before) a cell */
2204: if (treeStart->quadrants.elem_count) {
2205: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeStart->quadrants, myFineStart, p4est_quadrant_disjoint));
2206: if (overlapIndex < 0) {
2207: firstCell = 0;
2208: } else {
2209: firstCell = (PetscInt)(treeStart->quadrants_offset + overlapIndex);
2210: }
2211: } else {
2212: firstCell = 0;
2213: }
2214: if (treeEnd && treeEnd->quadrants.elem_count) {
2215: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeEnd->quadrants, myFineEnd, p4est_quadrant_disjoint));
2216: if (overlapIndex < 0) { /* all of this local section is overlapped */
2217: lastCell = p4estC->local_num_quadrants;
2218: } else {
2219: p4est_quadrant_t *container = &(((p4est_quadrant_t *)treeEnd->quadrants.array)[overlapIndex]);
2220: p4est_quadrant_t first_desc;
2221: int equal;
2223: PetscCallP4est(p4est_quadrant_first_descendant, (container, &first_desc, P4EST_QMAXLEVEL));
2224: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (myFineEnd, &first_desc));
2225: if (equal) {
2226: lastCell = (PetscInt)(treeEnd->quadrants_offset + overlapIndex);
2227: } else {
2228: lastCell = (PetscInt)(treeEnd->quadrants_offset + overlapIndex + 1);
2229: }
2230: }
2231: } else {
2232: lastCell = p4estC->local_num_quadrants;
2233: }
2234: send[2 * (p - startF)] = firstCell;
2235: send[2 * (p - startF) + 1] = lastCell - firstCell;
2236: PetscCallMPI(MPIU_Isend(&send[2 * (p - startF)], 2, MPIU_INT, p, tag, comm, &sendReqs[p - startF]));
2237: }
2238: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endC - startC), recvReqs, MPI_STATUSES_IGNORE));
2239: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
2240: PetscCall(PetscSectionSetChart(section, startC, endC));
2241: for (p = startC; p < endC; p++) {
2242: PetscInt numCells = recv[2 * (p - startC) + 1];
2243: PetscCall(PetscSectionSetDof(section, p, numCells));
2244: }
2245: PetscCall(PetscSectionSetUp(section));
2246: PetscCall(PetscSectionGetStorageSize(section, &nLeaves));
2247: PetscCall(PetscMalloc1(nLeaves, &leaves));
2248: for (p = startC; p < endC; p++) {
2249: PetscInt firstCell = recv[2 * (p - startC)];
2250: PetscInt numCells = recv[2 * (p - startC) + 1];
2251: PetscInt off, i;
2253: PetscCall(PetscSectionGetOffset(section, p, &off));
2254: for (i = 0; i < numCells; i++) {
2255: leaves[off + i].rank = p;
2256: leaves[off + i].index = firstCell + i;
2257: }
2258: }
2259: PetscCall(PetscSFCreate(comm, &sf));
2260: PetscCall(PetscSFSetGraph(sf, cEnd - cStart, nLeaves, NULL, PETSC_OWN_POINTER, leaves, PETSC_OWN_POINTER));
2261: PetscCall(PetscSectionDestroy(§ion));
2262: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endF - startF), sendReqs, MPI_STATUSES_IGNORE));
2263: PetscCall(PetscFree2(send, sendReqs));
2264: PetscCall(PetscFree2(recv, recvReqs));
2265: *coveringSF = sf;
2266: PetscFunctionReturn(PETSC_SUCCESS);
2267: }
2269: /* closure points for locally-owned cells */
2270: static PetscErrorCode DMPforestGetCellSFNodes(DM dm, PetscInt numClosureIndices, PetscInt *numClosurePoints, PetscSFNode **closurePoints, PetscBool redirect)
2271: {
2272: PetscInt cStart, cEnd;
2273: PetscInt count, c;
2274: PetscMPIInt rank;
2275: PetscInt closureSize = -1;
2276: PetscInt *closure = NULL;
2277: PetscSF pointSF;
2278: PetscInt nleaves, nroots;
2279: const PetscInt *ilocal;
2280: const PetscSFNode *iremote;
2281: DM plex;
2282: DM_Forest *forest;
2283: DM_Forest_pforest *pforest;
2285: PetscFunctionBegin;
2286: forest = (DM_Forest *)dm->data;
2287: pforest = (DM_Forest_pforest *)forest->data;
2288: cStart = pforest->cLocalStart;
2289: cEnd = pforest->cLocalEnd;
2290: PetscCall(DMPforestGetPlex(dm, &plex));
2291: PetscCall(DMGetPointSF(dm, &pointSF));
2292: PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
2293: nleaves = PetscMax(0, nleaves);
2294: nroots = PetscMax(0, nroots);
2295: *numClosurePoints = numClosureIndices * (cEnd - cStart);
2296: PetscCall(PetscMalloc1(*numClosurePoints, closurePoints));
2297: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2298: for (c = cStart, count = 0; c < cEnd; c++) {
2299: PetscInt i;
2300: PetscCall(DMPlexGetTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2302: for (i = 0; i < numClosureIndices; i++, count++) {
2303: PetscInt p = closure[2 * i];
2304: PetscInt loc = -1;
2306: PetscCall(PetscFindInt(p, nleaves, ilocal, &loc));
2307: if (redirect && loc >= 0) {
2308: (*closurePoints)[count].rank = iremote[loc].rank;
2309: (*closurePoints)[count].index = iremote[loc].index;
2310: } else {
2311: (*closurePoints)[count].rank = rank;
2312: (*closurePoints)[count].index = p;
2313: }
2314: }
2315: PetscCall(DMPlexRestoreTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2316: }
2317: PetscFunctionReturn(PETSC_SUCCESS);
2318: }
2320: static void MPIAPI DMPforestMaxSFNode(void *a, void *b, PetscMPIInt *len, MPI_Datatype *type)
2321: {
2322: PetscMPIInt i;
2324: for (i = 0; i < *len; i++) {
2325: PetscSFNode *A = (PetscSFNode *)a;
2326: PetscSFNode *B = (PetscSFNode *)b;
2328: if (B->rank < 0) *B = *A;
2329: }
2330: }
2332: static PetscErrorCode DMPforestGetTransferSF_Point(DM coarse, DM fine, PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2333: {
2334: MPI_Comm comm;
2335: PetscMPIInt rank, size;
2336: DM_Forest_pforest *pforestC, *pforestF;
2337: p4est_t *p4estC, *p4estF;
2338: PetscInt numClosureIndices;
2339: PetscInt numClosurePointsC, numClosurePointsF;
2340: PetscSFNode *closurePointsC, *closurePointsF;
2341: p4est_quadrant_t *coverQuads = NULL;
2342: p4est_quadrant_t **treeQuads;
2343: PetscInt *treeQuadCounts;
2344: MPI_Datatype nodeType;
2345: MPI_Datatype nodeClosureType;
2346: MPI_Op sfNodeReduce;
2347: p4est_topidx_t fltF, lltF, t;
2348: DM plexC, plexF;
2349: PetscInt pStartF, pEndF, pStartC, pEndC;
2350: PetscBool saveInCoarse = PETSC_FALSE;
2351: PetscBool saveInFine = PETSC_FALSE;
2352: PetscBool formCids = (childIds != NULL) ? PETSC_TRUE : PETSC_FALSE;
2353: PetscInt *cids = NULL;
2355: PetscFunctionBegin;
2356: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2357: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2358: p4estC = pforestC->forest;
2359: p4estF = pforestF->forest;
2360: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2361: comm = PetscObjectComm((PetscObject)coarse);
2362: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2363: PetscCallMPI(MPI_Comm_size(comm, &size));
2364: PetscCall(DMPforestGetPlex(fine, &plexF));
2365: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2366: PetscCall(DMPforestGetPlex(coarse, &plexC));
2367: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2368: { /* check if the results have been cached */
2369: DM adaptCoarse, adaptFine;
2371: PetscCall(DMForestGetAdaptivityForest(coarse, &adaptCoarse));
2372: PetscCall(DMForestGetAdaptivityForest(fine, &adaptFine));
2373: if (adaptCoarse && adaptCoarse->data == fine->data) { /* coarse is adapted from fine */
2374: if (pforestC->pointSelfToAdaptSF) {
2375: PetscCall(PetscObjectReference((PetscObject)pforestC->pointSelfToAdaptSF));
2376: *sf = pforestC->pointSelfToAdaptSF;
2377: if (childIds) {
2378: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2379: PetscCall(PetscArraycpy(cids, pforestC->pointSelfToAdaptCids, pEndF - pStartF));
2380: *childIds = cids;
2381: }
2382: PetscFunctionReturn(PETSC_SUCCESS);
2383: } else {
2384: saveInCoarse = PETSC_TRUE;
2385: formCids = PETSC_TRUE;
2386: }
2387: } else if (adaptFine && adaptFine->data == coarse->data) { /* fine is adapted from coarse */
2388: if (pforestF->pointAdaptToSelfSF) {
2389: PetscCall(PetscObjectReference((PetscObject)pforestF->pointAdaptToSelfSF));
2390: *sf = pforestF->pointAdaptToSelfSF;
2391: if (childIds) {
2392: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2393: PetscCall(PetscArraycpy(cids, pforestF->pointAdaptToSelfCids, pEndF - pStartF));
2394: *childIds = cids;
2395: }
2396: PetscFunctionReturn(PETSC_SUCCESS);
2397: } else {
2398: saveInFine = PETSC_TRUE;
2399: formCids = PETSC_TRUE;
2400: }
2401: }
2402: }
2404: /* count the number of closure points that have dofs and create a list */
2405: numClosureIndices = P4EST_INSUL;
2406: /* create the datatype */
2407: PetscCallMPI(MPI_Type_contiguous(2, MPIU_INT, &nodeType));
2408: PetscCallMPI(MPI_Type_commit(&nodeType));
2409: PetscCallMPI(MPI_Op_create(DMPforestMaxSFNode, PETSC_FALSE, &sfNodeReduce));
2410: PetscCallMPI(MPI_Type_contiguous(numClosureIndices * 2, MPIU_INT, &nodeClosureType));
2411: PetscCallMPI(MPI_Type_commit(&nodeClosureType));
2412: /* everything has to go through cells: for each cell, create a list of the sfnodes in its closure */
2413: /* get lists of closure point SF nodes for every cell */
2414: PetscCall(DMPforestGetCellSFNodes(coarse, numClosureIndices, &numClosurePointsC, &closurePointsC, PETSC_TRUE));
2415: PetscCall(DMPforestGetCellSFNodes(fine, numClosureIndices, &numClosurePointsF, &closurePointsF, PETSC_FALSE));
2416: /* create pointers for tree lists */
2417: fltF = p4estF->first_local_tree;
2418: lltF = p4estF->last_local_tree;
2419: PetscCall(PetscCalloc2(lltF + 1 - fltF, &treeQuads, lltF + 1 - fltF, &treeQuadCounts));
2420: /* if the partitions don't match, ship the coarse to cover the fine */
2421: if (size > 1) {
2422: PetscInt p;
2424: for (p = 0; p < size; p++) {
2425: int equal;
2427: PetscCallP4estReturn(equal, p4est_quadrant_is_equal_piggy, (&p4estC->global_first_position[p], &p4estF->global_first_position[p]));
2428: if (!equal) break;
2429: }
2430: if (p < size) { /* non-matching distribution: send the coarse to cover the fine */
2431: PetscInt cStartC, cEndC;
2432: PetscSF coveringSF;
2433: PetscInt nleaves;
2434: PetscInt count;
2435: PetscSFNode *newClosurePointsC;
2436: p4est_quadrant_t *coverQuadsSend;
2437: p4est_topidx_t fltC = p4estC->first_local_tree;
2438: p4est_topidx_t lltC = p4estC->last_local_tree;
2439: p4est_topidx_t t;
2440: PetscMPIInt blockSizes[4] = {P4EST_DIM, 2, 1, 1};
2441: MPI_Aint blockOffsets[4] = {offsetof(p4est_quadrant_t, x), offsetof(p4est_quadrant_t, level), offsetof(p4est_quadrant_t, pad16), offsetof(p4est_quadrant_t, p)};
2442: MPI_Datatype blockTypes[4] = {MPI_INT32_T, MPI_INT8_T, MPI_INT16_T, MPI_INT32_T /* p.which_tree */};
2443: MPI_Datatype quadStruct, quadType;
2445: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, 0, &cStartC, &cEndC));
2446: PetscCall(DMPforestGetCellCoveringSF(comm, p4estC, p4estF, pforestC->cLocalStart, pforestC->cLocalEnd, &coveringSF));
2447: PetscCall(PetscSFGetGraph(coveringSF, NULL, &nleaves, NULL, NULL));
2448: PetscCall(PetscMalloc1(numClosureIndices * nleaves, &newClosurePointsC));
2449: PetscCall(PetscMalloc1(nleaves, &coverQuads));
2450: PetscCall(PetscMalloc1(cEndC - cStartC, &coverQuadsSend));
2451: count = 0;
2452: for (t = fltC; t <= lltC; t++) { /* unfortunately, we need to pack a send array, since quads are not stored packed in p4est */
2453: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2454: PetscInt q;
2456: PetscCall(PetscMemcpy(&coverQuadsSend[count], tree->quadrants.array, tree->quadrants.elem_count * sizeof(p4est_quadrant_t)));
2457: for (q = 0; (size_t)q < tree->quadrants.elem_count; q++) coverQuadsSend[count + q].p.which_tree = t;
2458: count += tree->quadrants.elem_count;
2459: }
2460: /* p is of a union type p4est_quadrant_data, but only the p.which_tree field is active at this time. So, we
2461: have a simple blockTypes[] to use. Note that quadStruct does not count potential padding in array of
2462: p4est_quadrant_t. We have to call MPI_Type_create_resized() to change upper-bound of quadStruct.
2463: */
2464: PetscCallMPI(MPI_Type_create_struct(4, blockSizes, blockOffsets, blockTypes, &quadStruct));
2465: PetscCallMPI(MPI_Type_create_resized(quadStruct, 0, sizeof(p4est_quadrant_t), &quadType));
2466: PetscCallMPI(MPI_Type_commit(&quadType));
2467: PetscCall(PetscSFBcastBegin(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2468: PetscCall(PetscSFBcastBegin(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2469: PetscCall(PetscSFBcastEnd(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2470: PetscCall(PetscSFBcastEnd(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2471: PetscCallMPI(MPI_Type_free(&quadStruct));
2472: PetscCallMPI(MPI_Type_free(&quadType));
2473: PetscCall(PetscFree(coverQuadsSend));
2474: PetscCall(PetscFree(closurePointsC));
2475: PetscCall(PetscSFDestroy(&coveringSF));
2476: closurePointsC = newClosurePointsC;
2478: /* assign tree quads based on locations in coverQuads */
2479: {
2480: PetscInt q;
2481: for (q = 0; q < nleaves; q++) {
2482: p4est_locidx_t t = coverQuads[q].p.which_tree;
2483: if (!treeQuadCounts[t - fltF]++) treeQuads[t - fltF] = &coverQuads[q];
2484: }
2485: }
2486: }
2487: }
2488: if (!coverQuads) { /* matching partitions: assign tree quads based on locations in p4est native arrays */
2489: for (t = fltF; t <= lltF; t++) {
2490: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2492: treeQuadCounts[t - fltF] = (PetscInt)tree->quadrants.elem_count;
2493: treeQuads[t - fltF] = (p4est_quadrant_t *)tree->quadrants.array;
2494: }
2495: }
2497: {
2498: PetscInt p;
2499: PetscInt cLocalStartF;
2500: PetscSF pointSF;
2501: PetscSFNode *roots;
2502: PetscInt *rootType;
2503: DM refTree = NULL;
2504: DMLabel canonical;
2505: PetscInt *childClosures[P4EST_CHILDREN] = {NULL};
2506: PetscInt *rootClosure = NULL;
2507: PetscInt coarseOffset;
2508: PetscInt numCoarseQuads;
2510: PetscCall(PetscMalloc1(pEndF - pStartF, &roots));
2511: PetscCall(PetscMalloc1(pEndF - pStartF, &rootType));
2512: PetscCall(DMGetPointSF(fine, &pointSF));
2513: for (p = pStartF; p < pEndF; p++) {
2514: roots[p - pStartF].rank = -1;
2515: roots[p - pStartF].index = -1;
2516: rootType[p - pStartF] = -1;
2517: }
2518: if (formCids) {
2519: PetscInt child;
2521: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2522: for (p = pStartF; p < pEndF; p++) cids[p - pStartF] = -2;
2523: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2524: PetscCall(DMPlexGetTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2525: for (child = 0; child < P4EST_CHILDREN; child++) { /* get the closures of the child cells in the reference tree */
2526: PetscCall(DMPlexGetTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2527: }
2528: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2529: }
2530: cLocalStartF = pforestF->cLocalStart;
2531: for (t = fltF, coarseOffset = 0, numCoarseQuads = 0; t <= lltF; t++, coarseOffset += numCoarseQuads) {
2532: p4est_tree_t *tree = &(((p4est_tree_t *)p4estF->trees->array)[t]);
2533: PetscInt numFineQuads = (PetscInt)tree->quadrants.elem_count;
2534: p4est_quadrant_t *coarseQuads = treeQuads[t - fltF];
2535: p4est_quadrant_t *fineQuads = (p4est_quadrant_t *)tree->quadrants.array;
2536: PetscInt i, coarseCount = 0;
2537: PetscInt offset = tree->quadrants_offset;
2538: sc_array_t coarseQuadsArray;
2540: numCoarseQuads = treeQuadCounts[t - fltF];
2541: PetscCallP4est(sc_array_init_data, (&coarseQuadsArray, coarseQuads, sizeof(p4est_quadrant_t), (size_t)numCoarseQuads));
2542: for (i = 0; i < numFineQuads; i++) {
2543: PetscInt c = i + offset;
2544: p4est_quadrant_t *quad = &fineQuads[i];
2545: p4est_quadrant_t *quadCoarse = NULL;
2546: ssize_t disjoint = -1;
2548: while (disjoint < 0 && coarseCount < numCoarseQuads) {
2549: quadCoarse = &coarseQuads[coarseCount];
2550: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2551: if (disjoint < 0) coarseCount++;
2552: }
2553: PetscCheck(disjoint == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "did not find overlapping coarse quad");
2554: if (quadCoarse->level > quad->level || (quadCoarse->level == quad->level && !transferIdent)) { /* the "coarse" mesh is finer than the fine mesh at the point: continue */
2555: if (transferIdent) { /* find corners */
2556: PetscInt j = 0;
2558: do {
2559: if (j < P4EST_CHILDREN) {
2560: p4est_quadrant_t cornerQuad;
2561: int equal;
2563: PetscCallP4est(p4est_quadrant_corner_descendant, (quad, &cornerQuad, j, quadCoarse->level));
2564: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (&cornerQuad, quadCoarse));
2565: if (equal) {
2566: PetscInt petscJ = P4estVertToPetscVert[j];
2567: PetscInt p = closurePointsF[numClosureIndices * c + (P4EST_INSUL - P4EST_CHILDREN) + petscJ].index;
2568: PetscSFNode q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + (P4EST_INSUL - P4EST_CHILDREN) + petscJ];
2570: roots[p - pStartF] = q;
2571: rootType[p - pStartF] = PETSC_INT_MAX;
2572: cids[p - pStartF] = -1;
2573: j++;
2574: }
2575: }
2576: coarseCount++;
2577: disjoint = 1;
2578: if (coarseCount < numCoarseQuads) {
2579: quadCoarse = &coarseQuads[coarseCount];
2580: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2581: }
2582: } while (!disjoint);
2583: }
2584: continue;
2585: }
2586: if (quadCoarse->level == quad->level) { /* same quad present in coarse and fine mesh */
2587: PetscInt j;
2588: for (j = 0; j < numClosureIndices; j++) {
2589: PetscInt p = closurePointsF[numClosureIndices * c + j].index;
2591: roots[p - pStartF] = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + j];
2592: rootType[p - pStartF] = PETSC_INT_MAX; /* unconditionally accept */
2593: cids[p - pStartF] = -1;
2594: }
2595: } else {
2596: PetscInt levelDiff = quad->level - quadCoarse->level;
2597: PetscInt proposedCids[P4EST_INSUL] = {0};
2599: if (formCids) {
2600: PetscInt cl;
2601: PetscInt *pointClosure = NULL;
2602: int cid;
2604: PetscCheck(levelDiff <= 1, PETSC_COMM_SELF, PETSC_ERR_USER, "Recursive child ids not implemented");
2605: PetscCallP4estReturn(cid, p4est_quadrant_child_id, (quad));
2606: PetscCall(DMPlexGetTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2607: for (cl = 0; cl < P4EST_INSUL; cl++) {
2608: PetscInt p = pointClosure[2 * cl];
2609: PetscInt point = childClosures[cid][2 * cl];
2610: PetscInt ornt = childClosures[cid][2 * cl + 1];
2611: PetscInt newcid = -1;
2612: DMPolytopeType ct;
2614: if (rootType[p - pStartF] == PETSC_INT_MAX) continue;
2615: PetscCall(DMPlexGetCellType(refTree, point, &ct));
2616: ornt = DMPolytopeConvertNewOrientation_Internal(ct, ornt);
2617: if (!cl) {
2618: newcid = cid + 1;
2619: } else {
2620: PetscInt rcl, parent, parentOrnt = 0;
2622: PetscCall(DMPlexGetTreeParent(refTree, point, &parent, NULL));
2623: if (parent == point) {
2624: newcid = -1;
2625: } else if (!parent) { /* in the root */
2626: newcid = point;
2627: } else {
2628: DMPolytopeType rct = DM_POLYTOPE_UNKNOWN;
2630: for (rcl = 1; rcl < P4EST_INSUL; rcl++) {
2631: if (rootClosure[2 * rcl] == parent) {
2632: PetscCall(DMPlexGetCellType(refTree, parent, &rct));
2633: parentOrnt = DMPolytopeConvertNewOrientation_Internal(rct, rootClosure[2 * rcl + 1]);
2634: break;
2635: }
2636: }
2637: PetscCheck(rcl < P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Couldn't find parent in root closure");
2638: PetscCall(DMPlexReferenceTreeGetChildSymmetry(refTree, parent, parentOrnt, ornt, point, DMPolytopeConvertNewOrientation_Internal(rct, pointClosure[2 * rcl + 1]), NULL, &newcid));
2639: }
2640: }
2641: if (newcid >= 0) {
2642: if (canonical) PetscCall(DMLabelGetValue(canonical, newcid, &newcid));
2643: proposedCids[cl] = newcid;
2644: }
2645: }
2646: PetscCall(DMPlexRestoreTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2647: }
2648: p4est_qcoord_t coarseBound[2][P4EST_DIM] = {
2649: {quadCoarse->x, quadCoarse->y,
2650: #if defined(P4_TO_P8)
2651: quadCoarse->z
2652: #endif
2653: },
2654: {0}
2655: };
2656: p4est_qcoord_t fineBound[2][P4EST_DIM] = {
2657: {quad->x, quad->y,
2658: #if defined(P4_TO_P8)
2659: quad->z
2660: #endif
2661: },
2662: {0}
2663: };
2664: PetscInt j;
2665: for (j = 0; j < P4EST_DIM; j++) { /* get the coordinates of cell boundaries in each direction */
2666: coarseBound[1][j] = coarseBound[0][j] + P4EST_QUADRANT_LEN(quadCoarse->level);
2667: fineBound[1][j] = fineBound[0][j] + P4EST_QUADRANT_LEN(quad->level);
2668: }
2669: for (j = 0; j < numClosureIndices; j++) {
2670: PetscInt l, p;
2671: PetscSFNode q;
2673: p = closurePointsF[numClosureIndices * c + j].index;
2674: if (rootType[p - pStartF] == PETSC_INT_MAX) continue;
2675: if (j == 0) { /* volume: ancestor is volume */
2676: l = 0;
2677: } else if (j < 1 + P4EST_FACES) { /* facet */
2678: PetscInt face = PetscFaceToP4estFace[j - 1];
2679: PetscInt direction = face / 2;
2680: PetscInt coarseFace = -1;
2682: if (coarseBound[face % 2][direction] == fineBound[face % 2][direction]) {
2683: coarseFace = face;
2684: l = 1 + P4estFaceToPetscFace[coarseFace];
2685: } else {
2686: l = 0;
2687: }
2688: #if defined(P4_TO_P8)
2689: } else if (j < 1 + P4EST_FACES + P8EST_EDGES) {
2690: PetscInt edge = PetscEdgeToP4estEdge[j - (1 + P4EST_FACES)];
2691: PetscInt direction = edge / 4;
2692: PetscInt mod = edge % 4;
2693: PetscInt coarseEdge = -1, coarseFace = -1;
2694: PetscInt minDir = PetscMin((direction + 1) % 3, (direction + 2) % 3);
2695: PetscInt maxDir = PetscMax((direction + 1) % 3, (direction + 2) % 3);
2696: PetscBool dirTest[2];
2698: dirTest[0] = (PetscBool)(coarseBound[mod % 2][minDir] == fineBound[mod % 2][minDir]);
2699: dirTest[1] = (PetscBool)(coarseBound[mod / 2][maxDir] == fineBound[mod / 2][maxDir]);
2701: if (dirTest[0] && dirTest[1]) { /* fine edge falls on coarse edge */
2702: coarseEdge = edge;
2703: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2704: } else if (dirTest[0]) { /* fine edge falls on a coarse face in the minDir direction */
2705: coarseFace = 2 * minDir + (mod % 2);
2706: l = 1 + P4estFaceToPetscFace[coarseFace];
2707: } else if (dirTest[1]) { /* fine edge falls on a coarse face in the maxDir direction */
2708: coarseFace = 2 * maxDir + (mod / 2);
2709: l = 1 + P4estFaceToPetscFace[coarseFace];
2710: } else {
2711: l = 0;
2712: }
2713: #endif
2714: } else {
2715: PetscInt vertex = PetscVertToP4estVert[P4EST_CHILDREN - (P4EST_INSUL - j)];
2716: PetscBool dirTest[P4EST_DIM];
2717: PetscInt m;
2718: PetscInt numMatch = 0;
2719: PetscInt coarseVertex = -1, coarseFace = -1;
2720: #if defined(P4_TO_P8)
2721: PetscInt coarseEdge = -1;
2722: #endif
2724: for (m = 0; m < P4EST_DIM; m++) {
2725: dirTest[m] = (PetscBool)(coarseBound[(vertex >> m) & 1][m] == fineBound[(vertex >> m) & 1][m]);
2726: if (dirTest[m]) numMatch++;
2727: }
2728: if (numMatch == P4EST_DIM) { /* vertex on vertex */
2729: coarseVertex = vertex;
2730: l = P4EST_INSUL - (P4EST_CHILDREN - P4estVertToPetscVert[coarseVertex]);
2731: } else if (numMatch == 1) { /* vertex on face */
2732: for (m = 0; m < P4EST_DIM; m++) {
2733: if (dirTest[m]) {
2734: coarseFace = 2 * m + ((vertex >> m) & 1);
2735: break;
2736: }
2737: }
2738: l = 1 + P4estFaceToPetscFace[coarseFace];
2739: #if defined(P4_TO_P8)
2740: } else if (numMatch == 2) { /* vertex on edge */
2741: for (m = 0; m < P4EST_DIM; m++) {
2742: if (!dirTest[m]) {
2743: PetscInt otherDir1 = (m + 1) % 3;
2744: PetscInt otherDir2 = (m + 2) % 3;
2745: PetscInt minDir = PetscMin(otherDir1, otherDir2);
2746: PetscInt maxDir = PetscMax(otherDir1, otherDir2);
2748: coarseEdge = m * 4 + 2 * ((vertex >> maxDir) & 1) + ((vertex >> minDir) & 1);
2749: break;
2750: }
2751: }
2752: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2753: #endif
2754: } else { /* volume */
2755: l = 0;
2756: }
2757: }
2758: q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + l];
2759: if (l > rootType[p - pStartF]) {
2760: if (l >= P4EST_INSUL - P4EST_CHILDREN) { /* vertex on vertex: unconditional acceptance */
2761: if (transferIdent) {
2762: roots[p - pStartF] = q;
2763: rootType[p - pStartF] = PETSC_INT_MAX;
2764: if (formCids) cids[p - pStartF] = -1;
2765: }
2766: } else {
2767: PetscInt k, thisp = p, limit;
2769: roots[p - pStartF] = q;
2770: rootType[p - pStartF] = l;
2771: if (formCids) cids[p - pStartF] = proposedCids[j];
2772: limit = transferIdent ? levelDiff : (levelDiff - 1);
2773: for (k = 0; k < limit; k++) {
2774: PetscInt parent;
2776: PetscCall(DMPlexGetTreeParent(plexF, thisp, &parent, NULL));
2777: if (parent == thisp) break;
2779: roots[parent - pStartF] = q;
2780: rootType[parent - pStartF] = PETSC_INT_MAX;
2781: if (formCids) cids[parent - pStartF] = -1;
2782: thisp = parent;
2783: }
2784: }
2785: }
2786: }
2787: }
2788: }
2789: }
2791: /* 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 */
2792: if (size > 1) {
2793: PetscInt *rootTypeCopy, p;
2795: PetscCall(PetscMalloc1(pEndF - pStartF, &rootTypeCopy));
2796: PetscCall(PetscArraycpy(rootTypeCopy, rootType, pEndF - pStartF));
2797: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2798: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2799: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2800: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2801: for (p = pStartF; p < pEndF; p++) {
2802: 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 */
2803: roots[p - pStartF].rank = -1;
2804: roots[p - pStartF].index = -1;
2805: }
2806: if (formCids && rootTypeCopy[p - pStartF] == PETSC_INT_MAX) cids[p - pStartF] = -1; /* we have found an antecedent that is the same: no child id */
2807: }
2808: PetscCall(PetscFree(rootTypeCopy));
2809: PetscCall(PetscSFReduceBegin(pointSF, nodeType, roots, roots, sfNodeReduce));
2810: PetscCall(PetscSFReduceEnd(pointSF, nodeType, roots, roots, sfNodeReduce));
2811: PetscCall(PetscSFBcastBegin(pointSF, nodeType, roots, roots, MPI_REPLACE));
2812: PetscCall(PetscSFBcastEnd(pointSF, nodeType, roots, roots, MPI_REPLACE));
2813: }
2814: PetscCall(PetscFree(rootType));
2816: {
2817: PetscInt numRoots;
2818: PetscInt numLeaves;
2819: PetscInt *leaves;
2820: PetscSFNode *iremote;
2821: /* count leaves */
2823: numRoots = pEndC - pStartC;
2825: numLeaves = 0;
2826: for (p = pStartF; p < pEndF; p++) {
2827: if (roots[p - pStartF].index >= 0) numLeaves++;
2828: }
2829: PetscCall(PetscMalloc1(numLeaves, &leaves));
2830: PetscCall(PetscMalloc1(numLeaves, &iremote));
2831: numLeaves = 0;
2832: for (p = pStartF; p < pEndF; p++) {
2833: if (roots[p - pStartF].index >= 0) {
2834: leaves[numLeaves] = p - pStartF;
2835: iremote[numLeaves] = roots[p - pStartF];
2836: numLeaves++;
2837: }
2838: }
2839: PetscCall(PetscFree(roots));
2840: PetscCall(PetscSFCreate(comm, sf));
2841: if (numLeaves == (pEndF - pStartF)) {
2842: PetscCall(PetscFree(leaves));
2843: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2844: } else {
2845: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, leaves, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2846: }
2847: }
2848: if (formCids) {
2849: PetscSF pointSF;
2850: PetscInt child;
2852: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2853: PetscCall(DMGetPointSF(plexF, &pointSF));
2854: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2855: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2856: if (childIds) *childIds = cids;
2857: for (child = 0; child < P4EST_CHILDREN; child++) PetscCall(DMPlexRestoreTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2858: PetscCall(DMPlexRestoreTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2859: }
2860: }
2861: if (saveInCoarse) { /* cache results */
2862: PetscCall(PetscObjectReference((PetscObject)*sf));
2863: pforestC->pointSelfToAdaptSF = *sf;
2864: if (!childIds) {
2865: pforestC->pointSelfToAdaptCids = cids;
2866: } else {
2867: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestC->pointSelfToAdaptCids));
2868: PetscCall(PetscArraycpy(pforestC->pointSelfToAdaptCids, cids, pEndF - pStartF));
2869: }
2870: } else if (saveInFine) {
2871: PetscCall(PetscObjectReference((PetscObject)*sf));
2872: pforestF->pointAdaptToSelfSF = *sf;
2873: if (!childIds) {
2874: pforestF->pointAdaptToSelfCids = cids;
2875: } else {
2876: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestF->pointAdaptToSelfCids));
2877: PetscCall(PetscArraycpy(pforestF->pointAdaptToSelfCids, cids, pEndF - pStartF));
2878: }
2879: }
2880: PetscCall(PetscFree2(treeQuads, treeQuadCounts));
2881: PetscCall(PetscFree(coverQuads));
2882: PetscCall(PetscFree(closurePointsC));
2883: PetscCall(PetscFree(closurePointsF));
2884: PetscCallMPI(MPI_Type_free(&nodeClosureType));
2885: PetscCallMPI(MPI_Op_free(&sfNodeReduce));
2886: PetscCallMPI(MPI_Type_free(&nodeType));
2887: PetscFunctionReturn(PETSC_SUCCESS);
2888: }
2890: /* children are sf leaves of parents */
2891: static PetscErrorCode DMPforestGetTransferSF_Internal(DM coarse, DM fine, const PetscInt dofPerDim[], PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2892: {
2893: MPI_Comm comm;
2894: PetscMPIInt rank;
2895: DM_Forest_pforest *pforestC, *pforestF;
2896: DM plexC, plexF;
2897: PetscInt pStartC, pEndC, pStartF, pEndF;
2898: PetscSF pointTransferSF;
2899: PetscBool allOnes = PETSC_TRUE;
2901: PetscFunctionBegin;
2902: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2903: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2904: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2905: comm = PetscObjectComm((PetscObject)coarse);
2906: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2908: {
2909: PetscInt i;
2910: for (i = 0; i <= P4EST_DIM; i++) {
2911: if (dofPerDim[i] != 1) {
2912: allOnes = PETSC_FALSE;
2913: break;
2914: }
2915: }
2916: }
2917: PetscCall(DMPforestGetTransferSF_Point(coarse, fine, &pointTransferSF, transferIdent, childIds));
2918: if (allOnes) {
2919: *sf = pointTransferSF;
2920: PetscFunctionReturn(PETSC_SUCCESS);
2921: }
2923: PetscCall(DMPforestGetPlex(fine, &plexF));
2924: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2925: PetscCall(DMPforestGetPlex(coarse, &plexC));
2926: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2927: {
2928: PetscInt numRoots;
2929: PetscInt numLeaves;
2930: const PetscInt *leaves;
2931: const PetscSFNode *iremote;
2932: PetscInt d;
2933: PetscSection leafSection, rootSection;
2935: /* count leaves */
2936: PetscCall(PetscSFGetGraph(pointTransferSF, &numRoots, &numLeaves, &leaves, &iremote));
2937: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &rootSection));
2938: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &leafSection));
2939: PetscCall(PetscSectionSetChart(rootSection, pStartC, pEndC));
2940: PetscCall(PetscSectionSetChart(leafSection, pStartF, pEndF));
2942: for (d = 0; d <= P4EST_DIM; d++) {
2943: PetscInt startC, endC, e;
2945: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, P4EST_DIM - d, &startC, &endC));
2946: for (e = startC; e < endC; e++) PetscCall(PetscSectionSetDof(rootSection, e, dofPerDim[d]));
2947: }
2949: for (d = 0; d <= P4EST_DIM; d++) {
2950: PetscInt startF, endF, e;
2952: PetscCall(DMPlexGetSimplexOrBoxCells(plexF, P4EST_DIM - d, &startF, &endF));
2953: for (e = startF; e < endF; e++) PetscCall(PetscSectionSetDof(leafSection, e, dofPerDim[d]));
2954: }
2956: PetscCall(PetscSectionSetUp(rootSection));
2957: PetscCall(PetscSectionSetUp(leafSection));
2958: {
2959: PetscInt nroots, nleaves;
2960: PetscInt *mine, i, p;
2961: PetscInt *offsets, *offsetsRoot;
2962: PetscSFNode *remote;
2964: PetscCall(PetscMalloc1(pEndF - pStartF, &offsets));
2965: PetscCall(PetscMalloc1(pEndC - pStartC, &offsetsRoot));
2966: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionGetOffset(rootSection, p, &offsetsRoot[p - pStartC]));
2967: PetscCall(PetscSFBcastBegin(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2968: PetscCall(PetscSFBcastEnd(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2969: PetscCall(PetscSectionGetStorageSize(rootSection, &nroots));
2970: nleaves = 0;
2971: for (i = 0; i < numLeaves; i++) {
2972: PetscInt leaf = leaves ? leaves[i] : i;
2973: PetscInt dof;
2975: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2976: nleaves += dof;
2977: }
2978: PetscCall(PetscMalloc1(nleaves, &mine));
2979: PetscCall(PetscMalloc1(nleaves, &remote));
2980: nleaves = 0;
2981: for (i = 0; i < numLeaves; i++) {
2982: PetscInt leaf = leaves ? leaves[i] : i;
2983: PetscInt dof;
2984: PetscInt off, j;
2986: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2987: PetscCall(PetscSectionGetOffset(leafSection, leaf, &off));
2988: for (j = 0; j < dof; j++) {
2989: remote[nleaves].rank = iremote[i].rank;
2990: remote[nleaves].index = offsets[leaf] + j;
2991: mine[nleaves++] = off + j;
2992: }
2993: }
2994: PetscCall(PetscFree(offsetsRoot));
2995: PetscCall(PetscFree(offsets));
2996: PetscCall(PetscSFCreate(comm, sf));
2997: PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
2998: }
2999: PetscCall(PetscSectionDestroy(&leafSection));
3000: PetscCall(PetscSectionDestroy(&rootSection));
3001: PetscCall(PetscSFDestroy(&pointTransferSF));
3002: }
3003: PetscFunctionReturn(PETSC_SUCCESS);
3004: }
3006: static PetscErrorCode DMPforestGetTransferSF(DM dmA, DM dmB, const PetscInt dofPerDim[], PetscSF *sfAtoB, PetscSF *sfBtoA)
3007: {
3008: DM adaptA, adaptB;
3009: DMAdaptFlag purpose;
3011: PetscFunctionBegin;
3012: PetscCall(DMForestGetAdaptivityForest(dmA, &adaptA));
3013: PetscCall(DMForestGetAdaptivityForest(dmB, &adaptB));
3014: /* it is more efficient when the coarser mesh is the first argument: reorder if we know one is coarser than the other */
3015: if (adaptA && adaptA->data == dmB->data) { /* dmA was adapted from dmB */
3016: PetscCall(DMForestGetAdaptivityPurpose(dmA, &purpose));
3017: if (purpose == DM_ADAPT_REFINE) {
3018: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3019: PetscFunctionReturn(PETSC_SUCCESS);
3020: }
3021: } else if (adaptB && adaptB->data == dmA->data) { /* dmB was adapted from dmA */
3022: PetscCall(DMForestGetAdaptivityPurpose(dmB, &purpose));
3023: if (purpose == DM_ADAPT_COARSEN) {
3024: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3025: PetscFunctionReturn(PETSC_SUCCESS);
3026: }
3027: }
3028: if (sfAtoB) PetscCall(DMPforestGetTransferSF_Internal(dmA, dmB, dofPerDim, sfAtoB, PETSC_TRUE, NULL));
3029: if (sfBtoA) PetscCall(DMPforestGetTransferSF_Internal(dmB, dmA, dofPerDim, sfBtoA, (PetscBool)(sfAtoB == NULL), NULL));
3030: PetscFunctionReturn(PETSC_SUCCESS);
3031: }
3033: static PetscErrorCode DMPforestLabelsInitialize(DM dm, DM plex)
3034: {
3035: DM_Forest *forest = (DM_Forest *)dm->data;
3036: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
3037: PetscInt cLocalStart, cLocalEnd, cStart, cEnd, fStart, fEnd, eStart, eEnd, vStart, vEnd;
3038: PetscInt cStartBase, cEndBase, fStartBase, fEndBase, vStartBase, vEndBase, eStartBase, eEndBase;
3039: PetscInt pStart, pEnd, pStartBase, pEndBase, p;
3040: DM base;
3041: PetscInt *star = NULL, starSize;
3042: DMLabelLink next = dm->labels;
3043: PetscInt guess = 0;
3044: p4est_topidx_t num_trees = pforest->topo->conn->num_trees;
3046: PetscFunctionBegin;
3047: pforest->labelsFinalized = PETSC_TRUE;
3048: cLocalStart = pforest->cLocalStart;
3049: cLocalEnd = pforest->cLocalEnd;
3050: PetscCall(DMForestGetBaseDM(dm, &base));
3051: if (!base) {
3052: if (pforest->ghostName) { /* insert a label to make the boundaries, with stratum values denoting which face of the element touches the boundary */
3053: p4est_connectivity_t *conn = pforest->topo->conn;
3054: p4est_t *p4est = pforest->forest;
3055: p4est_tree_t *trees = (p4est_tree_t *)p4est->trees->array;
3056: p4est_topidx_t t, flt = p4est->first_local_tree;
3057: p4est_topidx_t llt = pforest->forest->last_local_tree;
3058: DMLabel ghostLabel;
3059: PetscInt c;
3061: PetscCall(DMCreateLabel(plex, pforest->ghostName));
3062: PetscCall(DMGetLabel(plex, pforest->ghostName, &ghostLabel));
3063: for (c = cLocalStart, t = flt; t <= llt; t++) {
3064: p4est_tree_t *tree = &trees[t];
3065: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3066: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3067: PetscInt q;
3069: for (q = 0; q < numQuads; q++, c++) {
3070: p4est_quadrant_t *quad = &quads[q];
3071: PetscInt f;
3073: for (f = 0; f < P4EST_FACES; f++) {
3074: p4est_quadrant_t neigh;
3075: int isOutside;
3077: PetscCallP4est(p4est_quadrant_face_neighbor, (quad, f, &neigh));
3078: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&neigh));
3079: if (isOutside) {
3080: p4est_topidx_t nt;
3081: PetscInt nf;
3083: nt = conn->tree_to_tree[t * P4EST_FACES + f];
3084: nf = (PetscInt)conn->tree_to_face[t * P4EST_FACES + f];
3085: nf = nf % P4EST_FACES;
3086: if (nt == t && nf == f) {
3087: PetscInt plexF = P4estFaceToPetscFace[f];
3088: const PetscInt *cone;
3090: PetscCall(DMPlexGetCone(plex, c, &cone));
3091: PetscCall(DMLabelSetValue(ghostLabel, cone[plexF], plexF + 1));
3092: }
3093: }
3094: }
3095: }
3096: }
3097: }
3098: PetscFunctionReturn(PETSC_SUCCESS);
3099: }
3100: PetscCall(DMPlexGetSimplexOrBoxCells(base, 0, &cStartBase, &cEndBase));
3101: PetscCall(DMPlexGetSimplexOrBoxCells(base, 1, &fStartBase, &fEndBase));
3102: PetscCall(DMPlexGetSimplexOrBoxCells(base, P4EST_DIM - 1, &eStartBase, &eEndBase));
3103: PetscCall(DMPlexGetDepthStratum(base, 0, &vStartBase, &vEndBase));
3105: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
3106: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 1, &fStart, &fEnd));
3107: PetscCall(DMPlexGetSimplexOrBoxCells(plex, P4EST_DIM - 1, &eStart, &eEnd));
3108: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3110: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3111: PetscCall(DMPlexGetChart(base, &pStartBase, &pEndBase));
3112: /* go through the mesh: use star to find a quadrant that borders a point. Use the closure to determine the
3113: * orientation of the quadrant relative to that point. Use that to relate the point to the numbering in the base
3114: * mesh, and extract a label value (since the base mesh is redundantly distributed, can be found locally). */
3115: while (next) {
3116: DMLabel baseLabel;
3117: DMLabel label = next->label;
3118: PetscBool isDepth, isCellType, isGhost, isVTK, isSpmap;
3119: const char *name;
3121: PetscCall(PetscObjectGetName((PetscObject)label, &name));
3122: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3123: if (isDepth) {
3124: next = next->next;
3125: continue;
3126: }
3127: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3128: if (isCellType) {
3129: next = next->next;
3130: continue;
3131: }
3132: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3133: if (isGhost) {
3134: next = next->next;
3135: continue;
3136: }
3137: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3138: if (isVTK) {
3139: next = next->next;
3140: continue;
3141: }
3142: PetscCall(PetscStrcmp(name, "_forest_base_subpoint_map", &isSpmap));
3143: if (!isSpmap) {
3144: PetscCall(DMGetLabel(base, name, &baseLabel));
3145: if (!baseLabel) {
3146: next = next->next;
3147: continue;
3148: }
3149: PetscCall(DMLabelCreateIndex(baseLabel, pStartBase, pEndBase));
3150: } else baseLabel = NULL;
3152: for (p = pStart; p < pEnd; p++) {
3153: PetscInt s, c = -1, l;
3154: PetscInt *closure = NULL, closureSize;
3155: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3156: p4est_tree_t *trees = (p4est_tree_t *)pforest->forest->trees->array;
3157: p4est_quadrant_t *q;
3158: PetscInt t, val;
3159: PetscBool zerosupportpoint = PETSC_FALSE;
3161: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3162: for (s = 0; s < starSize; s++) {
3163: PetscInt point = star[2 * s];
3165: if (cStart <= point && point < cEnd) {
3166: PetscCall(DMPlexGetTransitiveClosure(plex, point, PETSC_TRUE, &closureSize, &closure));
3167: for (l = 0; l < closureSize; l++) {
3168: PetscInt qParent = closure[2 * l], q, pp = p, pParent = p;
3169: do { /* check parents of q */
3170: q = qParent;
3171: if (q == p) {
3172: c = point;
3173: break;
3174: }
3175: PetscCall(DMPlexGetTreeParent(plex, q, &qParent, NULL));
3176: } while (qParent != q);
3177: if (c != -1) break;
3178: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3179: q = closure[2 * l];
3180: while (pParent != pp) { /* check parents of p */
3181: pp = pParent;
3182: if (pp == q) {
3183: c = point;
3184: break;
3185: }
3186: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3187: }
3188: if (c != -1) break;
3189: }
3190: PetscCall(DMPlexRestoreTransitiveClosure(plex, point, PETSC_TRUE, NULL, &closure));
3191: if (l < closureSize) break;
3192: } else {
3193: PetscInt supportSize;
3195: PetscCall(DMPlexGetSupportSize(plex, point, &supportSize));
3196: zerosupportpoint = (PetscBool)(zerosupportpoint || !supportSize);
3197: }
3198: }
3199: if (c < 0) {
3200: const char *prefix;
3201: PetscBool print = PETSC_FALSE;
3203: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
3204: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_forest_print_label_error", &print, NULL));
3205: if (print) {
3206: PetscInt i;
3208: 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));
3209: 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]));
3210: }
3211: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3212: if (zerosupportpoint) continue;
3213: else
3214: 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");
3215: }
3216: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3218: if (c < cLocalStart) {
3219: /* get from the beginning of the ghost layer */
3220: q = &ghosts[c];
3221: t = (PetscInt)q->p.which_tree;
3222: } else if (c < cLocalEnd) {
3223: PetscInt lo = 0, hi = num_trees;
3224: /* get from local quadrants: have to find the right tree */
3226: c -= cLocalStart;
3228: do {
3229: p4est_tree_t *tree;
3231: PetscCheck(guess >= lo && guess < num_trees && lo < hi, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed binary search");
3232: tree = &trees[guess];
3233: if (c < tree->quadrants_offset) {
3234: hi = guess;
3235: } else if (c < tree->quadrants_offset + (PetscInt)tree->quadrants.elem_count) {
3236: q = &((p4est_quadrant_t *)tree->quadrants.array)[c - (PetscInt)tree->quadrants_offset];
3237: t = guess;
3238: break;
3239: } else {
3240: lo = guess + 1;
3241: }
3242: guess = lo + (hi - lo) / 2;
3243: } while (1);
3244: } else {
3245: /* get from the end of the ghost layer */
3246: c -= (cLocalEnd - cLocalStart);
3248: q = &ghosts[c];
3249: t = (PetscInt)q->p.which_tree;
3250: }
3252: if (l == 0) { /* cell */
3253: if (baseLabel) {
3254: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3255: } else {
3256: val = t + cStartBase;
3257: }
3258: PetscCall(DMLabelSetValue(label, p, val));
3259: } else if (l >= 1 && l < 1 + P4EST_FACES) { /* facet */
3260: p4est_quadrant_t nq;
3261: int isInside;
3263: l = PetscFaceToP4estFace[l - 1];
3264: PetscCallP4est(p4est_quadrant_face_neighbor, (q, l, &nq));
3265: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3266: if (isInside) {
3267: /* this facet is in the interior of a tree, so it inherits the label of the tree */
3268: if (baseLabel) {
3269: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3270: } else {
3271: val = t + cStartBase;
3272: }
3273: PetscCall(DMLabelSetValue(label, p, val));
3274: } else {
3275: PetscInt f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + l];
3277: if (baseLabel) {
3278: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3279: } else {
3280: val = f + fStartBase;
3281: }
3282: PetscCall(DMLabelSetValue(label, p, val));
3283: }
3284: #if defined(P4_TO_P8)
3285: } else if (l >= 1 + P4EST_FACES && l < 1 + P4EST_FACES + P8EST_EDGES) { /* edge */
3286: p4est_quadrant_t nq;
3287: int isInside;
3289: l = PetscEdgeToP4estEdge[l - (1 + P4EST_FACES)];
3290: PetscCallP4est(p8est_quadrant_edge_neighbor, (q, l, &nq));
3291: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3292: if (isInside) {
3293: /* this edge is in the interior of a tree, so it inherits the label of the tree */
3294: if (baseLabel) {
3295: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3296: } else {
3297: val = t + cStartBase;
3298: }
3299: PetscCall(DMLabelSetValue(label, p, val));
3300: } else {
3301: int isOutsideFace;
3303: PetscCallP4estReturn(isOutsideFace, p4est_quadrant_is_outside_face, (&nq));
3304: if (isOutsideFace) {
3305: PetscInt f;
3307: if (nq.x < 0) {
3308: f = 0;
3309: } else if (nq.x >= P4EST_ROOT_LEN) {
3310: f = 1;
3311: } else if (nq.y < 0) {
3312: f = 2;
3313: } else if (nq.y >= P4EST_ROOT_LEN) {
3314: f = 3;
3315: } else if (nq.z < 0) {
3316: f = 4;
3317: } else {
3318: f = 5;
3319: }
3320: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3321: if (baseLabel) {
3322: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3323: } else {
3324: val = f + fStartBase;
3325: }
3326: PetscCall(DMLabelSetValue(label, p, val));
3327: } else { /* the quadrant edge corresponds to the tree edge */
3328: PetscInt e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + l];
3330: if (baseLabel) {
3331: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3332: } else {
3333: val = e + eStartBase;
3334: }
3335: PetscCall(DMLabelSetValue(label, p, val));
3336: }
3337: }
3338: #endif
3339: } else { /* vertex */
3340: p4est_quadrant_t nq;
3341: int isInside;
3343: #if defined(P4_TO_P8)
3344: l = PetscVertToP4estVert[l - (1 + P4EST_FACES + P8EST_EDGES)];
3345: #else
3346: l = PetscVertToP4estVert[l - (1 + P4EST_FACES)];
3347: #endif
3348: PetscCallP4est(p4est_quadrant_corner_neighbor, (q, l, &nq));
3349: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3350: if (isInside) {
3351: if (baseLabel) {
3352: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3353: } else {
3354: val = t + cStartBase;
3355: }
3356: PetscCall(DMLabelSetValue(label, p, val));
3357: } else {
3358: int isOutside;
3360: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&nq));
3361: if (isOutside) {
3362: PetscInt f = -1;
3364: if (nq.x < 0) {
3365: f = 0;
3366: } else if (nq.x >= P4EST_ROOT_LEN) {
3367: f = 1;
3368: } else if (nq.y < 0) {
3369: f = 2;
3370: } else if (nq.y >= P4EST_ROOT_LEN) {
3371: f = 3;
3372: #if defined(P4_TO_P8)
3373: } else if (nq.z < 0) {
3374: f = 4;
3375: } else {
3376: f = 5;
3377: #endif
3378: }
3379: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3380: if (baseLabel) {
3381: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3382: } else {
3383: val = f + fStartBase;
3384: }
3385: PetscCall(DMLabelSetValue(label, p, val));
3386: continue;
3387: }
3388: #if defined(P4_TO_P8)
3389: PetscCallP4estReturn(isOutside, p8est_quadrant_is_outside_edge, (&nq));
3390: if (isOutside) {
3391: /* outside edge */
3392: PetscInt e = -1;
3394: if (nq.x >= 0 && nq.x < P4EST_ROOT_LEN) {
3395: if (nq.z < 0) {
3396: if (nq.y < 0) {
3397: e = 0;
3398: } else {
3399: e = 1;
3400: }
3401: } else {
3402: if (nq.y < 0) {
3403: e = 2;
3404: } else {
3405: e = 3;
3406: }
3407: }
3408: } else if (nq.y >= 0 && nq.y < P4EST_ROOT_LEN) {
3409: if (nq.z < 0) {
3410: if (nq.x < 0) {
3411: e = 4;
3412: } else {
3413: e = 5;
3414: }
3415: } else {
3416: if (nq.x < 0) {
3417: e = 6;
3418: } else {
3419: e = 7;
3420: }
3421: }
3422: } else {
3423: if (nq.y < 0) {
3424: if (nq.x < 0) {
3425: e = 8;
3426: } else {
3427: e = 9;
3428: }
3429: } else {
3430: if (nq.x < 0) {
3431: e = 10;
3432: } else {
3433: e = 11;
3434: }
3435: }
3436: }
3438: e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + e];
3439: if (baseLabel) {
3440: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3441: } else {
3442: val = e + eStartBase;
3443: }
3444: PetscCall(DMLabelSetValue(label, p, val));
3445: continue;
3446: }
3447: #endif
3448: {
3449: /* outside vertex: same corner as quadrant corner */
3450: PetscInt v = pforest->topo->conn->tree_to_corner[P4EST_CHILDREN * t + l];
3452: if (baseLabel) {
3453: PetscCall(DMLabelGetValue(baseLabel, v + vStartBase, &val));
3454: } else {
3455: val = v + vStartBase;
3456: }
3457: PetscCall(DMLabelSetValue(label, p, val));
3458: }
3459: }
3460: }
3461: }
3462: next = next->next;
3463: }
3464: PetscFunctionReturn(PETSC_SUCCESS);
3465: }
3467: static PetscErrorCode DMPforestLabelsFinalize(DM dm, DM plex)
3468: {
3469: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
3470: DM adapt;
3472: PetscFunctionBegin;
3473: if (pforest->labelsFinalized) PetscFunctionReturn(PETSC_SUCCESS);
3474: pforest->labelsFinalized = PETSC_TRUE;
3475: PetscCall(DMForestGetAdaptivityForest(dm, &adapt));
3476: if (!adapt) {
3477: /* Initialize labels from the base dm */
3478: PetscCall(DMPforestLabelsInitialize(dm, plex));
3479: } else {
3480: PetscInt dofPerDim[4] = {1, 1, 1, 1};
3481: PetscSF transferForward, transferBackward, pointSF;
3482: PetscInt pStart, pEnd, pStartA, pEndA;
3483: PetscInt *values, *adaptValues;
3484: DMLabelLink next = adapt->labels;
3485: DMLabel adaptLabel;
3486: DM adaptPlex;
3488: PetscCall(DMForestGetAdaptivityLabel(dm, &adaptLabel));
3489: PetscCall(DMPforestGetPlex(adapt, &adaptPlex));
3490: PetscCall(DMPforestGetTransferSF(adapt, dm, dofPerDim, &transferForward, &transferBackward));
3491: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3492: PetscCall(DMPlexGetChart(adaptPlex, &pStartA, &pEndA));
3493: PetscCall(PetscMalloc2(pEnd - pStart, &values, pEndA - pStartA, &adaptValues));
3494: PetscCall(DMGetPointSF(plex, &pointSF));
3495: if (PetscDefined(USE_DEBUG)) {
3496: PetscInt p;
3497: for (p = pStartA; p < pEndA; p++) adaptValues[p - pStartA] = -1;
3498: for (p = pStart; p < pEnd; p++) values[p - pStart] = -2;
3499: if (transferForward) {
3500: PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3501: PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3502: }
3503: if (transferBackward) {
3504: PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3505: PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3506: }
3507: for (p = pStart; p < pEnd; p++) {
3508: PetscInt q = p, parent;
3510: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3511: while (parent != q) {
3512: if (values[parent] == -2) values[parent] = values[q];
3513: q = parent;
3514: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3515: }
3516: }
3517: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3518: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3519: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3520: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3521: for (p = pStart; p < pEnd; p++) PetscCheck(values[p - pStart] != -2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "uncovered point %" PetscInt_FMT, p);
3522: }
3523: while (next) {
3524: DMLabel nextLabel = next->label;
3525: const char *name;
3526: PetscBool isDepth, isCellType, isGhost, isVTK;
3527: DMLabel label;
3528: PetscInt p;
3530: PetscCall(PetscObjectGetName((PetscObject)nextLabel, &name));
3531: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3532: if (isDepth) {
3533: next = next->next;
3534: continue;
3535: }
3536: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3537: if (isCellType) {
3538: next = next->next;
3539: continue;
3540: }
3541: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3542: if (isGhost) {
3543: next = next->next;
3544: continue;
3545: }
3546: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3547: if (isVTK) {
3548: next = next->next;
3549: continue;
3550: }
3551: if (nextLabel == adaptLabel) {
3552: next = next->next;
3553: continue;
3554: }
3555: /* label was created earlier */
3556: PetscCall(DMGetLabel(dm, name, &label));
3557: for (p = pStartA; p < pEndA; p++) PetscCall(DMLabelGetValue(nextLabel, p, &adaptValues[p]));
3558: for (p = pStart; p < pEnd; p++) values[p] = PETSC_INT_MIN;
3560: if (transferForward) PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3561: if (transferBackward) PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3562: if (transferForward) PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3563: if (transferBackward) PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3564: for (p = pStart; p < pEnd; p++) {
3565: PetscInt q = p, parent;
3567: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3568: while (parent != q) {
3569: if (values[parent] == PETSC_INT_MIN) values[parent] = values[q];
3570: q = parent;
3571: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3572: }
3573: }
3574: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3575: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3576: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3577: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3579: for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(label, p, values[p]));
3580: next = next->next;
3581: }
3582: PetscCall(PetscFree2(values, adaptValues));
3583: PetscCall(PetscSFDestroy(&transferForward));
3584: PetscCall(PetscSFDestroy(&transferBackward));
3585: pforest->labelsFinalized = PETSC_TRUE;
3586: }
3587: PetscFunctionReturn(PETSC_SUCCESS);
3588: }
3590: 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)
3591: {
3592: PetscInt closureSize, c, coordStart, coordEnd, coordDim;
3593: PetscInt *closure = NULL;
3594: PetscSection coordSec;
3596: PetscFunctionBegin;
3597: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3598: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3599: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3600: PetscCall(DMPlexGetTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3601: for (c = 0; c < closureSize; c++) {
3602: PetscInt point = closure[2 * c];
3604: if (point >= coordStart && point < coordEnd) {
3605: PetscInt dof, off;
3606: PetscInt nCoords, i;
3607: PetscCall(PetscSectionGetDof(coordSec, point, &dof));
3608: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3609: nCoords = dof / coordDim;
3610: PetscCall(PetscSectionGetOffset(coordSec, point, &off));
3611: for (i = 0; i < nCoords; i++) {
3612: PetscScalar *coord = &coords[off + i * coordDim];
3613: double coordP4est[3] = {0.};
3614: double coordP4estMapped[3] = {0.};
3615: PetscInt j;
3616: PetscReal treeCoords[P4EST_CHILDREN][3] = {{0.}};
3617: PetscReal eta[3] = {0.};
3618: PetscInt numRounds = 10;
3619: PetscReal coordGuess[3] = {0.};
3621: eta[0] = (PetscReal)q->x / (PetscReal)P4EST_ROOT_LEN;
3622: eta[1] = (PetscReal)q->y / (PetscReal)P4EST_ROOT_LEN;
3623: #if defined(P4_TO_P8)
3624: eta[2] = (PetscReal)q->z / (PetscReal)P4EST_ROOT_LEN;
3625: #endif
3627: for (j = 0; j < P4EST_CHILDREN; j++) {
3628: PetscInt k;
3630: for (k = 0; k < 3; k++) treeCoords[j][k] = conn->vertices[3 * conn->tree_to_vertex[P4EST_CHILDREN * t + j] + k];
3631: }
3633: for (j = 0; j < P4EST_CHILDREN; j++) {
3634: PetscInt k;
3635: PetscReal prod = 1.;
3637: for (k = 0; k < P4EST_DIM; k++) prod *= (j & (1 << k)) ? eta[k] : (1. - eta[k]);
3638: for (k = 0; k < 3; k++) coordGuess[k] += prod * treeCoords[j][k];
3639: }
3641: for (j = 0; j < numRounds; j++) {
3642: PetscInt dir;
3644: for (dir = 0; dir < P4EST_DIM; dir++) {
3645: PetscInt k;
3646: PetscReal diff[3];
3647: PetscReal dXdeta[3] = {0.};
3648: PetscReal rhs, scale, update;
3650: for (k = 0; k < 3; k++) diff[k] = coordP4est[k] - coordGuess[k];
3651: for (k = 0; k < P4EST_CHILDREN; k++) {
3652: PetscInt l;
3653: PetscReal prod = 1.;
3655: for (l = 0; l < P4EST_DIM; l++) {
3656: if (l == dir) {
3657: prod *= (k & (1 << l)) ? 1. : -1.;
3658: } else {
3659: prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3660: }
3661: }
3662: for (l = 0; l < 3; l++) dXdeta[l] += prod * treeCoords[k][l];
3663: }
3664: rhs = 0.;
3665: scale = 0;
3666: for (k = 0; k < 3; k++) {
3667: rhs += diff[k] * dXdeta[k];
3668: scale += dXdeta[k] * dXdeta[k];
3669: }
3670: update = rhs / scale;
3671: eta[dir] += update;
3672: eta[dir] = PetscMin(eta[dir], 1.);
3673: eta[dir] = PetscMax(eta[dir], 0.);
3675: coordGuess[0] = coordGuess[1] = coordGuess[2] = 0.;
3676: for (k = 0; k < P4EST_CHILDREN; k++) {
3677: PetscInt l;
3678: PetscReal prod = 1.;
3680: for (l = 0; l < P4EST_DIM; l++) prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3681: for (l = 0; l < 3; l++) coordGuess[l] += prod * treeCoords[k][l];
3682: }
3683: }
3684: }
3685: for (j = 0; j < 3; j++) coordP4est[j] = (double)eta[j];
3687: PetscCheck(geom, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded");
3688: (geom->X)(geom, t, coordP4est, coordP4estMapped);
3689: for (j = 0; j < coordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3690: }
3691: }
3692: }
3693: PetscCall(DMPlexRestoreTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3694: PetscFunctionReturn(PETSC_SUCCESS);
3695: }
3697: static PetscErrorCode DMPforestMapCoordinates(DM dm, DM plex)
3698: {
3699: DM_Forest *forest;
3700: DM_Forest_pforest *pforest;
3701: p4est_geometry_t *geom;
3702: PetscInt cLocalStart, cLocalEnd;
3703: Vec coordLocalVec;
3704: PetscScalar *coords;
3705: p4est_topidx_t flt, llt, t;
3706: p4est_tree_t *trees;
3707: PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
3708: void *mapCtx;
3710: PetscFunctionBegin;
3711: forest = (DM_Forest *)dm->data;
3712: pforest = (DM_Forest_pforest *)forest->data;
3713: geom = pforest->topo->geom;
3714: PetscCall(DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx));
3715: if (!geom && !map) PetscFunctionReturn(PETSC_SUCCESS);
3716: PetscCall(DMGetCoordinatesLocal(plex, &coordLocalVec));
3717: PetscCall(VecGetArray(coordLocalVec, &coords));
3718: cLocalStart = pforest->cLocalStart;
3719: cLocalEnd = pforest->cLocalEnd;
3720: flt = pforest->forest->first_local_tree;
3721: llt = pforest->forest->last_local_tree;
3722: trees = (p4est_tree_t *)pforest->forest->trees->array;
3723: if (map) { /* apply the map directly to the existing coordinates */
3724: PetscSection coordSec;
3725: PetscInt coordStart, coordEnd, p, coordDim, p4estCoordDim, cStart, cEnd, cEndInterior;
3726: DM base;
3728: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3729: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3730: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3731: PetscCall(DMForestGetBaseDM(dm, &base));
3732: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3733: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3734: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3735: p4estCoordDim = PetscMin(coordDim, 3);
3736: for (p = coordStart; p < coordEnd; p++) {
3737: PetscInt *star = NULL, starSize;
3738: PetscInt dof, off, cell = -1, coarsePoint = -1;
3739: PetscInt nCoords, i;
3740: PetscCall(PetscSectionGetDof(coordSec, p, &dof));
3741: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3742: nCoords = dof / coordDim;
3743: PetscCall(PetscSectionGetOffset(coordSec, p, &off));
3744: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3745: for (i = 0; i < starSize; i++) {
3746: PetscInt point = star[2 * i];
3748: if (cStart <= point && point < cEnd) {
3749: cell = point;
3750: break;
3751: }
3752: }
3753: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3754: if (cell >= 0) {
3755: if (cell < cLocalStart) {
3756: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3758: coarsePoint = ghosts[cell].p.which_tree;
3759: } else if (cell < cLocalEnd) {
3760: cell -= cLocalStart;
3761: for (t = flt; t <= llt; t++) {
3762: p4est_tree_t *tree = &trees[t];
3764: if (cell >= tree->quadrants_offset && (size_t)cell < tree->quadrants_offset + tree->quadrants.elem_count) {
3765: coarsePoint = t;
3766: break;
3767: }
3768: }
3769: } else {
3770: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3772: coarsePoint = ghosts[cell - cLocalEnd].p.which_tree;
3773: }
3774: }
3775: for (i = 0; i < nCoords; i++) {
3776: PetscScalar *coord = &coords[off + i * coordDim];
3777: PetscReal coordP4est[3] = {0.};
3778: PetscReal coordP4estMapped[3] = {0.};
3779: PetscInt j;
3781: for (j = 0; j < p4estCoordDim; j++) coordP4est[j] = PetscRealPart(coord[j]);
3782: PetscCall((map)(base, coarsePoint, p4estCoordDim, coordP4est, coordP4estMapped, mapCtx));
3783: for (j = 0; j < p4estCoordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3784: }
3785: }
3786: } else { /* we have to transform coordinates back to the unit cube (where geom is defined), and then apply geom */
3787: PetscInt cStart, cEnd, cEndInterior;
3789: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3790: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3791: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3792: if (cLocalStart > 0) {
3793: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3794: PetscInt count;
3796: for (count = 0; count < cLocalStart; count++) {
3797: p4est_quadrant_t *quad = &ghosts[count];
3798: p4est_topidx_t t = quad->p.which_tree;
3800: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, quad, t, pforest->topo->conn, coords));
3801: }
3802: }
3803: for (t = flt; t <= llt; t++) {
3804: p4est_tree_t *tree = &trees[t];
3805: PetscInt offset = cLocalStart + tree->quadrants_offset, i;
3806: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3807: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3809: for (i = 0; i < numQuads; i++) {
3810: PetscInt count = i + offset;
3812: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, &quads[i], t, pforest->topo->conn, coords));
3813: }
3814: }
3815: if (cLocalEnd - cLocalStart < cEnd - cStart) {
3816: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3817: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
3818: PetscInt count;
3820: for (count = 0; count < numGhosts - cLocalStart; count++) {
3821: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
3822: p4est_topidx_t t = quad->p.which_tree;
3824: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count + cLocalEnd, quad, t, pforest->topo->conn, coords));
3825: }
3826: }
3827: }
3828: PetscCall(VecRestoreArray(coordLocalVec, &coords));
3829: PetscFunctionReturn(PETSC_SUCCESS);
3830: }
3832: static PetscErrorCode PforestQuadrantIsInterior(p4est_quadrant_t *quad, PetscBool *is_interior)
3833: {
3834: PetscFunctionBegin;
3835: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3836: if ((quad->x > 0) && (quad->x + h < P4EST_ROOT_LEN)
3837: #ifdef P4_TO_P8
3838: && (quad->z > 0) && (quad->z + h < P4EST_ROOT_LEN)
3839: #endif
3840: && (quad->y > 0) && (quad->y + h < P4EST_ROOT_LEN)) {
3841: *is_interior = PETSC_TRUE;
3842: } else {
3843: *is_interior = PETSC_FALSE;
3844: }
3845: PetscFunctionReturn(PETSC_SUCCESS);
3846: }
3848: /* We always use DG coordinates with p4est: if they do not match the vertex
3849: coordinates, add space for them in the section */
3850: 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)
3851: {
3852: PetscBool is_interior;
3854: PetscFunctionBegin;
3855: PetscCall(PforestQuadrantIsInterior(quad, &is_interior));
3856: if (is_interior) { // quads in the interior of a coarse cell can't touch periodic interfaces
3857: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3858: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3859: } else {
3860: PetscInt cSize;
3861: PetscScalar *values = NULL;
3862: PetscBool same_coords = PETSC_TRUE;
3864: PetscCall(DMPlexVecGetClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3865: PetscAssert(cSize == cDim * P4EST_CHILDREN, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected closure size");
3866: for (int c = 0; c < P4EST_CHILDREN; c++) {
3867: p4est_qcoord_t quad_coords[3];
3868: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3869: double corner_coords[3];
3870: double vert_coords[3];
3871: PetscInt corner = PetscVertToP4estVert[c];
3873: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) vert_coords[d] = PetscRealPart(values[c * cDim + d]);
3875: quad_coords[0] = quad->x;
3876: quad_coords[1] = quad->y;
3877: #ifdef P4_TO_P8
3878: quad_coords[2] = quad->z;
3879: #endif
3880: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3881: #ifndef P4_TO_P8
3882: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3883: #else
3884: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3885: #endif
3886: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) {
3887: if (fabs(vert_coords[d] - corner_coords[d]) > PETSC_SMALL) {
3888: same_coords = PETSC_FALSE;
3889: break;
3890: }
3891: }
3892: }
3893: if (same_coords) {
3894: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3895: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3896: } else {
3897: PetscCall(PetscSectionSetDof(newSection, cell, cSize));
3898: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, cSize));
3899: }
3900: PetscCall(DMPlexVecRestoreClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3901: }
3902: PetscFunctionReturn(PETSC_SUCCESS);
3903: }
3905: static PetscErrorCode PforestLocalizeCell(DM plex, PetscInt cDim, DM_Forest_pforest *pforest, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad, PetscScalar coords[])
3906: {
3907: PetscInt cdof, off;
3909: PetscFunctionBegin;
3910: PetscCall(PetscSectionGetDof(newSection, cell, &cdof));
3911: if (!cdof) PetscFunctionReturn(PETSC_SUCCESS);
3913: PetscCall(PetscSectionGetOffset(newSection, cell, &off));
3914: for (PetscInt c = 0, pos = off; c < P4EST_CHILDREN; c++) {
3915: p4est_qcoord_t quad_coords[3];
3916: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3917: double corner_coords[3];
3918: PetscInt corner = PetscVertToP4estVert[c];
3920: quad_coords[0] = quad->x;
3921: quad_coords[1] = quad->y;
3922: #ifdef P4_TO_P8
3923: quad_coords[2] = quad->z;
3924: #endif
3925: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3926: #ifndef P4_TO_P8
3927: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3928: #else
3929: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3930: #endif
3931: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) coords[pos++] = corner_coords[d];
3932: for (PetscInt d = PetscMin(cDim, 3); d < cDim; d++) coords[pos++] = 0.;
3933: }
3934: PetscFunctionReturn(PETSC_SUCCESS);
3935: }
3937: static PetscErrorCode DMPforestLocalizeCoordinates(DM dm, DM plex)
3938: {
3939: DM_Forest *forest;
3940: DM_Forest_pforest *pforest;
3941: DM base, cdm, cdmCell;
3942: Vec cVec, cVecOld;
3943: PetscSection oldSection, newSection;
3944: PetscScalar *coords2;
3945: const PetscReal *L;
3946: PetscInt cLocalStart, cLocalEnd, coarsePoint;
3947: PetscInt cDim, newStart, newEnd;
3948: PetscInt v, vStart, vEnd, cp, cStart, cEnd, cEndInterior;
3949: p4est_topidx_t flt, llt, t;
3950: p4est_tree_t *trees;
3951: PetscBool baseLocalized = PETSC_FALSE;
3953: PetscFunctionBegin;
3954: PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
3955: /* we localize on all cells if we don't have a base DM or the base DM coordinates have not been localized */
3956: PetscCall(DMGetCoordinateDim(dm, &cDim));
3957: PetscCall(DMForestGetBaseDM(dm, &base));
3958: if (base) PetscCall(DMGetCoordinatesLocalized(base, &baseLocalized));
3959: if (!baseLocalized) base = NULL;
3960: if (!baseLocalized && !L) PetscFunctionReturn(PETSC_SUCCESS);
3961: PetscCall(DMPlexGetChart(plex, &newStart, &newEnd));
3963: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newSection));
3964: PetscCall(PetscSectionSetNumFields(newSection, 1));
3965: PetscCall(PetscSectionSetFieldComponents(newSection, 0, cDim));
3966: PetscCall(PetscSectionSetChart(newSection, newStart, newEnd));
3968: PetscCall(DMGetCoordinateSection(plex, &oldSection));
3969: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3970: PetscCall(DMGetCoordinatesLocal(plex, &cVecOld));
3972: forest = (DM_Forest *)dm->data;
3973: pforest = (DM_Forest_pforest *)forest->data;
3974: cLocalStart = pforest->cLocalStart;
3975: cLocalEnd = pforest->cLocalEnd;
3976: flt = pforest->forest->first_local_tree;
3977: llt = pforest->forest->last_local_tree;
3978: trees = (p4est_tree_t *)pforest->forest->trees->array;
3980: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3981: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3982: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3983: cp = 0;
3984: if (cLocalStart > 0) {
3985: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3986: PetscInt cell;
3988: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
3989: p4est_quadrant_t *quad = &ghosts[cell];
3991: coarsePoint = quad->p.which_tree;
3992: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
3993: }
3994: }
3995: for (t = flt; t <= llt; t++) {
3996: p4est_tree_t *tree = &trees[t];
3997: PetscInt offset = cLocalStart + tree->quadrants_offset;
3998: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3999: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4000: PetscInt i;
4002: if (!numQuads) continue;
4003: coarsePoint = t;
4004: for (i = 0; i < numQuads; i++, cp++) {
4005: PetscInt cell = i + offset;
4006: p4est_quadrant_t *quad = &quads[i];
4008: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4009: }
4010: }
4011: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4012: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4013: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4014: PetscInt count;
4016: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4017: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4018: coarsePoint = quad->p.which_tree;
4019: PetscInt cell = count + cLocalEnd;
4021: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4022: }
4023: }
4024: PetscAssert(cp == cEnd - cStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of fine cells %" PetscInt_FMT " != %" PetscInt_FMT, cp, cEnd - cStart);
4026: PetscCall(PetscSectionSetUp(newSection));
4027: PetscCall(DMGetCoordinateDM(plex, &cdm));
4028: PetscCall(DMClone(cdm, &cdmCell));
4029: PetscCall(DMSetCellCoordinateDM(plex, cdmCell));
4030: PetscCall(DMDestroy(&cdmCell));
4031: PetscCall(DMSetCellCoordinateSection(plex, cDim, newSection));
4032: PetscCall(PetscSectionGetStorageSize(newSection, &v));
4033: PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
4034: PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
4035: PetscCall(VecSetBlockSize(cVec, cDim));
4036: PetscCall(VecSetSizes(cVec, v, PETSC_DETERMINE));
4037: PetscCall(VecSetType(cVec, VECSTANDARD));
4038: PetscCall(VecSet(cVec, PETSC_MIN_REAL));
4040: /* Localize coordinates on cells if needed */
4041: PetscCall(VecGetArray(cVec, &coords2));
4042: cp = 0;
4043: if (cLocalStart > 0) {
4044: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4045: PetscInt cell;
4047: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
4048: p4est_quadrant_t *quad = &ghosts[cell];
4050: coarsePoint = quad->p.which_tree;
4051: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4052: }
4053: }
4054: for (t = flt; t <= llt; t++) {
4055: p4est_tree_t *tree = &trees[t];
4056: PetscInt offset = cLocalStart + tree->quadrants_offset;
4057: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
4058: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4059: PetscInt i;
4061: if (!numQuads) continue;
4062: coarsePoint = t;
4063: for (i = 0; i < numQuads; i++, cp++) {
4064: PetscInt cell = i + offset;
4065: p4est_quadrant_t *quad = &quads[i];
4067: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4068: }
4069: }
4070: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4071: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4072: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4073: PetscInt count;
4075: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4076: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4077: coarsePoint = quad->p.which_tree;
4078: PetscInt cell = count + cLocalEnd;
4080: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4081: }
4082: }
4083: PetscCall(VecRestoreArray(cVec, &coords2));
4084: PetscCall(DMSetCellCoordinatesLocal(plex, cVec));
4085: PetscCall(VecDestroy(&cVec));
4086: PetscCall(PetscSectionDestroy(&newSection));
4087: PetscFunctionReturn(PETSC_SUCCESS);
4088: }
4090: #define DMForestClearAdaptivityForest_pforest _append_pforest(DMForestClearAdaptivityForest)
4091: static PetscErrorCode DMForestClearAdaptivityForest_pforest(DM dm)
4092: {
4093: DM_Forest *forest;
4094: DM_Forest_pforest *pforest;
4096: PetscFunctionBegin;
4097: forest = (DM_Forest *)dm->data;
4098: pforest = (DM_Forest_pforest *)forest->data;
4099: PetscCall(PetscSFDestroy(&pforest->pointAdaptToSelfSF));
4100: PetscCall(PetscSFDestroy(&pforest->pointSelfToAdaptSF));
4101: PetscCall(PetscFree(pforest->pointAdaptToSelfCids));
4102: PetscCall(PetscFree(pforest->pointSelfToAdaptCids));
4103: PetscFunctionReturn(PETSC_SUCCESS);
4104: }
4106: static PetscErrorCode DMConvert_pforest_plex(DM dm, DMType newtype, DM *plex)
4107: {
4108: DM_Forest *forest;
4109: DM_Forest_pforest *pforest;
4110: DM refTree, newPlex, base;
4111: PetscInt adjDim, adjCodim, coordDim;
4112: MPI_Comm comm;
4113: PetscBool isPforest;
4114: PetscInt dim;
4115: PetscInt overlap;
4116: p4est_connect_type_t ctype;
4117: p4est_locidx_t first_local_quad = -1;
4118: sc_array_t *points_per_dim, *cone_sizes, *cones, *cone_orientations, *coords, *children, *parents, *childids, *leaves, *remotes;
4119: PetscSection parentSection;
4120: PetscSF pointSF;
4121: size_t zz, count;
4122: PetscInt pStart, pEnd;
4123: DMLabel ghostLabelBase = NULL;
4125: PetscFunctionBegin;
4127: comm = PetscObjectComm((PetscObject)dm);
4128: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPFOREST, &isPforest));
4129: PetscCheck(isPforest, comm, PETSC_ERR_ARG_WRONG, "Expected DM type %s, got %s", DMPFOREST, ((PetscObject)dm)->type_name);
4130: PetscCall(DMSetUp(dm));
4131: PetscCall(DMGetDimension(dm, &dim));
4132: PetscCheck(dim == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Expected DM dimension %d, got %" PetscInt_FMT, P4EST_DIM, dim);
4133: forest = (DM_Forest *)dm->data;
4134: pforest = (DM_Forest_pforest *)forest->data;
4135: PetscCall(DMForestGetBaseDM(dm, &base));
4136: if (base) PetscCall(DMGetLabel(base, "ghost", &ghostLabelBase));
4137: if (!pforest->plex) {
4138: PetscMPIInt size;
4139: const char *name;
4141: PetscCallMPI(MPI_Comm_size(comm, &size));
4142: PetscCall(DMCreate(comm, &newPlex));
4143: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
4144: PetscCall(PetscObjectSetName((PetscObject)newPlex, name));
4145: PetscCall(DMSetType(newPlex, DMPLEX));
4146: PetscCall(DMSetMatType(newPlex, dm->mattype));
4147: /* share labels */
4148: PetscCall(DMCopyLabels(dm, newPlex, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4149: PetscCall(DMForestGetAdjacencyDimension(dm, &adjDim));
4150: PetscCall(DMForestGetAdjacencyCodimension(dm, &adjCodim));
4151: PetscCall(DMGetCoordinateDim(dm, &coordDim));
4152: if (adjDim == 0) {
4153: ctype = P4EST_CONNECT_FULL;
4154: } else if (adjCodim == 1) {
4155: ctype = P4EST_CONNECT_FACE;
4156: #if defined(P4_TO_P8)
4157: } else if (adjDim == 1) {
4158: ctype = P8EST_CONNECT_EDGE;
4159: #endif
4160: } else {
4161: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid adjacency dimension %" PetscInt_FMT, adjDim);
4162: }
4163: PetscCheck(ctype == P4EST_CONNECT_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Adjacency dimension %" PetscInt_FMT " / codimension %" PetscInt_FMT " not supported yet", adjDim, adjCodim);
4164: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
4165: PetscCall(DMPlexSetOverlap_Plex(newPlex, NULL, overlap));
4167: points_per_dim = sc_array_new(sizeof(p4est_locidx_t));
4168: cone_sizes = sc_array_new(sizeof(p4est_locidx_t));
4169: cones = sc_array_new(sizeof(p4est_locidx_t));
4170: cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
4171: coords = sc_array_new(3 * sizeof(double));
4172: children = sc_array_new(sizeof(p4est_locidx_t));
4173: parents = sc_array_new(sizeof(p4est_locidx_t));
4174: childids = sc_array_new(sizeof(p4est_locidx_t));
4175: leaves = sc_array_new(sizeof(p4est_locidx_t));
4176: remotes = sc_array_new(2 * sizeof(p4est_locidx_t));
4178: 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));
4180: pforest->cLocalStart = (PetscInt)first_local_quad;
4181: pforest->cLocalEnd = pforest->cLocalStart + (PetscInt)pforest->forest->local_num_quadrants;
4182: PetscCall(locidx_to_PetscInt(points_per_dim));
4183: PetscCall(locidx_to_PetscInt(cone_sizes));
4184: PetscCall(locidx_to_PetscInt(cones));
4185: PetscCall(locidx_to_PetscInt(cone_orientations));
4186: PetscCall(coords_double_to_PetscScalar(coords, coordDim));
4187: PetscCall(locidx_to_PetscInt(children));
4188: PetscCall(locidx_to_PetscInt(parents));
4189: PetscCall(locidx_to_PetscInt(childids));
4190: PetscCall(locidx_to_PetscInt(leaves));
4191: PetscCall(locidx_pair_to_PetscSFNode(remotes));
4193: PetscCall(DMSetDimension(newPlex, P4EST_DIM));
4194: PetscCall(DMSetCoordinateDim(newPlex, coordDim));
4195: PetscCall(DMPlexSetMaxProjectionHeight(newPlex, P4EST_DIM - 1));
4196: PetscCall(DMPlexCreateFromDAG(newPlex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array));
4197: PetscCall(DMPlexConvertOldOrientations_Internal(newPlex));
4198: PetscCall(DMCreateReferenceTree_pforest(comm, &refTree));
4199: PetscCall(DMPlexSetReferenceTree(newPlex, refTree));
4200: PetscCall(PetscSectionCreate(comm, &parentSection));
4201: PetscCall(DMPlexGetChart(newPlex, &pStart, &pEnd));
4202: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
4203: count = children->elem_count;
4204: for (zz = 0; zz < count; zz++) {
4205: PetscInt child = *((PetscInt *)sc_array_index(children, zz));
4207: PetscCall(PetscSectionSetDof(parentSection, child, 1));
4208: }
4209: PetscCall(PetscSectionSetUp(parentSection));
4210: PetscCall(DMPlexSetTree(newPlex, parentSection, (PetscInt *)parents->array, (PetscInt *)childids->array));
4211: PetscCall(PetscSectionDestroy(&parentSection));
4212: PetscCall(PetscSFCreate(comm, &pointSF));
4213: /*
4214: These arrays defining the sf are from the p4est library, but the code there shows the leaves being populated in increasing order.
4215: https://gitlab.com/petsc/petsc/merge_requests/2248#note_240186391
4216: */
4217: PetscCall(PetscSFSetGraph(pointSF, pEnd - pStart, (PetscInt)leaves->elem_count, (PetscInt *)leaves->array, PETSC_COPY_VALUES, (PetscSFNode *)remotes->array, PETSC_COPY_VALUES));
4218: PetscCall(DMSetPointSF(newPlex, pointSF));
4219: PetscCall(DMSetPointSF(dm, pointSF));
4220: {
4221: DM coordDM;
4223: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4224: PetscCall(DMSetPointSF(coordDM, pointSF));
4225: }
4226: PetscCall(PetscSFDestroy(&pointSF));
4227: sc_array_destroy(points_per_dim);
4228: sc_array_destroy(cone_sizes);
4229: sc_array_destroy(cones);
4230: sc_array_destroy(cone_orientations);
4231: sc_array_destroy(coords);
4232: sc_array_destroy(children);
4233: sc_array_destroy(parents);
4234: sc_array_destroy(childids);
4235: sc_array_destroy(leaves);
4236: sc_array_destroy(remotes);
4238: {
4239: const PetscReal *maxCell, *Lstart, *L;
4241: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
4242: PetscCall(DMSetPeriodicity(newPlex, maxCell, Lstart, L));
4243: PetscCall(DMPforestLocalizeCoordinates(dm, newPlex));
4244: }
4246: if (overlap > 0) { /* the p4est routine can't set all of the coordinates in its routine if there is overlap */
4247: Vec coordsGlobal, coordsLocal;
4248: const PetscScalar *globalArray;
4249: PetscScalar *localArray;
4250: PetscSF coordSF;
4251: DM coordDM;
4253: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4254: PetscCall(DMGetSectionSF(coordDM, &coordSF));
4255: PetscCall(DMGetCoordinates(newPlex, &coordsGlobal));
4256: PetscCall(DMGetCoordinatesLocal(newPlex, &coordsLocal));
4257: PetscCall(VecGetArrayRead(coordsGlobal, &globalArray));
4258: PetscCall(VecGetArray(coordsLocal, &localArray));
4259: PetscCall(PetscSFBcastBegin(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4260: PetscCall(PetscSFBcastEnd(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4261: PetscCall(VecRestoreArray(coordsLocal, &localArray));
4262: PetscCall(VecRestoreArrayRead(coordsGlobal, &globalArray));
4263: PetscCall(DMSetCoordinatesLocal(newPlex, coordsLocal));
4264: }
4265: PetscCall(DMPforestMapCoordinates(dm, newPlex));
4267: pforest->plex = newPlex;
4269: /* copy labels */
4270: PetscCall(DMPforestLabelsFinalize(dm, newPlex));
4272: if (ghostLabelBase || pforest->ghostName) { /* we have to do this after copying labels because the labels drive the construction of ghost cells */
4273: PetscInt numAdded;
4274: DM newPlexGhosted;
4275: void *ctx;
4277: PetscCall(DMPlexConstructGhostCells(newPlex, pforest->ghostName, &numAdded, &newPlexGhosted));
4278: PetscCall(DMGetApplicationContext(newPlex, &ctx));
4279: PetscCall(DMSetApplicationContext(newPlexGhosted, ctx));
4280: /* we want the sf for the ghost dm to be the one for the p4est dm as well */
4281: PetscCall(DMGetPointSF(newPlexGhosted, &pointSF));
4282: PetscCall(DMSetPointSF(dm, pointSF));
4283: PetscCall(DMDestroy(&newPlex));
4284: PetscCall(DMPlexSetReferenceTree(newPlexGhosted, refTree));
4285: PetscCall(DMForestClearAdaptivityForest_pforest(dm));
4286: newPlex = newPlexGhosted;
4288: /* share the labels back */
4289: PetscCall(DMDestroyLabelLinkList_Internal(dm));
4290: PetscCall(DMCopyLabels(newPlex, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4291: pforest->plex = newPlex;
4292: }
4293: PetscCall(DMDestroy(&refTree));
4294: if (dm->setfromoptionscalled) {
4295: PetscObjectOptionsBegin((PetscObject)newPlex);
4296: PetscCall(DMSetFromOptions_NonRefinement_Plex(newPlex, PetscOptionsObject));
4297: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)newPlex, PetscOptionsObject));
4298: PetscOptionsEnd();
4299: }
4300: PetscCall(DMViewFromOptions(newPlex, NULL, "-dm_p4est_plex_view"));
4301: {
4302: DM cdm;
4303: PetscSection coordsSec;
4304: Vec coords;
4305: PetscInt cDim;
4307: PetscCall(DMGetCoordinateDim(newPlex, &cDim));
4308: PetscCall(DMGetCoordinateSection(newPlex, &coordsSec));
4309: PetscCall(DMSetCoordinateSection(dm, cDim, coordsSec));
4310: PetscCall(DMGetCoordinatesLocal(newPlex, &coords));
4311: PetscCall(DMSetCoordinatesLocal(dm, coords));
4312: PetscCall(DMGetCoordinateDM(newPlex, &cdm));
4313: if (cdm) {
4314: PetscFE fe;
4315: #if !defined(P4_TO_P8)
4316: DMPolytopeType celltype = DM_POLYTOPE_QUADRILATERAL;
4317: #else
4318: DMPolytopeType celltype = DM_POLYTOPE_HEXAHEDRON;
4319: #endif
4321: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dim, celltype, 1, PETSC_DEFAULT, &fe));
4322: PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe));
4323: PetscCall(PetscFEDestroy(&fe));
4324: PetscCall(DMCreateDS(cdm));
4325: }
4326: PetscCall(DMGetCellCoordinateDM(newPlex, &cdm));
4327: if (cdm) PetscCall(DMSetCellCoordinateDM(dm, cdm));
4328: PetscCall(DMGetCellCoordinateSection(newPlex, &coordsSec));
4329: if (coordsSec) PetscCall(DMSetCellCoordinateSection(dm, cDim, coordsSec));
4330: PetscCall(DMGetCellCoordinatesLocal(newPlex, &coords));
4331: if (coords) PetscCall(DMSetCellCoordinatesLocal(dm, coords));
4332: }
4333: } else {
4334: PetscCall(DMCopyLabels(dm, pforest->plex, PETSC_OWN_POINTER, PETSC_FALSE, DM_COPY_LABELS_REPLACE));
4335: }
4336: newPlex = pforest->plex;
4337: if (plex) {
4338: PetscCall(DMClone(newPlex, plex));
4339: #if 0
4340: PetscCall(DMGetCoordinateDM(newPlex,&coordDM));
4341: PetscCall(DMSetCoordinateDM(*plex,coordDM));
4342: PetscCall(DMGetCellCoordinateDM(newPlex,&coordDM));
4343: PetscCall(DMSetCellCoordinateDM(*plex,coordDM));
4344: #endif
4345: PetscCall(DMShareDiscretization(dm, *plex));
4346: }
4347: PetscFunctionReturn(PETSC_SUCCESS);
4348: }
4350: static PetscErrorCode DMSetFromOptions_pforest(DM dm, PetscOptionItems PetscOptionsObject)
4351: {
4352: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4353: char stringBuffer[256];
4354: PetscBool flg;
4356: PetscFunctionBegin;
4357: PetscCall(DMSetFromOptions_Forest(dm, PetscOptionsObject));
4358: PetscOptionsHeadBegin(PetscOptionsObject, "DM" P4EST_STRING " options");
4359: PetscCall(PetscOptionsBool("-dm_p4est_partition_for_coarsening", "partition forest to allow for coarsening", "DMP4estSetPartitionForCoarsening", pforest->partition_for_coarsening, &pforest->partition_for_coarsening, NULL));
4360: PetscCall(PetscOptionsString("-dm_p4est_ghost_label_name", "the name of the ghost label when converting from a DMPlex", NULL, NULL, stringBuffer, sizeof(stringBuffer), &flg));
4361: PetscOptionsHeadEnd();
4362: if (flg) {
4363: PetscCall(PetscFree(pforest->ghostName));
4364: PetscCall(PetscStrallocpy(stringBuffer, &pforest->ghostName));
4365: }
4366: PetscFunctionReturn(PETSC_SUCCESS);
4367: }
4369: #if !defined(P4_TO_P8)
4370: #define DMPforestGetPartitionForCoarsening DMP4estGetPartitionForCoarsening
4371: #define DMPforestSetPartitionForCoarsening DMP4estSetPartitionForCoarsening
4372: #else
4373: #define DMPforestGetPartitionForCoarsening DMP8estGetPartitionForCoarsening
4374: #define DMPforestSetPartitionForCoarsening DMP8estSetPartitionForCoarsening
4375: #endif
4377: PETSC_EXTERN PetscErrorCode DMPforestGetPartitionForCoarsening(DM dm, PetscBool *flg)
4378: {
4379: DM_Forest_pforest *pforest;
4381: PetscFunctionBegin;
4383: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4384: *flg = pforest->partition_for_coarsening;
4385: PetscFunctionReturn(PETSC_SUCCESS);
4386: }
4388: PETSC_EXTERN PetscErrorCode DMPforestSetPartitionForCoarsening(DM dm, PetscBool flg)
4389: {
4390: DM_Forest_pforest *pforest;
4392: PetscFunctionBegin;
4394: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4395: pforest->partition_for_coarsening = flg;
4396: PetscFunctionReturn(PETSC_SUCCESS);
4397: }
4399: static PetscErrorCode DMPforestGetPlex(DM dm, DM *plex)
4400: {
4401: DM_Forest_pforest *pforest;
4403: PetscFunctionBegin;
4404: if (plex) *plex = NULL;
4405: PetscCall(DMSetUp(dm));
4406: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4407: if (!pforest->plex) PetscCall(DMConvert_pforest_plex(dm, DMPLEX, NULL));
4408: PetscCall(DMShareDiscretization(dm, pforest->plex));
4409: if (plex) *plex = pforest->plex;
4410: PetscFunctionReturn(PETSC_SUCCESS);
4411: }
4413: #define DMCreateInterpolation_pforest _append_pforest(DMCreateInterpolation)
4414: static PetscErrorCode DMCreateInterpolation_pforest(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
4415: {
4416: PetscSection gsc, gsf;
4417: PetscInt m, n;
4418: DM cdm;
4420: PetscFunctionBegin;
4421: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4422: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
4423: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4424: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &n));
4426: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), interpolation));
4427: PetscCall(MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4428: PetscCall(MatSetType(*interpolation, MATAIJ));
4430: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4431: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only interpolation from coarse DM for now");
4433: {
4434: DM plexF, plexC;
4435: PetscSF sf;
4436: PetscInt *cids;
4437: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4439: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4440: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4441: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4442: PetscCall(PetscSFSetUp(sf));
4443: PetscCall(DMPlexComputeInterpolatorTree(plexC, plexF, sf, cids, *interpolation));
4444: PetscCall(PetscSFDestroy(&sf));
4445: PetscCall(PetscFree(cids));
4446: }
4447: PetscCall(MatViewFromOptions(*interpolation, NULL, "-interp_mat_view"));
4448: /* Use naive scaling */
4449: PetscCall(DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling));
4450: PetscFunctionReturn(PETSC_SUCCESS);
4451: }
4453: #define DMCreateInjection_pforest _append_pforest(DMCreateInjection)
4454: static PetscErrorCode DMCreateInjection_pforest(DM dmCoarse, DM dmFine, Mat *injection)
4455: {
4456: PetscSection gsc, gsf;
4457: PetscInt m, n;
4458: DM cdm;
4460: PetscFunctionBegin;
4461: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4462: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &n));
4463: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4464: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &m));
4466: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), injection));
4467: PetscCall(MatSetSizes(*injection, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4468: PetscCall(MatSetType(*injection, MATAIJ));
4470: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4471: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only injection to coarse DM for now");
4473: {
4474: DM plexF, plexC;
4475: PetscSF sf;
4476: PetscInt *cids;
4477: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4479: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4480: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4481: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4482: PetscCall(PetscSFSetUp(sf));
4483: PetscCall(DMPlexComputeInjectorTree(plexC, plexF, sf, cids, *injection));
4484: PetscCall(PetscSFDestroy(&sf));
4485: PetscCall(PetscFree(cids));
4486: }
4487: PetscCall(MatViewFromOptions(*injection, NULL, "-inject_mat_view"));
4488: /* Use naive scaling */
4489: PetscFunctionReturn(PETSC_SUCCESS);
4490: }
4492: #define DMForestTransferVecFromBase_pforest _append_pforest(DMForestTransferVecFromBase)
4493: static PetscErrorCode DMForestTransferVecFromBase_pforest(DM dm, Vec vecIn, Vec vecOut)
4494: {
4495: DM dmIn, dmVecIn, base, basec, plex, coarseDM;
4496: DM *hierarchy;
4497: PetscSF sfRed = NULL;
4498: PetscDS ds;
4499: Vec vecInLocal, vecOutLocal;
4500: DMLabel subpointMap;
4501: PetscInt minLevel, mh, n_hi, i;
4502: PetscBool hiforest, *hierarchy_forest;
4504: PetscFunctionBegin;
4505: PetscCall(VecGetDM(vecIn, &dmVecIn));
4506: PetscCall(DMGetDS(dmVecIn, &ds));
4507: PetscCheck(ds, PetscObjectComm((PetscObject)dmVecIn), PETSC_ERR_SUP, "Cannot transfer without a PetscDS object");
4508: { /* we cannot stick user contexts into function callbacks for DMProjectFieldLocal! */
4509: PetscSection section;
4510: PetscInt Nf;
4512: PetscCall(DMGetLocalSection(dmVecIn, §ion));
4513: PetscCall(PetscSectionGetNumFields(section, &Nf));
4514: 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);
4515: }
4516: PetscCall(DMForestGetMinimumRefinement(dm, &minLevel));
4517: PetscCheck(!minLevel, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot transfer with minimum refinement set to %" PetscInt_FMT ". Rerun with DMForestSetMinimumRefinement(dm,0)", minLevel);
4518: PetscCall(DMForestGetBaseDM(dm, &base));
4519: PetscCheck(base, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing base DM");
4521: PetscCall(VecSet(vecOut, 0.0));
4522: if (dmVecIn == base) { /* sequential runs */
4523: PetscCall(PetscObjectReference((PetscObject)vecIn));
4524: } else {
4525: PetscSection secIn, secInRed;
4526: Vec vecInRed, vecInLocal;
4528: PetscCall(PetscObjectQuery((PetscObject)base, "_base_migration_sf", (PetscObject *)&sfRed));
4529: PetscCheck(sfRed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not the DM set with DMForestSetBaseDM()");
4530: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmVecIn), &secInRed));
4531: PetscCall(VecCreate(PETSC_COMM_SELF, &vecInRed));
4532: PetscCall(DMGetLocalSection(dmVecIn, &secIn));
4533: PetscCall(DMGetLocalVector(dmVecIn, &vecInLocal));
4534: PetscCall(DMGlobalToLocalBegin(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4535: PetscCall(DMGlobalToLocalEnd(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4536: PetscCall(DMPlexDistributeField(dmVecIn, sfRed, secIn, vecInLocal, secInRed, vecInRed));
4537: PetscCall(DMRestoreLocalVector(dmVecIn, &vecInLocal));
4538: PetscCall(PetscSectionDestroy(&secInRed));
4539: vecIn = vecInRed;
4540: }
4542: /* we first search through the AdaptivityForest hierarchy
4543: once we found the first disconnected forest, we upsweep the DM hierarchy */
4544: hiforest = PETSC_TRUE;
4546: /* upsweep to the coarsest DM */
4547: n_hi = 0;
4548: coarseDM = dm;
4549: do {
4550: PetscBool isforest;
4552: dmIn = coarseDM;
4553: /* need to call DMSetUp to have the hierarchy recursively setup */
4554: PetscCall(DMSetUp(dmIn));
4555: PetscCall(DMIsForest(dmIn, &isforest));
4556: PetscCheck(isforest, PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Cannot currently transfer through a mixed hierarchy! Found DM type %s", ((PetscObject)dmIn)->type_name);
4557: coarseDM = NULL;
4558: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4559: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4560: hiforest = PETSC_FALSE;
4561: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4562: }
4563: n_hi++;
4564: } while (coarseDM);
4566: PetscCall(PetscMalloc2(n_hi, &hierarchy, n_hi, &hierarchy_forest));
4568: i = 0;
4569: hiforest = PETSC_TRUE;
4570: coarseDM = dm;
4571: do {
4572: dmIn = coarseDM;
4573: coarseDM = NULL;
4574: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4575: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4576: hiforest = PETSC_FALSE;
4577: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4578: }
4579: i++;
4580: hierarchy[n_hi - i] = dmIn;
4581: } while (coarseDM);
4583: /* project base vector on the coarsest forest (minimum refinement = 0) */
4584: PetscCall(DMPforestGetPlex(dmIn, &plex));
4586: /* Check this plex is compatible with the base */
4587: {
4588: IS gnum[2];
4589: PetscInt ncells[2], gncells[2];
4591: PetscCall(DMPlexGetCellNumbering(base, &gnum[0]));
4592: PetscCall(DMPlexGetCellNumbering(plex, &gnum[1]));
4593: PetscCall(ISGetMinMax(gnum[0], NULL, &ncells[0]));
4594: PetscCall(ISGetMinMax(gnum[1], NULL, &ncells[1]));
4595: PetscCallMPI(MPIU_Allreduce(ncells, gncells, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
4596: 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);
4597: }
4599: PetscCall(DMGetLabel(dmIn, "_forest_base_subpoint_map", &subpointMap));
4600: PetscCheck(subpointMap, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing _forest_base_subpoint_map label");
4602: PetscCall(DMPlexGetMaxProjectionHeight(base, &mh));
4603: PetscCall(DMPlexSetMaxProjectionHeight(plex, mh));
4605: PetscCall(DMClone(base, &basec));
4606: PetscCall(DMCopyDisc(dmVecIn, basec));
4607: if (sfRed) {
4608: PetscCall(PetscObjectReference((PetscObject)vecIn));
4609: vecInLocal = vecIn;
4610: } else {
4611: PetscCall(DMCreateLocalVector(basec, &vecInLocal));
4612: PetscCall(DMGlobalToLocalBegin(basec, vecIn, INSERT_VALUES, vecInLocal));
4613: PetscCall(DMGlobalToLocalEnd(basec, vecIn, INSERT_VALUES, vecInLocal));
4614: }
4616: PetscCall(DMGetLocalVector(dmIn, &vecOutLocal));
4617: { /* get degrees of freedom ordered onto dmIn */
4618: PetscSF basetocoarse;
4619: PetscInt bStart, bEnd, nroots;
4620: PetscInt iStart, iEnd, nleaves, leaf;
4621: PetscMPIInt rank;
4622: PetscSFNode *remotes;
4623: PetscSection secIn, secOut;
4624: PetscInt *remoteOffsets;
4625: PetscSF transferSF;
4626: const PetscScalar *inArray;
4627: PetscScalar *outArray;
4629: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)basec), &rank));
4630: PetscCall(DMPlexGetChart(basec, &bStart, &bEnd));
4631: nroots = PetscMax(bEnd - bStart, 0);
4632: PetscCall(DMPlexGetChart(plex, &iStart, &iEnd));
4633: nleaves = PetscMax(iEnd - iStart, 0);
4635: PetscCall(PetscMalloc1(nleaves, &remotes));
4636: for (leaf = iStart; leaf < iEnd; leaf++) {
4637: PetscInt index;
4639: remotes[leaf - iStart].rank = rank;
4640: PetscCall(DMLabelGetValue(subpointMap, leaf, &index));
4641: remotes[leaf - iStart].index = index;
4642: }
4644: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)basec), &basetocoarse));
4645: PetscCall(PetscSFSetGraph(basetocoarse, nroots, nleaves, NULL, PETSC_OWN_POINTER, remotes, PETSC_OWN_POINTER));
4646: PetscCall(PetscSFSetUp(basetocoarse));
4647: PetscCall(DMGetLocalSection(basec, &secIn));
4648: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmIn), &secOut));
4649: PetscCall(PetscSFDistributeSection(basetocoarse, secIn, &remoteOffsets, secOut));
4650: PetscCall(PetscSFCreateSectionSF(basetocoarse, secIn, remoteOffsets, secOut, &transferSF));
4651: PetscCall(PetscFree(remoteOffsets));
4652: PetscCall(VecGetArrayWrite(vecOutLocal, &outArray));
4653: PetscCall(VecGetArrayRead(vecInLocal, &inArray));
4654: PetscCall(PetscSFBcastBegin(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4655: PetscCall(PetscSFBcastEnd(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4656: PetscCall(VecRestoreArrayRead(vecInLocal, &inArray));
4657: PetscCall(VecRestoreArrayWrite(vecOutLocal, &outArray));
4658: PetscCall(PetscSFDestroy(&transferSF));
4659: PetscCall(PetscSectionDestroy(&secOut));
4660: PetscCall(PetscSFDestroy(&basetocoarse));
4661: }
4662: PetscCall(VecDestroy(&vecInLocal));
4663: PetscCall(DMDestroy(&basec));
4664: PetscCall(VecDestroy(&vecIn));
4666: /* output */
4667: if (n_hi > 1) { /* downsweep the stored hierarchy */
4668: Vec vecOut1, vecOut2;
4669: DM fineDM;
4671: PetscCall(DMGetGlobalVector(dmIn, &vecOut1));
4672: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut1));
4673: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4674: for (i = 1; i < n_hi - 1; i++) {
4675: fineDM = hierarchy[i];
4676: PetscCall(DMGetGlobalVector(fineDM, &vecOut2));
4677: PetscCall(DMForestTransferVec(dmIn, vecOut1, fineDM, vecOut2, PETSC_TRUE, 0.0));
4678: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4679: vecOut1 = vecOut2;
4680: dmIn = fineDM;
4681: }
4682: PetscCall(DMForestTransferVec(dmIn, vecOut1, dm, vecOut, PETSC_TRUE, 0.0));
4683: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4684: } else {
4685: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut));
4686: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4687: }
4688: PetscCall(PetscFree2(hierarchy, hierarchy_forest));
4689: PetscFunctionReturn(PETSC_SUCCESS);
4690: }
4692: #define DMForestTransferVec_pforest _append_pforest(DMForestTransferVec)
4693: static PetscErrorCode DMForestTransferVec_pforest(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscBool useBCs, PetscReal time)
4694: {
4695: DM adaptIn, adaptOut, plexIn, plexOut;
4696: DM_Forest *forestIn, *forestOut, *forestAdaptIn, *forestAdaptOut;
4697: PetscInt dofPerDim[] = {1, 1, 1, 1};
4698: PetscSF inSF = NULL, outSF = NULL;
4699: PetscInt *inCids = NULL, *outCids = NULL;
4700: DMAdaptFlag purposeIn, purposeOut;
4702: PetscFunctionBegin;
4703: forestOut = (DM_Forest *)dmOut->data;
4704: forestIn = (DM_Forest *)dmIn->data;
4706: PetscCall(DMForestGetAdaptivityForest(dmOut, &adaptOut));
4707: PetscCall(DMForestGetAdaptivityPurpose(dmOut, &purposeOut));
4708: forestAdaptOut = adaptOut ? (DM_Forest *)adaptOut->data : NULL;
4710: PetscCall(DMForestGetAdaptivityForest(dmIn, &adaptIn));
4711: PetscCall(DMForestGetAdaptivityPurpose(dmIn, &purposeIn));
4712: forestAdaptIn = adaptIn ? (DM_Forest *)adaptIn->data : NULL;
4714: if (forestAdaptOut == forestIn) {
4715: switch (purposeOut) {
4716: case DM_ADAPT_REFINE:
4717: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4718: PetscCall(PetscSFSetUp(inSF));
4719: break;
4720: case DM_ADAPT_COARSEN:
4721: case DM_ADAPT_COARSEN_LAST:
4722: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &outCids));
4723: PetscCall(PetscSFSetUp(outSF));
4724: break;
4725: default:
4726: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4727: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4728: PetscCall(PetscSFSetUp(inSF));
4729: PetscCall(PetscSFSetUp(outSF));
4730: }
4731: } else if (forestAdaptIn == forestOut) {
4732: switch (purposeIn) {
4733: case DM_ADAPT_REFINE:
4734: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &inCids));
4735: PetscCall(PetscSFSetUp(outSF));
4736: break;
4737: case DM_ADAPT_COARSEN:
4738: case DM_ADAPT_COARSEN_LAST:
4739: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4740: PetscCall(PetscSFSetUp(inSF));
4741: break;
4742: default:
4743: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4744: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4745: PetscCall(PetscSFSetUp(inSF));
4746: PetscCall(PetscSFSetUp(outSF));
4747: }
4748: } else SETERRQ(PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Only support transfer from pre-adaptivity to post-adaptivity right now");
4749: PetscCall(DMPforestGetPlex(dmIn, &plexIn));
4750: PetscCall(DMPforestGetPlex(dmOut, &plexOut));
4752: PetscCall(DMPlexTransferVecTree(plexIn, vecIn, plexOut, vecOut, inSF, outSF, inCids, outCids, useBCs, time));
4753: PetscCall(PetscFree(inCids));
4754: PetscCall(PetscFree(outCids));
4755: PetscCall(PetscSFDestroy(&inSF));
4756: PetscCall(PetscSFDestroy(&outSF));
4757: PetscCall(PetscFree(inCids));
4758: PetscCall(PetscFree(outCids));
4759: PetscFunctionReturn(PETSC_SUCCESS);
4760: }
4762: #define DMCreateCoordinateDM_pforest _append_pforest(DMCreateCoordinateDM)
4763: static PetscErrorCode DMCreateCoordinateDM_pforest(DM dm, DM *cdm)
4764: {
4765: DM plex;
4767: PetscFunctionBegin;
4769: PetscCall(DMPforestGetPlex(dm, &plex));
4770: PetscCall(DMGetCoordinateDM(plex, cdm));
4771: PetscCall(PetscObjectReference((PetscObject)*cdm));
4772: PetscFunctionReturn(PETSC_SUCCESS);
4773: }
4775: #define VecViewLocal_pforest _append_pforest(VecViewLocal)
4776: static PetscErrorCode VecViewLocal_pforest(Vec vec, PetscViewer viewer)
4777: {
4778: DM dm, plex;
4780: PetscFunctionBegin;
4781: PetscCall(VecGetDM(vec, &dm));
4782: PetscCall(PetscObjectReference((PetscObject)dm));
4783: PetscCall(DMPforestGetPlex(dm, &plex));
4784: PetscCall(VecSetDM(vec, plex));
4785: PetscCall(VecView_Plex_Local(vec, viewer));
4786: PetscCall(VecSetDM(vec, dm));
4787: PetscCall(DMDestroy(&dm));
4788: PetscFunctionReturn(PETSC_SUCCESS);
4789: }
4791: #define VecView_pforest _append_pforest(VecView)
4792: static PetscErrorCode VecView_pforest(Vec vec, PetscViewer viewer)
4793: {
4794: DM dm, plex;
4796: PetscFunctionBegin;
4797: PetscCall(VecGetDM(vec, &dm));
4798: PetscCall(PetscObjectReference((PetscObject)dm));
4799: PetscCall(DMPforestGetPlex(dm, &plex));
4800: PetscCall(VecSetDM(vec, plex));
4801: PetscCall(VecView_Plex(vec, viewer));
4802: PetscCall(VecSetDM(vec, dm));
4803: PetscCall(DMDestroy(&dm));
4804: PetscFunctionReturn(PETSC_SUCCESS);
4805: }
4807: #define VecView_pforest_Native _infix_pforest(VecView, _Native)
4808: static PetscErrorCode VecView_pforest_Native(Vec vec, PetscViewer viewer)
4809: {
4810: DM dm, plex;
4812: PetscFunctionBegin;
4813: PetscCall(VecGetDM(vec, &dm));
4814: PetscCall(PetscObjectReference((PetscObject)dm));
4815: PetscCall(DMPforestGetPlex(dm, &plex));
4816: PetscCall(VecSetDM(vec, plex));
4817: PetscCall(VecView_Plex_Native(vec, viewer));
4818: PetscCall(VecSetDM(vec, dm));
4819: PetscCall(DMDestroy(&dm));
4820: PetscFunctionReturn(PETSC_SUCCESS);
4821: }
4823: #define VecLoad_pforest _append_pforest(VecLoad)
4824: static PetscErrorCode VecLoad_pforest(Vec vec, PetscViewer viewer)
4825: {
4826: DM dm, plex;
4828: PetscFunctionBegin;
4829: PetscCall(VecGetDM(vec, &dm));
4830: PetscCall(PetscObjectReference((PetscObject)dm));
4831: PetscCall(DMPforestGetPlex(dm, &plex));
4832: PetscCall(VecSetDM(vec, plex));
4833: PetscCall(VecLoad_Plex(vec, viewer));
4834: PetscCall(VecSetDM(vec, dm));
4835: PetscCall(DMDestroy(&dm));
4836: PetscFunctionReturn(PETSC_SUCCESS);
4837: }
4839: #define VecLoad_pforest_Native _infix_pforest(VecLoad, _Native)
4840: static PetscErrorCode VecLoad_pforest_Native(Vec vec, PetscViewer viewer)
4841: {
4842: DM dm, plex;
4844: PetscFunctionBegin;
4845: PetscCall(VecGetDM(vec, &dm));
4846: PetscCall(PetscObjectReference((PetscObject)dm));
4847: PetscCall(DMPforestGetPlex(dm, &plex));
4848: PetscCall(VecSetDM(vec, plex));
4849: PetscCall(VecLoad_Plex_Native(vec, viewer));
4850: PetscCall(VecSetDM(vec, dm));
4851: PetscCall(DMDestroy(&dm));
4852: PetscFunctionReturn(PETSC_SUCCESS);
4853: }
4855: #define DMCreateGlobalVector_pforest _append_pforest(DMCreateGlobalVector)
4856: static PetscErrorCode DMCreateGlobalVector_pforest(DM dm, Vec *vec)
4857: {
4858: PetscFunctionBegin;
4859: PetscCall(DMCreateGlobalVector_Section_Private(dm, vec));
4860: /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */
4861: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_pforest));
4862: PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void))VecView_pforest_Native));
4863: PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_pforest));
4864: PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void))VecLoad_pforest_Native));
4865: PetscFunctionReturn(PETSC_SUCCESS);
4866: }
4868: #define DMCreateLocalVector_pforest _append_pforest(DMCreateLocalVector)
4869: static PetscErrorCode DMCreateLocalVector_pforest(DM dm, Vec *vec)
4870: {
4871: PetscFunctionBegin;
4872: PetscCall(DMCreateLocalVector_Section_Private(dm, vec));
4873: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecViewLocal_pforest));
4874: PetscFunctionReturn(PETSC_SUCCESS);
4875: }
4877: #define DMCreateMatrix_pforest _append_pforest(DMCreateMatrix)
4878: static PetscErrorCode DMCreateMatrix_pforest(DM dm, Mat *mat)
4879: {
4880: DM plex;
4882: PetscFunctionBegin;
4884: PetscCall(DMPforestGetPlex(dm, &plex));
4885: if (plex->prealloc_only != dm->prealloc_only) plex->prealloc_only = dm->prealloc_only; /* maybe this should go into forest->plex */
4886: PetscCall(DMSetMatType(plex, dm->mattype));
4887: PetscCall(DMCreateMatrix(plex, mat));
4888: PetscCall(MatSetDM(*mat, dm));
4889: PetscFunctionReturn(PETSC_SUCCESS);
4890: }
4892: #define DMProjectFunctionLocal_pforest _append_pforest(DMProjectFunctionLocal)
4893: static PetscErrorCode DMProjectFunctionLocal_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4894: {
4895: DM plex;
4897: PetscFunctionBegin;
4899: PetscCall(DMPforestGetPlex(dm, &plex));
4900: PetscCall(DMProjectFunctionLocal(plex, time, funcs, ctxs, mode, localX));
4901: PetscFunctionReturn(PETSC_SUCCESS);
4902: }
4904: #define DMProjectFunctionLabelLocal_pforest _append_pforest(DMProjectFunctionLabelLocal)
4905: 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)
4906: {
4907: DM plex;
4909: PetscFunctionBegin;
4911: PetscCall(DMPforestGetPlex(dm, &plex));
4912: PetscCall(DMProjectFunctionLabelLocal(plex, time, label, numIds, ids, Ncc, comps, funcs, ctxs, mode, localX));
4913: PetscFunctionReturn(PETSC_SUCCESS);
4914: }
4916: #define DMProjectFieldLocal_pforest _append_pforest(DMProjectFieldLocal)
4917: 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)
4918: {
4919: DM plex;
4921: PetscFunctionBegin;
4923: PetscCall(DMPforestGetPlex(dm, &plex));
4924: PetscCall(DMProjectFieldLocal(plex, time, localU, funcs, mode, localX));
4925: PetscFunctionReturn(PETSC_SUCCESS);
4926: }
4928: #define DMComputeL2Diff_pforest _append_pforest(DMComputeL2Diff)
4929: PetscErrorCode DMComputeL2Diff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
4930: {
4931: DM plex;
4933: PetscFunctionBegin;
4935: PetscCall(DMPforestGetPlex(dm, &plex));
4936: PetscCall(DMComputeL2Diff(plex, time, funcs, ctxs, X, diff));
4937: PetscFunctionReturn(PETSC_SUCCESS);
4938: }
4940: #define DMComputeL2FieldDiff_pforest _append_pforest(DMComputeL2FieldDiff)
4941: PetscErrorCode DMComputeL2FieldDiff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
4942: {
4943: DM plex;
4945: PetscFunctionBegin;
4947: PetscCall(DMPforestGetPlex(dm, &plex));
4948: PetscCall(DMComputeL2FieldDiff(plex, time, funcs, ctxs, X, diff));
4949: PetscFunctionReturn(PETSC_SUCCESS);
4950: }
4952: #define DMCreatelocalsection_pforest _append_pforest(DMCreatelocalsection)
4953: static PetscErrorCode DMCreatelocalsection_pforest(DM dm)
4954: {
4955: DM plex;
4956: PetscSection section;
4958: PetscFunctionBegin;
4960: PetscCall(DMPforestGetPlex(dm, &plex));
4961: PetscCall(DMGetLocalSection(plex, §ion));
4962: PetscCall(DMSetLocalSection(dm, section));
4963: PetscFunctionReturn(PETSC_SUCCESS);
4964: }
4966: #define DMCreateDefaultConstraints_pforest _append_pforest(DMCreateDefaultConstraints)
4967: static PetscErrorCode DMCreateDefaultConstraints_pforest(DM dm)
4968: {
4969: DM plex;
4970: Mat mat;
4971: Vec bias;
4972: PetscSection section;
4974: PetscFunctionBegin;
4976: PetscCall(DMPforestGetPlex(dm, &plex));
4977: PetscCall(DMGetDefaultConstraints(plex, §ion, &mat, &bias));
4978: PetscCall(DMSetDefaultConstraints(dm, section, mat, bias));
4979: PetscFunctionReturn(PETSC_SUCCESS);
4980: }
4982: #define DMGetDimPoints_pforest _append_pforest(DMGetDimPoints)
4983: static PetscErrorCode DMGetDimPoints_pforest(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
4984: {
4985: DM plex;
4987: PetscFunctionBegin;
4989: PetscCall(DMPforestGetPlex(dm, &plex));
4990: PetscCall(DMGetDimPoints(plex, dim, cStart, cEnd));
4991: PetscFunctionReturn(PETSC_SUCCESS);
4992: }
4994: /* Need to forward declare */
4995: #define DMInitialize_pforest _append_pforest(DMInitialize)
4996: static PetscErrorCode DMInitialize_pforest(DM dm);
4998: #define DMClone_pforest _append_pforest(DMClone)
4999: static PetscErrorCode DMClone_pforest(DM dm, DM *newdm)
5000: {
5001: PetscFunctionBegin;
5002: PetscCall(DMClone_Forest(dm, newdm));
5003: PetscCall(DMInitialize_pforest(*newdm));
5004: PetscFunctionReturn(PETSC_SUCCESS);
5005: }
5007: #define DMForestCreateCellChart_pforest _append_pforest(DMForestCreateCellChart)
5008: static PetscErrorCode DMForestCreateCellChart_pforest(DM dm, PetscInt *cStart, PetscInt *cEnd)
5009: {
5010: DM_Forest *forest;
5011: DM_Forest_pforest *pforest;
5012: PetscInt overlap;
5014: PetscFunctionBegin;
5015: PetscCall(DMSetUp(dm));
5016: forest = (DM_Forest *)dm->data;
5017: pforest = (DM_Forest_pforest *)forest->data;
5018: *cStart = 0;
5019: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
5020: if (overlap && pforest->ghost) {
5021: *cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[pforest->forest->mpisize];
5022: } else {
5023: *cEnd = pforest->forest->local_num_quadrants;
5024: }
5025: PetscFunctionReturn(PETSC_SUCCESS);
5026: }
5028: #define DMForestCreateCellSF_pforest _append_pforest(DMForestCreateCellSF)
5029: static PetscErrorCode DMForestCreateCellSF_pforest(DM dm, PetscSF *cellSF)
5030: {
5031: DM_Forest *forest;
5032: DM_Forest_pforest *pforest;
5033: PetscMPIInt rank;
5034: PetscInt overlap;
5035: PetscInt cStart, cEnd, cLocalStart, cLocalEnd;
5036: PetscInt nRoots, nLeaves, *mine = NULL;
5037: PetscSFNode *remote = NULL;
5038: PetscSF sf;
5040: PetscFunctionBegin;
5041: PetscCall(DMForestGetCellChart(dm, &cStart, &cEnd));
5042: forest = (DM_Forest *)dm->data;
5043: pforest = (DM_Forest_pforest *)forest->data;
5044: nRoots = cEnd - cStart;
5045: cLocalStart = pforest->cLocalStart;
5046: cLocalEnd = pforest->cLocalEnd;
5047: nLeaves = 0;
5048: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
5049: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
5050: if (overlap && pforest->ghost) {
5051: PetscSFNode *mirror;
5052: p4est_quadrant_t *mirror_array;
5053: PetscInt nMirror, nGhostPre, nSelf, q;
5054: void **mirrorPtrs;
5056: nMirror = (PetscInt)pforest->ghost->mirrors.elem_count;
5057: nSelf = cLocalEnd - cLocalStart;
5058: nLeaves = nRoots - nSelf;
5059: nGhostPre = (PetscInt)pforest->ghost->proc_offsets[rank];
5060: PetscCall(PetscMalloc1(nLeaves, &mine));
5061: PetscCall(PetscMalloc1(nLeaves, &remote));
5062: PetscCall(PetscMalloc2(nMirror, &mirror, nMirror, &mirrorPtrs));
5063: mirror_array = (p4est_quadrant_t *)pforest->ghost->mirrors.array;
5064: for (q = 0; q < nMirror; q++) {
5065: p4est_quadrant_t *mir = &mirror_array[q];
5067: mirror[q].rank = rank;
5068: mirror[q].index = (PetscInt)mir->p.piggy3.local_num + cLocalStart;
5069: mirrorPtrs[q] = (void *)&mirror[q];
5070: }
5071: PetscCallP4est(p4est_ghost_exchange_custom, (pforest->forest, pforest->ghost, sizeof(PetscSFNode), mirrorPtrs, remote));
5072: PetscCall(PetscFree2(mirror, mirrorPtrs));
5073: for (q = 0; q < nGhostPre; q++) mine[q] = q;
5074: for (; q < nLeaves; q++) mine[q] = (q - nGhostPre) + cLocalEnd;
5075: }
5076: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf));
5077: PetscCall(PetscSFSetGraph(sf, nRoots, nLeaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
5078: *cellSF = sf;
5079: PetscFunctionReturn(PETSC_SUCCESS);
5080: }
5082: static PetscErrorCode DMCreateNeumannOverlap_pforest(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx)
5083: {
5084: DM plex;
5086: PetscFunctionBegin;
5087: PetscCall(DMPforestGetPlex(dm, &plex));
5088: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5089: PetscCall(DMCreateNeumannOverlap_Plex(plex, ovl, J, setup, setup_ctx));
5090: PetscCall(DMClearAuxiliaryVec(plex));
5091: if (!*setup) {
5092: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup));
5093: if (*setup) PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm));
5094: }
5095: PetscFunctionReturn(PETSC_SUCCESS);
5096: }
5098: #define DMCreateDomainDecomposition_pforest _append_pforest(DMCreateDomainDecomposition)
5099: static PetscErrorCode DMCreateDomainDecomposition_pforest(DM dm, PetscInt *nsub, char ***names, IS **innerises, IS **outerises, DM **dms)
5100: {
5101: DM plex;
5103: PetscFunctionBegin;
5104: PetscCall(DMPforestGetPlex(dm, &plex));
5105: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5106: PetscCall(DMCreateDomainDecomposition(plex, nsub, names, innerises, outerises, dms));
5107: PetscCall(DMClearAuxiliaryVec(plex));
5108: PetscFunctionReturn(PETSC_SUCCESS);
5109: }
5111: #define DMCreateDomainDecompositionScatters_pforest _append_pforest(DMCreateDomainDecompositionScatters)
5112: static PetscErrorCode DMCreateDomainDecompositionScatters_pforest(DM dm, PetscInt n, DM *subdms, VecScatter **iscat, VecScatter **oscat, VecScatter **lscat)
5113: {
5114: DM plex;
5116: PetscFunctionBegin;
5117: PetscCall(DMPforestGetPlex(dm, &plex));
5118: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5119: PetscCall(DMCreateDomainDecompositionScatters(plex, n, subdms, iscat, oscat, lscat));
5120: PetscFunctionReturn(PETSC_SUCCESS);
5121: }
5123: static PetscErrorCode DMInitialize_pforest(DM dm)
5124: {
5125: PetscFunctionBegin;
5126: dm->ops->setup = DMSetUp_pforest;
5127: dm->ops->view = DMView_pforest;
5128: dm->ops->clone = DMClone_pforest;
5129: dm->ops->createinterpolation = DMCreateInterpolation_pforest;
5130: dm->ops->createinjection = DMCreateInjection_pforest;
5131: dm->ops->setfromoptions = DMSetFromOptions_pforest;
5132: dm->ops->createcoordinatedm = DMCreateCoordinateDM_pforest;
5133: dm->ops->createcellcoordinatedm = NULL;
5134: dm->ops->createglobalvector = DMCreateGlobalVector_pforest;
5135: dm->ops->createlocalvector = DMCreateLocalVector_pforest;
5136: dm->ops->creatematrix = DMCreateMatrix_pforest;
5137: dm->ops->projectfunctionlocal = DMProjectFunctionLocal_pforest;
5138: dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_pforest;
5139: dm->ops->projectfieldlocal = DMProjectFieldLocal_pforest;
5140: dm->ops->createlocalsection = DMCreatelocalsection_pforest;
5141: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_pforest;
5142: dm->ops->computel2diff = DMComputeL2Diff_pforest;
5143: dm->ops->computel2fielddiff = DMComputeL2FieldDiff_pforest;
5144: dm->ops->getdimpoints = DMGetDimPoints_pforest;
5145: dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_pforest;
5146: dm->ops->createddscatters = DMCreateDomainDecompositionScatters_pforest;
5148: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", DMConvert_plex_pforest));
5149: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", DMConvert_pforest_plex));
5150: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_pforest));
5151: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMForestGetPartitionOverlap));
5152: PetscFunctionReturn(PETSC_SUCCESS);
5153: }
5155: #define DMCreate_pforest _append_pforest(DMCreate)
5156: PETSC_EXTERN PetscErrorCode DMCreate_pforest(DM dm)
5157: {
5158: DM_Forest *forest;
5159: DM_Forest_pforest *pforest;
5161: PetscFunctionBegin;
5162: PetscCall(PetscP4estInitialize());
5163: PetscCall(DMCreate_Forest(dm));
5164: PetscCall(DMInitialize_pforest(dm));
5165: PetscCall(DMSetDimension(dm, P4EST_DIM));
5167: /* set forest defaults */
5168: PetscCall(DMForestSetTopology(dm, "unit"));
5169: PetscCall(DMForestSetMinimumRefinement(dm, 0));
5170: PetscCall(DMForestSetInitialRefinement(dm, 0));
5171: PetscCall(DMForestSetMaximumRefinement(dm, P4EST_QMAXLEVEL));
5172: PetscCall(DMForestSetGradeFactor(dm, 2));
5173: PetscCall(DMForestSetAdjacencyDimension(dm, 0));
5174: PetscCall(DMForestSetPartitionOverlap(dm, 0));
5176: /* create p4est data */
5177: PetscCall(PetscNew(&pforest));
5179: forest = (DM_Forest *)dm->data;
5180: forest->data = pforest;
5181: forest->destroy = DMForestDestroy_pforest;
5182: forest->ftemplate = DMForestTemplate_pforest;
5183: forest->transfervec = DMForestTransferVec_pforest;
5184: forest->transfervecfrombase = DMForestTransferVecFromBase_pforest;
5185: forest->createcellchart = DMForestCreateCellChart_pforest;
5186: forest->createcellsf = DMForestCreateCellSF_pforest;
5187: forest->clearadaptivityforest = DMForestClearAdaptivityForest_pforest;
5188: forest->getadaptivitysuccess = DMForestGetAdaptivitySuccess_pforest;
5189: pforest->topo = NULL;
5190: pforest->forest = NULL;
5191: pforest->ghost = NULL;
5192: pforest->lnodes = NULL;
5193: pforest->partition_for_coarsening = PETSC_TRUE;
5194: pforest->coarsen_hierarchy = PETSC_FALSE;
5195: pforest->cLocalStart = -1;
5196: pforest->cLocalEnd = -1;
5197: pforest->labelsFinalized = PETSC_FALSE;
5198: pforest->ghostName = NULL;
5199: PetscFunctionReturn(PETSC_SUCCESS);
5200: }
5202: #endif /* defined(PETSC_HAVE_P4EST) */