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: PetscCall(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_MAX_INT;
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: PetscCall(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: if (!forest->cellWeights && forest->weightCapacity == 1. && forest->weightsFactor == 1.) {
1117: PetscCallP4estReturn(shipped, p4est_partition_ext, (pforest->forest, (int)pforest->partition_for_coarsening, NULL));
1118: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Non-uniform partition cases not implemented yet");
1119: if (shipped) ctx.anyChange = PETSC_TRUE;
1120: if (forest_copy) {
1121: if (preCoarseToFine || coarseToPreFine) {
1122: PetscSF repartSF; /* repartSF has roots in the old partition */
1123: PetscInt pStart = -1, pEnd = -1, p;
1124: PetscInt numRoots, numLeaves;
1125: PetscSFNode *repartRoots;
1126: p4est_gloidx_t postStart = pforest->forest->global_first_quadrant[rank];
1127: p4est_gloidx_t postEnd = pforest->forest->global_first_quadrant[rank + 1];
1128: p4est_gloidx_t partOffset = postStart;
1130: numRoots = (PetscInt)(forest_copy->global_first_quadrant[rank + 1] - forest_copy->global_first_quadrant[rank]);
1131: numLeaves = (PetscInt)(postEnd - postStart);
1132: PetscCall(DMPforestComputeOverlappingRanks(size, rank, pforest->forest, forest_copy, &pStart, &pEnd));
1133: PetscCall(PetscMalloc1((PetscInt)pforest->forest->local_num_quadrants, &repartRoots));
1134: for (p = pStart; p < pEnd; p++) {
1135: p4est_gloidx_t preStart = forest_copy->global_first_quadrant[p];
1136: p4est_gloidx_t preEnd = forest_copy->global_first_quadrant[p + 1];
1137: PetscInt q;
1139: if (preEnd == preStart) continue;
1140: PetscCheck(preStart <= postStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad partition overlap computation");
1141: preEnd = preEnd > postEnd ? postEnd : preEnd;
1142: for (q = partOffset; q < preEnd; q++) {
1143: repartRoots[q - postStart].rank = p;
1144: repartRoots[q - postStart].index = partOffset - preStart;
1145: }
1146: partOffset = preEnd;
1147: }
1148: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &repartSF));
1149: PetscCall(PetscSFSetGraph(repartSF, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, repartRoots, PETSC_OWN_POINTER));
1150: PetscCall(PetscSFSetUp(repartSF));
1151: if (preCoarseToFine) {
1152: PetscSF repartSFembed, preCoarseToFineNew;
1153: PetscInt nleaves;
1154: const PetscInt *leaves;
1156: PetscCall(PetscSFSetUp(preCoarseToFine));
1157: PetscCall(PetscSFGetGraph(preCoarseToFine, NULL, &nleaves, &leaves, NULL));
1158: if (leaves) {
1159: PetscCall(PetscSFCreateEmbeddedRootSF(repartSF, nleaves, leaves, &repartSFembed));
1160: } else {
1161: repartSFembed = repartSF;
1162: PetscCall(PetscObjectReference((PetscObject)repartSFembed));
1163: }
1164: PetscCall(PetscSFCompose(preCoarseToFine, repartSFembed, &preCoarseToFineNew));
1165: PetscCall(PetscSFDestroy(&preCoarseToFine));
1166: PetscCall(PetscSFDestroy(&repartSFembed));
1167: preCoarseToFine = preCoarseToFineNew;
1168: }
1169: if (coarseToPreFine) {
1170: PetscSF repartSFinv, coarseToPreFineNew;
1172: PetscCall(PetscSFCreateInverseSF(repartSF, &repartSFinv));
1173: PetscCall(PetscSFCompose(repartSFinv, coarseToPreFine, &coarseToPreFineNew));
1174: PetscCall(PetscSFDestroy(&coarseToPreFine));
1175: PetscCall(PetscSFDestroy(&repartSFinv));
1176: coarseToPreFine = coarseToPreFineNew;
1177: }
1178: PetscCall(PetscSFDestroy(&repartSF));
1179: }
1180: PetscCallP4est(p4est_destroy, (forest_copy));
1181: }
1182: }
1183: if (size > 1) {
1184: PetscInt overlap;
1186: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
1188: if (adaptFrom) {
1189: PetscInt aoverlap;
1191: PetscCall(DMForestGetPartitionOverlap(adaptFrom, &aoverlap));
1192: if (aoverlap != overlap) ctx.anyChange = PETSC_TRUE;
1193: }
1195: if (overlap > 0) {
1196: PetscInt i, cLocalStart;
1197: PetscInt cEnd;
1198: PetscSF preCellSF = NULL, cellSF = NULL;
1200: PetscCallP4estReturn(pforest->ghost, p4est_ghost_new, (pforest->forest, P4EST_CONNECT_FULL));
1201: PetscCallP4estReturn(pforest->lnodes, p4est_lnodes_new, (pforest->forest, pforest->ghost, -P4EST_DIM));
1202: PetscCallP4est(p4est_ghost_support_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));
1203: for (i = 1; i < overlap; i++) PetscCallP4est(p4est_ghost_expand_by_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost));
1205: cLocalStart = pforest->cLocalStart = pforest->ghost->proc_offsets[rank];
1206: cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[size];
1208: /* shift sfs by cLocalStart, expand by cell SFs */
1209: if (preCoarseToFine || coarseToPreFine) {
1210: if (adaptFrom) PetscCall(DMForestGetCellSF(adaptFrom, &preCellSF));
1211: dm->setupcalled = PETSC_TRUE;
1212: PetscCall(DMForestGetCellSF(dm, &cellSF));
1213: }
1214: if (preCoarseToFine) {
1215: PetscSF preCoarseToFineNew;
1216: PetscInt nleaves, nroots, *leavesNew, i, nleavesNew;
1217: const PetscInt *leaves;
1218: const PetscSFNode *remotes;
1219: PetscSFNode *remotesAll;
1221: PetscCall(PetscSFSetUp(preCoarseToFine));
1222: PetscCall(PetscSFGetGraph(preCoarseToFine, &nroots, &nleaves, &leaves, &remotes));
1223: PetscCall(PetscMalloc1(cEnd, &remotesAll));
1224: for (i = 0; i < cEnd; i++) {
1225: remotesAll[i].rank = -1;
1226: remotesAll[i].index = -1;
1227: }
1228: for (i = 0; i < nleaves; i++) remotesAll[(leaves ? leaves[i] : i) + cLocalStart] = remotes[i];
1229: PetscCall(PetscSFSetUp(cellSF));
1230: PetscCall(PetscSFBcastBegin(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE));
1231: PetscCall(PetscSFBcastEnd(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE));
1232: nleavesNew = 0;
1233: for (i = 0; i < nleaves; i++) {
1234: if (remotesAll[i].rank >= 0) nleavesNew++;
1235: }
1236: PetscCall(PetscMalloc1(nleavesNew, &leavesNew));
1237: nleavesNew = 0;
1238: for (i = 0; i < nleaves; i++) {
1239: if (remotesAll[i].rank >= 0) {
1240: leavesNew[nleavesNew] = i;
1241: if (i > nleavesNew) remotesAll[nleavesNew] = remotesAll[i];
1242: nleavesNew++;
1243: }
1244: }
1245: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &preCoarseToFineNew));
1246: if (nleavesNew < cEnd) {
1247: PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, leavesNew, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES));
1248: } else { /* all cells are leaves */
1249: PetscCall(PetscFree(leavesNew));
1250: PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, NULL, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES));
1251: }
1252: PetscCall(PetscFree(remotesAll));
1253: PetscCall(PetscSFDestroy(&preCoarseToFine));
1254: preCoarseToFine = preCoarseToFineNew;
1255: preCoarseToFine = preCoarseToFineNew;
1256: }
1257: if (coarseToPreFine) {
1258: PetscSF coarseToPreFineNew;
1259: PetscInt nleaves, nroots, i, nleavesCellSF, nleavesExpanded, *leavesNew;
1260: const PetscInt *leaves;
1261: const PetscSFNode *remotes;
1262: PetscSFNode *remotesNew, *remotesNewRoot, *remotesExpanded;
1264: PetscCall(PetscSFSetUp(coarseToPreFine));
1265: PetscCall(PetscSFGetGraph(coarseToPreFine, &nroots, &nleaves, &leaves, &remotes));
1266: PetscCall(PetscSFGetGraph(preCellSF, NULL, &nleavesCellSF, NULL, NULL));
1267: PetscCall(PetscMalloc1(nroots, &remotesNewRoot));
1268: PetscCall(PetscMalloc1(nleaves, &remotesNew));
1269: for (i = 0; i < nroots; i++) {
1270: remotesNewRoot[i].rank = rank;
1271: remotesNewRoot[i].index = i + cLocalStart;
1272: }
1273: PetscCall(PetscSFBcastBegin(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE));
1274: PetscCall(PetscSFBcastEnd(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE));
1275: PetscCall(PetscFree(remotesNewRoot));
1276: PetscCall(PetscMalloc1(nleavesCellSF, &remotesExpanded));
1277: for (i = 0; i < nleavesCellSF; i++) {
1278: remotesExpanded[i].rank = -1;
1279: remotesExpanded[i].index = -1;
1280: }
1281: for (i = 0; i < nleaves; i++) remotesExpanded[leaves ? leaves[i] : i] = remotesNew[i];
1282: PetscCall(PetscFree(remotesNew));
1283: PetscCall(PetscSFBcastBegin(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE));
1284: PetscCall(PetscSFBcastEnd(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE));
1286: nleavesExpanded = 0;
1287: for (i = 0; i < nleavesCellSF; i++) {
1288: if (remotesExpanded[i].rank >= 0) nleavesExpanded++;
1289: }
1290: PetscCall(PetscMalloc1(nleavesExpanded, &leavesNew));
1291: nleavesExpanded = 0;
1292: for (i = 0; i < nleavesCellSF; i++) {
1293: if (remotesExpanded[i].rank >= 0) {
1294: leavesNew[nleavesExpanded] = i;
1295: if (i > nleavesExpanded) remotesExpanded[nleavesExpanded] = remotes[i];
1296: nleavesExpanded++;
1297: }
1298: }
1299: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &coarseToPreFineNew));
1300: if (nleavesExpanded < nleavesCellSF) {
1301: PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, leavesNew, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES));
1302: } else {
1303: PetscCall(PetscFree(leavesNew));
1304: PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, NULL, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES));
1305: }
1306: PetscCall(PetscFree(remotesExpanded));
1307: PetscCall(PetscSFDestroy(&coarseToPreFine));
1308: coarseToPreFine = coarseToPreFineNew;
1309: }
1310: }
1311: }
1312: }
1313: forest->preCoarseToFine = preCoarseToFine;
1314: forest->coarseToPreFine = coarseToPreFine;
1315: dm->setupcalled = PETSC_TRUE;
1316: PetscCall(MPIU_Allreduce(&ctx.anyChange, &pforest->adaptivitySuccess, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
1317: PetscCall(DMPforestGetPlex(dm, NULL));
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: #define DMForestGetAdaptivitySuccess_pforest _append_pforest(DMForestGetAdaptivitySuccess)
1322: static PetscErrorCode DMForestGetAdaptivitySuccess_pforest(DM dm, PetscBool *success)
1323: {
1324: DM_Forest *forest;
1325: DM_Forest_pforest *pforest;
1327: PetscFunctionBegin;
1328: forest = (DM_Forest *)dm->data;
1329: pforest = (DM_Forest_pforest *)forest->data;
1330: *success = pforest->adaptivitySuccess;
1331: PetscFunctionReturn(PETSC_SUCCESS);
1332: }
1334: #define DMView_ASCII_pforest _append_pforest(DMView_ASCII)
1335: static PetscErrorCode DMView_ASCII_pforest(PetscObject odm, PetscViewer viewer)
1336: {
1337: DM dm = (DM)odm;
1339: PetscFunctionBegin;
1342: PetscCall(DMSetUp(dm));
1343: switch (viewer->format) {
1344: case PETSC_VIEWER_DEFAULT:
1345: case PETSC_VIEWER_ASCII_INFO: {
1346: PetscInt dim;
1347: const char *name;
1349: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1350: PetscCall(DMGetDimension(dm, &dim));
1351: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "Forest %s in %" PetscInt_FMT " dimensions:\n", name, dim));
1352: else PetscCall(PetscViewerASCIIPrintf(viewer, "Forest in %" PetscInt_FMT " dimensions:\n", dim));
1353: } /* fall through */
1354: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1355: case PETSC_VIEWER_LOAD_BALANCE: {
1356: DM plex;
1358: PetscCall(DMPforestGetPlex(dm, &plex));
1359: PetscCall(DMView(plex, viewer));
1360: } break;
1361: default:
1362: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1363: }
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: #define DMView_VTK_pforest _append_pforest(DMView_VTK)
1368: static PetscErrorCode DMView_VTK_pforest(PetscObject odm, PetscViewer viewer)
1369: {
1370: DM dm = (DM)odm;
1371: DM_Forest *forest = (DM_Forest *)dm->data;
1372: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
1373: PetscBool isvtk;
1374: PetscReal vtkScale = 1. - PETSC_MACHINE_EPSILON;
1375: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
1376: const char *name;
1377: char *filenameStrip = NULL;
1378: PetscBool hasExt;
1379: size_t len;
1380: p4est_geometry_t *geom;
1382: PetscFunctionBegin;
1385: PetscCall(DMSetUp(dm));
1386: geom = pforest->topo->geom;
1387: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1388: PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
1389: switch (viewer->format) {
1390: case PETSC_VIEWER_VTK_VTU:
1391: PetscCheck(pforest->forest, PetscObjectComm(odm), PETSC_ERR_ARG_WRONG, "DM has not been setup with a valid forest");
1392: name = vtk->filename;
1393: PetscCall(PetscStrlen(name, &len));
1394: PetscCall(PetscStrcasecmp(name + len - 4, ".vtu", &hasExt));
1395: if (hasExt) {
1396: PetscCall(PetscStrallocpy(name, &filenameStrip));
1397: filenameStrip[len - 4] = '\0';
1398: name = filenameStrip;
1399: }
1400: if (!pforest->topo->geom) PetscCallP4estReturn(geom, p4est_geometry_new_connectivity, (pforest->topo->conn));
1401: {
1402: p4est_vtk_context_t *pvtk;
1403: int footerr;
1405: PetscCallP4estReturn(pvtk, p4est_vtk_context_new, (pforest->forest, name));
1406: PetscCallP4est(p4est_vtk_context_set_geom, (pvtk, geom));
1407: PetscCallP4est(p4est_vtk_context_set_scale, (pvtk, (double)vtkScale));
1408: PetscCallP4estReturn(pvtk, p4est_vtk_write_header, (pvtk));
1409: PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_header() failed");
1410: PetscCallP4estReturn(pvtk, p4est_vtk_write_cell_dataf,
1411: (pvtk, 1, /* write tree */
1412: 1, /* write level */
1413: 1, /* write rank */
1414: 0, /* do not wrap rank */
1415: 0, /* no scalar fields */
1416: 0, /* no vector fields */
1417: pvtk));
1418: PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_cell_dataf() failed");
1419: PetscCallP4estReturn(footerr, p4est_vtk_write_footer, (pvtk));
1420: PetscCheck(!footerr, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_footer() failed");
1421: }
1422: if (!pforest->topo->geom) PetscCallP4est(p4est_geometry_destroy, (geom));
1423: PetscCall(PetscFree(filenameStrip));
1424: break;
1425: default:
1426: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
1427: }
1428: PetscFunctionReturn(PETSC_SUCCESS);
1429: }
1431: #define DMView_HDF5_pforest _append_pforest(DMView_HDF5)
1432: static PetscErrorCode DMView_HDF5_pforest(DM dm, PetscViewer viewer)
1433: {
1434: DM plex;
1436: PetscFunctionBegin;
1437: PetscCall(DMSetUp(dm));
1438: PetscCall(DMPforestGetPlex(dm, &plex));
1439: PetscCall(DMView(plex, viewer));
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: #define DMView_GLVis_pforest _append_pforest(DMView_GLVis)
1444: static PetscErrorCode DMView_GLVis_pforest(DM dm, PetscViewer viewer)
1445: {
1446: DM plex;
1448: PetscFunctionBegin;
1449: PetscCall(DMSetUp(dm));
1450: PetscCall(DMPforestGetPlex(dm, &plex));
1451: PetscCall(DMView(plex, viewer));
1452: PetscFunctionReturn(PETSC_SUCCESS);
1453: }
1455: #define DMView_pforest _append_pforest(DMView)
1456: static PetscErrorCode DMView_pforest(DM dm, PetscViewer viewer)
1457: {
1458: PetscBool isascii, isvtk, ishdf5, isglvis;
1460: PetscFunctionBegin;
1463: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1464: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1465: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1466: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
1467: if (isascii) {
1468: PetscCall(DMView_ASCII_pforest((PetscObject)dm, viewer));
1469: } else if (isvtk) {
1470: PetscCall(DMView_VTK_pforest((PetscObject)dm, viewer));
1471: } else if (ishdf5) {
1472: PetscCall(DMView_HDF5_pforest(dm, viewer));
1473: } else if (isglvis) {
1474: PetscCall(DMView_GLVis_pforest(dm, viewer));
1475: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Viewer not supported (not VTK, HDF5, or GLVis)");
1476: PetscFunctionReturn(PETSC_SUCCESS);
1477: }
1479: static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *conn, PetscInt **tree_face_to_uniq)
1480: {
1481: PetscInt *ttf, f, t, g, count;
1482: PetscInt numFacets;
1484: PetscFunctionBegin;
1485: numFacets = conn->num_trees * P4EST_FACES;
1486: PetscCall(PetscMalloc1(numFacets, &ttf));
1487: for (f = 0; f < numFacets; f++) ttf[f] = -1;
1488: for (g = 0, count = 0, t = 0; t < conn->num_trees; t++) {
1489: for (f = 0; f < P4EST_FACES; f++, g++) {
1490: if (ttf[g] == -1) {
1491: PetscInt ng;
1493: ttf[g] = count++;
1494: ng = conn->tree_to_tree[g] * P4EST_FACES + (conn->tree_to_face[g] % P4EST_FACES);
1495: ttf[ng] = ttf[g];
1496: }
1497: }
1498: }
1499: *tree_face_to_uniq = ttf;
1500: PetscFunctionReturn(PETSC_SUCCESS);
1501: }
1503: static PetscErrorCode DMPlexCreateConnectivity_pforest(DM dm, p4est_connectivity_t **connOut, PetscInt **tree_face_to_uniq)
1504: {
1505: p4est_topidx_t numTrees, numVerts, numCorns, numCtt;
1506: PetscSection ctt;
1507: #if defined(P4_TO_P8)
1508: p4est_topidx_t numEdges, numEtt;
1509: PetscSection ett;
1510: PetscInt eStart, eEnd, e, ettSize;
1511: PetscInt vertOff = 1 + P4EST_FACES + P8EST_EDGES;
1512: PetscInt edgeOff = 1 + P4EST_FACES;
1513: #else
1514: PetscInt vertOff = 1 + P4EST_FACES;
1515: #endif
1516: p4est_connectivity_t *conn;
1517: PetscInt cStart, cEnd, c, vStart, vEnd, v, fStart, fEnd, f;
1518: PetscInt *star = NULL, *closure = NULL, closureSize, starSize, cttSize;
1519: PetscInt *ttf;
1521: PetscFunctionBegin;
1522: /* 1: count objects, allocate */
1523: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1524: PetscCall(P4estTopidxCast(cEnd - cStart, &numTrees));
1525: numVerts = P4EST_CHILDREN * numTrees;
1526: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1527: PetscCall(P4estTopidxCast(vEnd - vStart, &numCorns));
1528: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ctt));
1529: PetscCall(PetscSectionSetChart(ctt, vStart, vEnd));
1530: for (v = vStart; v < vEnd; v++) {
1531: PetscInt s;
1533: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1534: for (s = 0; s < starSize; s++) {
1535: PetscInt p = star[2 * s];
1537: if (p >= cStart && p < cEnd) {
1538: /* we want to count every time cell p references v, so we see how many times it comes up in the closure. This
1539: * only protects against periodicity problems */
1540: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1541: 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);
1542: for (c = 0; c < P4EST_CHILDREN; c++) {
1543: PetscInt cellVert = closure[2 * (c + vertOff)];
1545: PetscCheck(cellVert >= vStart && cellVert < vEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: vertices");
1546: if (cellVert == v) PetscCall(PetscSectionAddDof(ctt, v, 1));
1547: }
1548: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1549: }
1550: }
1551: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1552: }
1553: PetscCall(PetscSectionSetUp(ctt));
1554: PetscCall(PetscSectionGetStorageSize(ctt, &cttSize));
1555: PetscCall(P4estTopidxCast(cttSize, &numCtt));
1556: #if defined(P4_TO_P8)
1557: PetscCall(DMPlexGetSimplexOrBoxCells(dm, P4EST_DIM - 1, &eStart, &eEnd));
1558: PetscCall(P4estTopidxCast(eEnd - eStart, &numEdges));
1559: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ett));
1560: PetscCall(PetscSectionSetChart(ett, eStart, eEnd));
1561: for (e = eStart; e < eEnd; e++) {
1562: PetscInt s;
1564: PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1565: for (s = 0; s < starSize; s++) {
1566: PetscInt p = star[2 * s];
1568: if (p >= cStart && p < cEnd) {
1569: /* we want to count every time cell p references e, so we see how many times it comes up in the closure. This
1570: * only protects against periodicity problems */
1571: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1572: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell with wrong closure size");
1573: for (c = 0; c < P8EST_EDGES; c++) {
1574: PetscInt cellEdge = closure[2 * (c + edgeOff)];
1576: PetscCheck(cellEdge >= eStart && cellEdge < eEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: edges");
1577: if (cellEdge == e) PetscCall(PetscSectionAddDof(ett, e, 1));
1578: }
1579: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1580: }
1581: }
1582: PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1583: }
1584: PetscCall(PetscSectionSetUp(ett));
1585: PetscCall(PetscSectionGetStorageSize(ett, &ettSize));
1586: PetscCall(P4estTopidxCast(ettSize, &numEtt));
1588: /* This routine allocates space for the arrays, which we fill below */
1589: PetscCallP4estReturn(conn, p8est_connectivity_new, (numVerts, numTrees, numEdges, numEtt, numCorns, numCtt));
1590: #else
1591: PetscCallP4estReturn(conn, p4est_connectivity_new, (numVerts, numTrees, numCorns, numCtt));
1592: #endif
1594: /* 2: visit every face, determine neighboring cells(trees) */
1595: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 1, &fStart, &fEnd));
1596: PetscCall(PetscMalloc1((cEnd - cStart) * P4EST_FACES, &ttf));
1597: for (f = fStart; f < fEnd; f++) {
1598: PetscInt numSupp, s;
1599: PetscInt myFace[2] = {-1, -1};
1600: PetscInt myOrnt[2] = {PETSC_MIN_INT, PETSC_MIN_INT};
1601: const PetscInt *supp;
1603: PetscCall(DMPlexGetSupportSize(dm, f, &numSupp));
1604: 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);
1605: PetscCall(DMPlexGetSupport(dm, f, &supp));
1607: for (s = 0; s < numSupp; s++) {
1608: PetscInt p = supp[s];
1610: if (p >= cEnd) {
1611: numSupp--;
1612: if (s) supp = &supp[1 - s];
1613: break;
1614: }
1615: }
1616: for (s = 0; s < numSupp; s++) {
1617: PetscInt p = supp[s], i;
1618: PetscInt numCone;
1619: DMPolytopeType ct;
1620: const PetscInt *cone;
1621: const PetscInt *ornt;
1622: PetscInt orient = PETSC_MIN_INT;
1624: PetscCall(DMPlexGetConeSize(dm, p, &numCone));
1625: PetscCheck(numCone == P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " has %" PetscInt_FMT " facets, expect %d", p, numCone, P4EST_FACES);
1626: PetscCall(DMPlexGetCone(dm, p, &cone));
1627: PetscCall(DMPlexGetCellType(dm, cone[0], &ct));
1628: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
1629: for (i = 0; i < P4EST_FACES; i++) {
1630: if (cone[i] == f) {
1631: orient = DMPolytopeConvertNewOrientation_Internal(ct, ornt[i]);
1632: break;
1633: }
1634: }
1635: PetscCheck(i < P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " faced %" PetscInt_FMT " mismatch", p, f);
1636: if (p < cStart || p >= cEnd) {
1637: DMPolytopeType ct;
1638: PetscCall(DMPlexGetCellType(dm, p, &ct));
1639: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " (%s) should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], cStart, cEnd);
1640: }
1641: ttf[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = f - fStart;
1642: if (numSupp == 1) {
1643: /* boundary faces indicated by self reference */
1644: conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = p - cStart;
1645: conn->tree_to_face[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = (int8_t)PetscFaceToP4estFace[i];
1646: } else {
1647: const PetscInt N = P4EST_CHILDREN / 2;
1649: conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = supp[1 - s] - cStart;
1650: myFace[s] = PetscFaceToP4estFace[i];
1651: /* get the orientation of cell p in p4est-type closure to facet f, by composing the p4est-closure to
1652: * petsc-closure permutation and the petsc-closure to facet orientation */
1653: myOrnt[s] = DihedralCompose(N, orient, DMPolytopeConvertNewOrientation_Internal(ct, P4estFaceToPetscOrnt[myFace[s]]));
1654: }
1655: }
1656: if (numSupp == 2) {
1657: for (s = 0; s < numSupp; s++) {
1658: PetscInt p = supp[s];
1659: PetscInt orntAtoB;
1660: PetscInt p4estOrient;
1661: const PetscInt N = P4EST_CHILDREN / 2;
1663: /* composing the forward permutation with the other cell's inverse permutation gives the self-to-neighbor
1664: * permutation of this cell-facet's cone */
1665: orntAtoB = DihedralCompose(N, DihedralInvert(N, myOrnt[1 - s]), myOrnt[s]);
1667: /* convert cone-description permutation (i.e., edges around facet) to cap-description permutation (i.e.,
1668: * vertices around facet) */
1669: #if !defined(P4_TO_P8)
1670: p4estOrient = orntAtoB < 0 ? -(orntAtoB + 1) : orntAtoB;
1671: #else
1672: {
1673: PetscInt firstVert = orntAtoB < 0 ? ((-orntAtoB) % N) : orntAtoB;
1674: PetscInt p4estFirstVert = firstVert < 2 ? firstVert : (firstVert ^ 1);
1676: /* swap bits */
1677: p4estOrient = ((myFace[s] <= myFace[1 - s]) || (orntAtoB < 0)) ? p4estFirstVert : ((p4estFirstVert >> 1) | ((p4estFirstVert & 1) << 1));
1678: }
1679: #endif
1680: /* encode neighbor face and orientation in tree_to_face per p4est_connectivity standard (see
1681: * p4est_connectivity.h, p8est_connectivity.h) */
1682: conn->tree_to_face[P4EST_FACES * (p - cStart) + myFace[s]] = (int8_t)myFace[1 - s] + p4estOrient * P4EST_FACES;
1683: }
1684: }
1685: }
1687: #if defined(P4_TO_P8)
1688: /* 3: visit every edge */
1689: conn->ett_offset[0] = 0;
1690: for (e = eStart; e < eEnd; e++) {
1691: PetscInt off, s;
1693: PetscCall(PetscSectionGetOffset(ett, e, &off));
1694: conn->ett_offset[e - eStart] = (p4est_topidx_t)off;
1695: PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1696: for (s = 0; s < starSize; s++) {
1697: PetscInt p = star[2 * s];
1699: if (p >= cStart && p < cEnd) {
1700: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1701: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure");
1702: for (c = 0; c < P8EST_EDGES; c++) {
1703: PetscInt cellEdge = closure[2 * (c + edgeOff)];
1704: PetscInt cellOrnt = closure[2 * (c + edgeOff) + 1];
1705: DMPolytopeType ct;
1707: PetscCall(DMPlexGetCellType(dm, cellEdge, &ct));
1708: cellOrnt = DMPolytopeConvertNewOrientation_Internal(ct, cellOrnt);
1709: if (cellEdge == e) {
1710: PetscInt p4estEdge = PetscEdgeToP4estEdge[c];
1711: PetscInt totalOrient;
1713: /* compose p4est-closure to petsc-closure permutation and petsc-closure to edge orientation */
1714: totalOrient = DihedralCompose(2, cellOrnt, DMPolytopeConvertNewOrientation_Internal(DM_POLYTOPE_SEGMENT, P4estEdgeToPetscOrnt[p4estEdge]));
1715: /* p4est orientations are positive: -2 => 1, -1 => 0 */
1716: totalOrient = (totalOrient < 0) ? -(totalOrient + 1) : totalOrient;
1717: conn->edge_to_tree[off] = (p4est_locidx_t)(p - cStart);
1718: /* encode cell-edge and orientation in edge_to_edge per p8est_connectivity standard (see
1719: * p8est_connectivity.h) */
1720: conn->edge_to_edge[off++] = (int8_t)p4estEdge + P8EST_EDGES * totalOrient;
1721: conn->tree_to_edge[P8EST_EDGES * (p - cStart) + p4estEdge] = e - eStart;
1722: }
1723: }
1724: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1725: }
1726: }
1727: PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star));
1728: }
1729: PetscCall(PetscSectionDestroy(&ett));
1730: #endif
1732: /* 4: visit every vertex */
1733: conn->ctt_offset[0] = 0;
1734: for (v = vStart; v < vEnd; v++) {
1735: PetscInt off, s;
1737: PetscCall(PetscSectionGetOffset(ctt, v, &off));
1738: conn->ctt_offset[v - vStart] = (p4est_topidx_t)off;
1739: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1740: for (s = 0; s < starSize; s++) {
1741: PetscInt p = star[2 * s];
1743: if (p >= cStart && p < cEnd) {
1744: PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1745: PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure");
1746: for (c = 0; c < P4EST_CHILDREN; c++) {
1747: PetscInt cellVert = closure[2 * (c + vertOff)];
1749: if (cellVert == v) {
1750: PetscInt p4estVert = PetscVertToP4estVert[c];
1752: conn->corner_to_tree[off] = (p4est_locidx_t)(p - cStart);
1753: conn->corner_to_corner[off++] = (int8_t)p4estVert;
1754: conn->tree_to_corner[P4EST_CHILDREN * (p - cStart) + p4estVert] = v - vStart;
1755: }
1756: }
1757: PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
1758: }
1759: }
1760: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
1761: }
1762: PetscCall(PetscSectionDestroy(&ctt));
1764: /* 5: Compute the coordinates */
1765: {
1766: PetscInt coordDim;
1768: PetscCall(DMGetCoordinateDim(dm, &coordDim));
1769: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1770: for (c = cStart; c < cEnd; c++) {
1771: PetscInt dof;
1772: PetscBool isDG;
1773: PetscScalar *cellCoords = NULL;
1774: const PetscScalar *array;
1776: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords));
1777: 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);
1778: for (v = 0; v < P4EST_CHILDREN; v++) {
1779: PetscInt i, lim = PetscMin(3, coordDim);
1780: PetscInt p4estVert = PetscVertToP4estVert[v];
1782: conn->tree_to_vertex[P4EST_CHILDREN * (c - cStart) + v] = P4EST_CHILDREN * (c - cStart) + v;
1783: /* p4est vertices are always embedded in R^3 */
1784: for (i = 0; i < 3; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = 0.;
1785: for (i = 0; i < lim; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = PetscRealPart(cellCoords[v * coordDim + i]);
1786: }
1787: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords));
1788: }
1789: }
1791: #if defined(P4EST_ENABLE_DEBUG)
1792: PetscCheck(p4est_connectivity_is_valid(conn), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Plex to p4est conversion failed");
1793: #endif
1795: *connOut = conn;
1797: *tree_face_to_uniq = ttf;
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: }
1801: static PetscErrorCode locidx_to_PetscInt(sc_array_t *array)
1802: {
1803: sc_array_t *newarray;
1804: size_t zz, count = array->elem_count;
1806: PetscFunctionBegin;
1807: PetscCheck(array->elem_size == sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size");
1809: if (sizeof(p4est_locidx_t) == sizeof(PetscInt)) PetscFunctionReturn(PETSC_SUCCESS);
1811: newarray = sc_array_new_size(sizeof(PetscInt), array->elem_count);
1812: for (zz = 0; zz < count; zz++) {
1813: p4est_locidx_t il = *((p4est_locidx_t *)sc_array_index(array, zz));
1814: PetscInt *ip = (PetscInt *)sc_array_index(newarray, zz);
1816: *ip = (PetscInt)il;
1817: }
1819: sc_array_reset(array);
1820: sc_array_init_size(array, sizeof(PetscInt), count);
1821: sc_array_copy(array, newarray);
1822: sc_array_destroy(newarray);
1823: PetscFunctionReturn(PETSC_SUCCESS);
1824: }
1826: static PetscErrorCode coords_double_to_PetscScalar(sc_array_t *array, PetscInt dim)
1827: {
1828: sc_array_t *newarray;
1829: size_t zz, count = array->elem_count;
1831: PetscFunctionBegin;
1832: PetscCheck(array->elem_size == 3 * sizeof(double), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong coordinate size");
1833: #if !defined(PETSC_USE_COMPLEX)
1834: if (sizeof(double) == sizeof(PetscScalar) && dim == 3) PetscFunctionReturn(PETSC_SUCCESS);
1835: #endif
1837: newarray = sc_array_new_size(dim * sizeof(PetscScalar), array->elem_count);
1838: for (zz = 0; zz < count; zz++) {
1839: int i;
1840: double *id = (double *)sc_array_index(array, zz);
1841: PetscScalar *ip = (PetscScalar *)sc_array_index(newarray, zz);
1843: for (i = 0; i < dim; i++) ip[i] = 0.;
1844: for (i = 0; i < PetscMin(dim, 3); i++) ip[i] = (PetscScalar)id[i];
1845: }
1847: sc_array_reset(array);
1848: sc_array_init_size(array, dim * sizeof(PetscScalar), count);
1849: sc_array_copy(array, newarray);
1850: sc_array_destroy(newarray);
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: static PetscErrorCode locidx_pair_to_PetscSFNode(sc_array_t *array)
1855: {
1856: sc_array_t *newarray;
1857: size_t zz, count = array->elem_count;
1859: PetscFunctionBegin;
1860: PetscCheck(array->elem_size == 2 * sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size");
1862: newarray = sc_array_new_size(sizeof(PetscSFNode), array->elem_count);
1863: for (zz = 0; zz < count; zz++) {
1864: p4est_locidx_t *il = (p4est_locidx_t *)sc_array_index(array, zz);
1865: PetscSFNode *ip = (PetscSFNode *)sc_array_index(newarray, zz);
1867: ip->rank = (PetscInt)il[0];
1868: ip->index = (PetscInt)il[1];
1869: }
1871: sc_array_reset(array);
1872: sc_array_init_size(array, sizeof(PetscSFNode), count);
1873: sc_array_copy(array, newarray);
1874: sc_array_destroy(newarray);
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: static PetscErrorCode P4estToPlex_Local(p4est_t *p4est, DM *plex)
1879: {
1880: PetscFunctionBegin;
1881: {
1882: sc_array_t *points_per_dim = sc_array_new(sizeof(p4est_locidx_t));
1883: sc_array_t *cone_sizes = sc_array_new(sizeof(p4est_locidx_t));
1884: sc_array_t *cones = sc_array_new(sizeof(p4est_locidx_t));
1885: sc_array_t *cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
1886: sc_array_t *coords = sc_array_new(3 * sizeof(double));
1887: sc_array_t *children = sc_array_new(sizeof(p4est_locidx_t));
1888: sc_array_t *parents = sc_array_new(sizeof(p4est_locidx_t));
1889: sc_array_t *childids = sc_array_new(sizeof(p4est_locidx_t));
1890: sc_array_t *leaves = sc_array_new(sizeof(p4est_locidx_t));
1891: sc_array_t *remotes = sc_array_new(2 * sizeof(p4est_locidx_t));
1892: p4est_locidx_t first_local_quad;
1894: 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));
1896: PetscCall(locidx_to_PetscInt(points_per_dim));
1897: PetscCall(locidx_to_PetscInt(cone_sizes));
1898: PetscCall(locidx_to_PetscInt(cones));
1899: PetscCall(locidx_to_PetscInt(cone_orientations));
1900: PetscCall(coords_double_to_PetscScalar(coords, P4EST_DIM));
1902: PetscCall(DMPlexCreate(PETSC_COMM_SELF, plex));
1903: PetscCall(DMSetDimension(*plex, P4EST_DIM));
1904: PetscCall(DMPlexCreateFromDAG(*plex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array));
1905: PetscCall(DMPlexConvertOldOrientations_Internal(*plex));
1906: sc_array_destroy(points_per_dim);
1907: sc_array_destroy(cone_sizes);
1908: sc_array_destroy(cones);
1909: sc_array_destroy(cone_orientations);
1910: sc_array_destroy(coords);
1911: sc_array_destroy(children);
1912: sc_array_destroy(parents);
1913: sc_array_destroy(childids);
1914: sc_array_destroy(leaves);
1915: sc_array_destroy(remotes);
1916: }
1917: PetscFunctionReturn(PETSC_SUCCESS);
1918: }
1920: #define DMReferenceTreeGetChildSymmetry_pforest _append_pforest(DMReferenceTreeGetChildSymmetry)
1921: static PetscErrorCode DMReferenceTreeGetChildSymmetry_pforest(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
1922: {
1923: PetscInt coneSize, dStart, dEnd, vStart, vEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
1925: PetscFunctionBegin;
1926: if (parentOrientA == parentOrientB) {
1927: if (childOrientB) *childOrientB = childOrientA;
1928: if (childB) *childB = childA;
1929: PetscFunctionReturn(PETSC_SUCCESS);
1930: }
1931: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1932: if (childA >= vStart && childA < vEnd) { /* vertices (always in the middle) are invariant under rotation */
1933: if (childOrientB) *childOrientB = 0;
1934: if (childB) *childB = childA;
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1937: for (dim = 0; dim < 3; dim++) {
1938: PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd));
1939: if (parent >= dStart && parent <= dEnd) break;
1940: }
1941: PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim);
1942: PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children");
1943: if (childA < dStart || childA >= dEnd) { /* a 1-cell in a 2-cell */
1944: /* this is a lower-dimensional child: bootstrap */
1945: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
1946: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
1948: PetscCall(DMPlexGetSupportSize(dm, childA, &size));
1949: PetscCall(DMPlexGetSupport(dm, childA, &supp));
1951: /* find a point sA in supp(childA) that has the same parent */
1952: for (i = 0; i < size; i++) {
1953: PetscInt sParent;
1955: sA = supp[i];
1956: if (sA == parent) continue;
1957: PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL));
1958: if (sParent == parent) break;
1959: }
1960: PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children");
1961: /* find out which point sB is in an equivalent position to sA under
1962: * parentOrientB */
1963: PetscCall(DMReferenceTreeGetChildSymmetry_pforest(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB));
1964: PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize));
1965: PetscCall(DMPlexGetCone(dm, sA, &coneA));
1966: PetscCall(DMPlexGetCone(dm, sB, &coneB));
1967: PetscCall(DMPlexGetConeOrientation(dm, sA, &oA));
1968: PetscCall(DMPlexGetConeOrientation(dm, sB, &oB));
1969: /* step through the cone of sA in natural order */
1970: for (i = 0; i < sConeSize; i++) {
1971: if (coneA[i] == childA) {
1972: /* if childA is at position i in coneA,
1973: * then we want the point that is at sOrientB*i in coneB */
1974: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
1975: if (childB) *childB = coneB[j];
1976: if (childOrientB) {
1977: DMPolytopeType ct;
1978: PetscInt oBtrue;
1980: PetscCall(DMPlexGetConeSize(dm, childA, &coneSize));
1981: /* compose sOrientB and oB[j] */
1982: PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge");
1983: ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1984: /* we may have to flip an edge */
1985: oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientationInv(ct, -1, oB[j]);
1986: oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
1987: ABswap = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
1988: *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
1989: }
1990: break;
1991: }
1992: }
1993: PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch");
1994: PetscFunctionReturn(PETSC_SUCCESS);
1995: }
1996: /* get the cone size and symmetry swap */
1997: PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
1998: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
1999: if (dim == 2) {
2000: /* orientations refer to cones: we want them to refer to vertices:
2001: * if it's a rotation, they are the same, but if the order is reversed, a
2002: * permutation that puts side i first does *not* put vertex i first */
2003: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
2004: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
2005: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
2006: } else {
2007: oAvert = parentOrientA;
2008: oBvert = parentOrientB;
2009: ABswapVert = ABswap;
2010: }
2011: if (childB) {
2012: /* assume that each child corresponds to a vertex, in the same order */
2013: PetscInt p, posA = -1, numChildren, i;
2014: const PetscInt *children;
2016: /* count which position the child is in */
2017: PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children));
2018: for (i = 0; i < numChildren; i++) {
2019: p = children[i];
2020: if (p == childA) {
2021: if (dim == 1) {
2022: posA = i;
2023: } else { /* 2D Morton to rotation */
2024: posA = (i & 2) ? (i ^ 1) : i;
2025: }
2026: break;
2027: }
2028: }
2029: if (posA >= coneSize) {
2030: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find childA in children of parent");
2031: } else {
2032: /* figure out position B by applying ABswapVert */
2033: PetscInt posB, childIdB;
2035: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
2036: if (dim == 1) {
2037: childIdB = posB;
2038: } else { /* 2D rotation to Morton */
2039: childIdB = (posB & 2) ? (posB ^ 1) : posB;
2040: }
2041: if (childB) *childB = children[childIdB];
2042: }
2043: }
2044: if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
2045: PetscFunctionReturn(PETSC_SUCCESS);
2046: }
2048: #define DMCreateReferenceTree_pforest _append_pforest(DMCreateReferenceTree)
2049: static PetscErrorCode DMCreateReferenceTree_pforest(MPI_Comm comm, DM *dm)
2050: {
2051: p4est_connectivity_t *refcube;
2052: p4est_t *root, *refined;
2053: DM dmRoot, dmRefined;
2054: DM_Plex *mesh;
2055: PetscMPIInt rank;
2056: #if defined(PETSC_HAVE_MPIUNI)
2057: sc_MPI_Comm comm_self = sc_MPI_COMM_SELF;
2058: #else
2059: MPI_Comm comm_self = PETSC_COMM_SELF;
2060: #endif
2062: PetscFunctionBegin;
2063: PetscCallP4estReturn(refcube, p4est_connectivity_new_byname, ("unit"));
2064: { /* [-1,1]^d geometry */
2065: PetscInt i, j;
2067: for (i = 0; i < P4EST_CHILDREN; i++) {
2068: for (j = 0; j < 3; j++) {
2069: refcube->vertices[3 * i + j] *= 2.;
2070: refcube->vertices[3 * i + j] -= 1.;
2071: }
2072: }
2073: }
2074: PetscCallP4estReturn(root, p4est_new, (comm_self, refcube, 0, NULL, NULL));
2075: PetscCallP4estReturn(refined, p4est_new_ext, (comm_self, refcube, 0, 1, 1, 0, NULL, NULL));
2076: PetscCall(P4estToPlex_Local(root, &dmRoot));
2077: PetscCall(P4estToPlex_Local(refined, &dmRefined));
2078: {
2079: #if !defined(P4_TO_P8)
2080: PetscInt nPoints = 25;
2081: 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};
2082: 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};
2083: #else
2084: PetscInt nPoints = 125;
2085: 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,
2086: 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,
2087: 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};
2088: 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,
2089: 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};
2091: #endif
2092: IS permIS;
2093: DM dmPerm;
2095: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nPoints, perm, PETSC_USE_POINTER, &permIS));
2096: PetscCall(DMPlexPermute(dmRefined, permIS, &dmPerm));
2097: if (dmPerm) {
2098: PetscCall(DMDestroy(&dmRefined));
2099: dmRefined = dmPerm;
2100: }
2101: PetscCall(ISDestroy(&permIS));
2102: {
2103: PetscInt p;
2104: PetscCall(DMCreateLabel(dmRoot, "identity"));
2105: PetscCall(DMCreateLabel(dmRefined, "identity"));
2106: for (p = 0; p < P4EST_INSUL; p++) PetscCall(DMSetLabelValue(dmRoot, "identity", p, p));
2107: for (p = 0; p < nPoints; p++) PetscCall(DMSetLabelValue(dmRefined, "identity", p, ident[p]));
2108: }
2109: }
2110: PetscCall(DMPlexCreateReferenceTree_Union(dmRoot, dmRefined, "identity", dm));
2111: mesh = (DM_Plex *)(*dm)->data;
2112: mesh->getchildsymmetry = DMReferenceTreeGetChildSymmetry_pforest;
2113: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2114: if (rank == 0) {
2115: PetscCall(DMViewFromOptions(dmRoot, NULL, "-dm_p4est_ref_root_view"));
2116: PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_refined_view"));
2117: PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_tree_view"));
2118: }
2119: PetscCall(DMDestroy(&dmRefined));
2120: PetscCall(DMDestroy(&dmRoot));
2121: PetscCallP4est(p4est_destroy, (refined));
2122: PetscCallP4est(p4est_destroy, (root));
2123: PetscCallP4est(p4est_connectivity_destroy, (refcube));
2124: PetscFunctionReturn(PETSC_SUCCESS);
2125: }
2127: static PetscErrorCode DMShareDiscretization(DM dmA, DM dmB)
2128: {
2129: void *ctx;
2130: PetscInt num;
2131: PetscReal val;
2133: PetscFunctionBegin;
2134: PetscCall(DMGetApplicationContext(dmA, &ctx));
2135: PetscCall(DMSetApplicationContext(dmB, ctx));
2136: PetscCall(DMCopyDisc(dmA, dmB));
2137: PetscCall(DMGetOutputSequenceNumber(dmA, &num, &val));
2138: PetscCall(DMSetOutputSequenceNumber(dmB, num, val));
2139: if (dmB->localSection != dmA->localSection || dmB->globalSection != dmA->globalSection) {
2140: PetscCall(DMClearLocalVectors(dmB));
2141: PetscCall(PetscObjectReference((PetscObject)dmA->localSection));
2142: PetscCall(PetscSectionDestroy(&dmB->localSection));
2143: dmB->localSection = dmA->localSection;
2144: PetscCall(DMClearGlobalVectors(dmB));
2145: PetscCall(PetscObjectReference((PetscObject)dmA->globalSection));
2146: PetscCall(PetscSectionDestroy(&dmB->globalSection));
2147: dmB->globalSection = dmA->globalSection;
2148: PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.section));
2149: PetscCall(PetscSectionDestroy(&dmB->defaultConstraint.section));
2150: dmB->defaultConstraint.section = dmA->defaultConstraint.section;
2151: PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.mat));
2152: PetscCall(MatDestroy(&dmB->defaultConstraint.mat));
2153: dmB->defaultConstraint.mat = dmA->defaultConstraint.mat;
2154: if (dmA->map) PetscCall(PetscLayoutReference(dmA->map, &dmB->map));
2155: }
2156: if (dmB->sectionSF != dmA->sectionSF) {
2157: PetscCall(PetscObjectReference((PetscObject)dmA->sectionSF));
2158: PetscCall(PetscSFDestroy(&dmB->sectionSF));
2159: dmB->sectionSF = dmA->sectionSF;
2160: }
2161: PetscFunctionReturn(PETSC_SUCCESS);
2162: }
2164: /* Get an SF that broadcasts a coarse-cell covering of the local fine cells */
2165: static PetscErrorCode DMPforestGetCellCoveringSF(MPI_Comm comm, p4est_t *p4estC, p4est_t *p4estF, PetscInt cStart, PetscInt cEnd, PetscSF *coveringSF)
2166: {
2167: PetscInt startF, endF, startC, endC, p, nLeaves;
2168: PetscSFNode *leaves;
2169: PetscSF sf;
2170: PetscInt *recv, *send;
2171: PetscMPIInt tag;
2172: MPI_Request *recvReqs, *sendReqs;
2173: PetscSection section;
2175: PetscFunctionBegin;
2176: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estF, p4estC, &startC, &endC));
2177: PetscCall(PetscMalloc2(2 * (endC - startC), &recv, endC - startC, &recvReqs));
2178: PetscCall(PetscCommGetNewTag(comm, &tag));
2179: for (p = startC; p < endC; p++) {
2180: recvReqs[p - startC] = MPI_REQUEST_NULL; /* just in case we don't initiate a receive */
2181: if (p4estC->global_first_quadrant[p] == p4estC->global_first_quadrant[p + 1]) { /* empty coarse partition */
2182: recv[2 * (p - startC)] = 0;
2183: recv[2 * (p - startC) + 1] = 0;
2184: continue;
2185: }
2187: PetscCallMPI(MPI_Irecv(&recv[2 * (p - startC)], 2, MPIU_INT, p, tag, comm, &recvReqs[p - startC]));
2188: }
2189: PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estC, p4estF, &startF, &endF));
2190: PetscCall(PetscMalloc2(2 * (endF - startF), &send, endF - startF, &sendReqs));
2191: /* count the quadrants rank will send to each of [startF,endF) */
2192: for (p = startF; p < endF; p++) {
2193: p4est_quadrant_t *myFineStart = &p4estF->global_first_position[p];
2194: p4est_quadrant_t *myFineEnd = &p4estF->global_first_position[p + 1];
2195: PetscInt tStart = (PetscInt)myFineStart->p.which_tree;
2196: PetscInt tEnd = (PetscInt)myFineEnd->p.which_tree;
2197: PetscInt firstCell = -1, lastCell = -1;
2198: p4est_tree_t *treeStart = &(((p4est_tree_t *)p4estC->trees->array)[tStart]);
2199: p4est_tree_t *treeEnd = (size_t)tEnd < p4estC->trees->elem_count ? &(((p4est_tree_t *)p4estC->trees->array)[tEnd]) : NULL;
2200: ssize_t overlapIndex;
2202: sendReqs[p - startF] = MPI_REQUEST_NULL; /* just in case we don't initiate a send */
2203: if (p4estF->global_first_quadrant[p] == p4estF->global_first_quadrant[p + 1]) continue;
2205: /* locate myFineStart in (or before) a cell */
2206: if (treeStart->quadrants.elem_count) {
2207: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeStart->quadrants, myFineStart, p4est_quadrant_disjoint));
2208: if (overlapIndex < 0) {
2209: firstCell = 0;
2210: } else {
2211: firstCell = treeStart->quadrants_offset + overlapIndex;
2212: }
2213: } else {
2214: firstCell = 0;
2215: }
2216: if (treeEnd && treeEnd->quadrants.elem_count) {
2217: PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeEnd->quadrants, myFineEnd, p4est_quadrant_disjoint));
2218: if (overlapIndex < 0) { /* all of this local section is overlapped */
2219: lastCell = p4estC->local_num_quadrants;
2220: } else {
2221: p4est_quadrant_t *container = &(((p4est_quadrant_t *)treeEnd->quadrants.array)[overlapIndex]);
2222: p4est_quadrant_t first_desc;
2223: int equal;
2225: PetscCallP4est(p4est_quadrant_first_descendant, (container, &first_desc, P4EST_QMAXLEVEL));
2226: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (myFineEnd, &first_desc));
2227: if (equal) {
2228: lastCell = treeEnd->quadrants_offset + overlapIndex;
2229: } else {
2230: lastCell = treeEnd->quadrants_offset + overlapIndex + 1;
2231: }
2232: }
2233: } else {
2234: lastCell = p4estC->local_num_quadrants;
2235: }
2236: send[2 * (p - startF)] = firstCell;
2237: send[2 * (p - startF) + 1] = lastCell - firstCell;
2238: PetscCallMPI(MPI_Isend(&send[2 * (p - startF)], 2, MPIU_INT, p, tag, comm, &sendReqs[p - startF]));
2239: }
2240: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endC - startC), recvReqs, MPI_STATUSES_IGNORE));
2241: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
2242: PetscCall(PetscSectionSetChart(section, startC, endC));
2243: for (p = startC; p < endC; p++) {
2244: PetscInt numCells = recv[2 * (p - startC) + 1];
2245: PetscCall(PetscSectionSetDof(section, p, numCells));
2246: }
2247: PetscCall(PetscSectionSetUp(section));
2248: PetscCall(PetscSectionGetStorageSize(section, &nLeaves));
2249: PetscCall(PetscMalloc1(nLeaves, &leaves));
2250: for (p = startC; p < endC; p++) {
2251: PetscInt firstCell = recv[2 * (p - startC)];
2252: PetscInt numCells = recv[2 * (p - startC) + 1];
2253: PetscInt off, i;
2255: PetscCall(PetscSectionGetOffset(section, p, &off));
2256: for (i = 0; i < numCells; i++) {
2257: leaves[off + i].rank = p;
2258: leaves[off + i].index = firstCell + i;
2259: }
2260: }
2261: PetscCall(PetscSFCreate(comm, &sf));
2262: PetscCall(PetscSFSetGraph(sf, cEnd - cStart, nLeaves, NULL, PETSC_OWN_POINTER, leaves, PETSC_OWN_POINTER));
2263: PetscCall(PetscSectionDestroy(§ion));
2264: PetscCallMPI(MPI_Waitall((PetscMPIInt)(endF - startF), sendReqs, MPI_STATUSES_IGNORE));
2265: PetscCall(PetscFree2(send, sendReqs));
2266: PetscCall(PetscFree2(recv, recvReqs));
2267: *coveringSF = sf;
2268: PetscFunctionReturn(PETSC_SUCCESS);
2269: }
2271: /* closure points for locally-owned cells */
2272: static PetscErrorCode DMPforestGetCellSFNodes(DM dm, PetscInt numClosureIndices, PetscInt *numClosurePoints, PetscSFNode **closurePoints, PetscBool redirect)
2273: {
2274: PetscInt cStart, cEnd;
2275: PetscInt count, c;
2276: PetscMPIInt rank;
2277: PetscInt closureSize = -1;
2278: PetscInt *closure = NULL;
2279: PetscSF pointSF;
2280: PetscInt nleaves, nroots;
2281: const PetscInt *ilocal;
2282: const PetscSFNode *iremote;
2283: DM plex;
2284: DM_Forest *forest;
2285: DM_Forest_pforest *pforest;
2287: PetscFunctionBegin;
2288: forest = (DM_Forest *)dm->data;
2289: pforest = (DM_Forest_pforest *)forest->data;
2290: cStart = pforest->cLocalStart;
2291: cEnd = pforest->cLocalEnd;
2292: PetscCall(DMPforestGetPlex(dm, &plex));
2293: PetscCall(DMGetPointSF(dm, &pointSF));
2294: PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
2295: nleaves = PetscMax(0, nleaves);
2296: nroots = PetscMax(0, nroots);
2297: *numClosurePoints = numClosureIndices * (cEnd - cStart);
2298: PetscCall(PetscMalloc1(*numClosurePoints, closurePoints));
2299: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2300: for (c = cStart, count = 0; c < cEnd; c++) {
2301: PetscInt i;
2302: PetscCall(DMPlexGetTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2304: for (i = 0; i < numClosureIndices; i++, count++) {
2305: PetscInt p = closure[2 * i];
2306: PetscInt loc = -1;
2308: PetscCall(PetscFindInt(p, nleaves, ilocal, &loc));
2309: if (redirect && loc >= 0) {
2310: (*closurePoints)[count].rank = iremote[loc].rank;
2311: (*closurePoints)[count].index = iremote[loc].index;
2312: } else {
2313: (*closurePoints)[count].rank = rank;
2314: (*closurePoints)[count].index = p;
2315: }
2316: }
2317: PetscCall(DMPlexRestoreTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure));
2318: }
2319: PetscFunctionReturn(PETSC_SUCCESS);
2320: }
2322: static void MPIAPI DMPforestMaxSFNode(void *a, void *b, PetscMPIInt *len, MPI_Datatype *type)
2323: {
2324: PetscMPIInt i;
2326: for (i = 0; i < *len; i++) {
2327: PetscSFNode *A = (PetscSFNode *)a;
2328: PetscSFNode *B = (PetscSFNode *)b;
2330: if (B->rank < 0) *B = *A;
2331: }
2332: }
2334: static PetscErrorCode DMPforestGetTransferSF_Point(DM coarse, DM fine, PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2335: {
2336: MPI_Comm comm;
2337: PetscMPIInt rank, size;
2338: DM_Forest_pforest *pforestC, *pforestF;
2339: p4est_t *p4estC, *p4estF;
2340: PetscInt numClosureIndices;
2341: PetscInt numClosurePointsC, numClosurePointsF;
2342: PetscSFNode *closurePointsC, *closurePointsF;
2343: p4est_quadrant_t *coverQuads = NULL;
2344: p4est_quadrant_t **treeQuads;
2345: PetscInt *treeQuadCounts;
2346: MPI_Datatype nodeType;
2347: MPI_Datatype nodeClosureType;
2348: MPI_Op sfNodeReduce;
2349: p4est_topidx_t fltF, lltF, t;
2350: DM plexC, plexF;
2351: PetscInt pStartF, pEndF, pStartC, pEndC;
2352: PetscBool saveInCoarse = PETSC_FALSE;
2353: PetscBool saveInFine = PETSC_FALSE;
2354: PetscBool formCids = (childIds != NULL) ? PETSC_TRUE : PETSC_FALSE;
2355: PetscInt *cids = NULL;
2357: PetscFunctionBegin;
2358: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2359: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2360: p4estC = pforestC->forest;
2361: p4estF = pforestF->forest;
2362: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2363: comm = PetscObjectComm((PetscObject)coarse);
2364: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2365: PetscCallMPI(MPI_Comm_size(comm, &size));
2366: PetscCall(DMPforestGetPlex(fine, &plexF));
2367: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2368: PetscCall(DMPforestGetPlex(coarse, &plexC));
2369: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2370: { /* check if the results have been cached */
2371: DM adaptCoarse, adaptFine;
2373: PetscCall(DMForestGetAdaptivityForest(coarse, &adaptCoarse));
2374: PetscCall(DMForestGetAdaptivityForest(fine, &adaptFine));
2375: if (adaptCoarse && adaptCoarse->data == fine->data) { /* coarse is adapted from fine */
2376: if (pforestC->pointSelfToAdaptSF) {
2377: PetscCall(PetscObjectReference((PetscObject)pforestC->pointSelfToAdaptSF));
2378: *sf = pforestC->pointSelfToAdaptSF;
2379: if (childIds) {
2380: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2381: PetscCall(PetscArraycpy(cids, pforestC->pointSelfToAdaptCids, pEndF - pStartF));
2382: *childIds = cids;
2383: }
2384: PetscFunctionReturn(PETSC_SUCCESS);
2385: } else {
2386: saveInCoarse = PETSC_TRUE;
2387: formCids = PETSC_TRUE;
2388: }
2389: } else if (adaptFine && adaptFine->data == coarse->data) { /* fine is adapted from coarse */
2390: if (pforestF->pointAdaptToSelfSF) {
2391: PetscCall(PetscObjectReference((PetscObject)pforestF->pointAdaptToSelfSF));
2392: *sf = pforestF->pointAdaptToSelfSF;
2393: if (childIds) {
2394: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2395: PetscCall(PetscArraycpy(cids, pforestF->pointAdaptToSelfCids, pEndF - pStartF));
2396: *childIds = cids;
2397: }
2398: PetscFunctionReturn(PETSC_SUCCESS);
2399: } else {
2400: saveInFine = PETSC_TRUE;
2401: formCids = PETSC_TRUE;
2402: }
2403: }
2404: }
2406: /* count the number of closure points that have dofs and create a list */
2407: numClosureIndices = P4EST_INSUL;
2408: /* create the datatype */
2409: PetscCallMPI(MPI_Type_contiguous(2, MPIU_INT, &nodeType));
2410: PetscCallMPI(MPI_Type_commit(&nodeType));
2411: PetscCallMPI(MPI_Op_create(DMPforestMaxSFNode, PETSC_FALSE, &sfNodeReduce));
2412: PetscCallMPI(MPI_Type_contiguous(numClosureIndices * 2, MPIU_INT, &nodeClosureType));
2413: PetscCallMPI(MPI_Type_commit(&nodeClosureType));
2414: /* everything has to go through cells: for each cell, create a list of the sfnodes in its closure */
2415: /* get lists of closure point SF nodes for every cell */
2416: PetscCall(DMPforestGetCellSFNodes(coarse, numClosureIndices, &numClosurePointsC, &closurePointsC, PETSC_TRUE));
2417: PetscCall(DMPforestGetCellSFNodes(fine, numClosureIndices, &numClosurePointsF, &closurePointsF, PETSC_FALSE));
2418: /* create pointers for tree lists */
2419: fltF = p4estF->first_local_tree;
2420: lltF = p4estF->last_local_tree;
2421: PetscCall(PetscCalloc2(lltF + 1 - fltF, &treeQuads, lltF + 1 - fltF, &treeQuadCounts));
2422: /* if the partitions don't match, ship the coarse to cover the fine */
2423: if (size > 1) {
2424: PetscInt p;
2426: for (p = 0; p < size; p++) {
2427: int equal;
2429: PetscCallP4estReturn(equal, p4est_quadrant_is_equal_piggy, (&p4estC->global_first_position[p], &p4estF->global_first_position[p]));
2430: if (!equal) break;
2431: }
2432: if (p < size) { /* non-matching distribution: send the coarse to cover the fine */
2433: PetscInt cStartC, cEndC;
2434: PetscSF coveringSF;
2435: PetscInt nleaves;
2436: PetscInt count;
2437: PetscSFNode *newClosurePointsC;
2438: p4est_quadrant_t *coverQuadsSend;
2439: p4est_topidx_t fltC = p4estC->first_local_tree;
2440: p4est_topidx_t lltC = p4estC->last_local_tree;
2441: p4est_topidx_t t;
2442: PetscMPIInt blockSizes[4] = {P4EST_DIM, 2, 1, 1};
2443: MPI_Aint blockOffsets[4] = {offsetof(p4est_quadrant_t, x), offsetof(p4est_quadrant_t, level), offsetof(p4est_quadrant_t, pad16), offsetof(p4est_quadrant_t, p)};
2444: MPI_Datatype blockTypes[4] = {MPI_INT32_T, MPI_INT8_T, MPI_INT16_T, MPI_INT32_T /* p.which_tree */};
2445: MPI_Datatype quadStruct, quadType;
2447: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, 0, &cStartC, &cEndC));
2448: PetscCall(DMPforestGetCellCoveringSF(comm, p4estC, p4estF, pforestC->cLocalStart, pforestC->cLocalEnd, &coveringSF));
2449: PetscCall(PetscSFGetGraph(coveringSF, NULL, &nleaves, NULL, NULL));
2450: PetscCall(PetscMalloc1(numClosureIndices * nleaves, &newClosurePointsC));
2451: PetscCall(PetscMalloc1(nleaves, &coverQuads));
2452: PetscCall(PetscMalloc1(cEndC - cStartC, &coverQuadsSend));
2453: count = 0;
2454: for (t = fltC; t <= lltC; t++) { /* unfortunately, we need to pack a send array, since quads are not stored packed in p4est */
2455: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2456: PetscInt q;
2458: PetscCall(PetscMemcpy(&coverQuadsSend[count], tree->quadrants.array, tree->quadrants.elem_count * sizeof(p4est_quadrant_t)));
2459: for (q = 0; (size_t)q < tree->quadrants.elem_count; q++) coverQuadsSend[count + q].p.which_tree = t;
2460: count += tree->quadrants.elem_count;
2461: }
2462: /* p is of a union type p4est_quadrant_data, but only the p.which_tree field is active at this time. So, we
2463: have a simple blockTypes[] to use. Note that quadStruct does not count potential padding in array of
2464: p4est_quadrant_t. We have to call MPI_Type_create_resized() to change upper-bound of quadStruct.
2465: */
2466: PetscCallMPI(MPI_Type_create_struct(4, blockSizes, blockOffsets, blockTypes, &quadStruct));
2467: PetscCallMPI(MPI_Type_create_resized(quadStruct, 0, sizeof(p4est_quadrant_t), &quadType));
2468: PetscCallMPI(MPI_Type_commit(&quadType));
2469: PetscCall(PetscSFBcastBegin(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2470: PetscCall(PetscSFBcastBegin(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2471: PetscCall(PetscSFBcastEnd(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE));
2472: PetscCall(PetscSFBcastEnd(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE));
2473: PetscCallMPI(MPI_Type_free(&quadStruct));
2474: PetscCallMPI(MPI_Type_free(&quadType));
2475: PetscCall(PetscFree(coverQuadsSend));
2476: PetscCall(PetscFree(closurePointsC));
2477: PetscCall(PetscSFDestroy(&coveringSF));
2478: closurePointsC = newClosurePointsC;
2480: /* assign tree quads based on locations in coverQuads */
2481: {
2482: PetscInt q;
2483: for (q = 0; q < nleaves; q++) {
2484: p4est_locidx_t t = coverQuads[q].p.which_tree;
2485: if (!treeQuadCounts[t - fltF]++) treeQuads[t - fltF] = &coverQuads[q];
2486: }
2487: }
2488: }
2489: }
2490: if (!coverQuads) { /* matching partitions: assign tree quads based on locations in p4est native arrays */
2491: for (t = fltF; t <= lltF; t++) {
2492: p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]);
2494: treeQuadCounts[t - fltF] = tree->quadrants.elem_count;
2495: treeQuads[t - fltF] = (p4est_quadrant_t *)tree->quadrants.array;
2496: }
2497: }
2499: {
2500: PetscInt p;
2501: PetscInt cLocalStartF;
2502: PetscSF pointSF;
2503: PetscSFNode *roots;
2504: PetscInt *rootType;
2505: DM refTree = NULL;
2506: DMLabel canonical;
2507: PetscInt *childClosures[P4EST_CHILDREN] = {NULL};
2508: PetscInt *rootClosure = NULL;
2509: PetscInt coarseOffset;
2510: PetscInt numCoarseQuads;
2512: PetscCall(PetscMalloc1(pEndF - pStartF, &roots));
2513: PetscCall(PetscMalloc1(pEndF - pStartF, &rootType));
2514: PetscCall(DMGetPointSF(fine, &pointSF));
2515: for (p = pStartF; p < pEndF; p++) {
2516: roots[p - pStartF].rank = -1;
2517: roots[p - pStartF].index = -1;
2518: rootType[p - pStartF] = -1;
2519: }
2520: if (formCids) {
2521: PetscInt child;
2523: PetscCall(PetscMalloc1(pEndF - pStartF, &cids));
2524: for (p = pStartF; p < pEndF; p++) cids[p - pStartF] = -2;
2525: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2526: PetscCall(DMPlexGetTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2527: for (child = 0; child < P4EST_CHILDREN; child++) { /* get the closures of the child cells in the reference tree */
2528: PetscCall(DMPlexGetTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2529: }
2530: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2531: }
2532: cLocalStartF = pforestF->cLocalStart;
2533: for (t = fltF, coarseOffset = 0, numCoarseQuads = 0; t <= lltF; t++, coarseOffset += numCoarseQuads) {
2534: p4est_tree_t *tree = &(((p4est_tree_t *)p4estF->trees->array)[t]);
2535: PetscInt numFineQuads = tree->quadrants.elem_count;
2536: p4est_quadrant_t *coarseQuads = treeQuads[t - fltF];
2537: p4est_quadrant_t *fineQuads = (p4est_quadrant_t *)tree->quadrants.array;
2538: PetscInt i, coarseCount = 0;
2539: PetscInt offset = tree->quadrants_offset;
2540: sc_array_t coarseQuadsArray;
2542: numCoarseQuads = treeQuadCounts[t - fltF];
2543: PetscCallP4est(sc_array_init_data, (&coarseQuadsArray, coarseQuads, sizeof(p4est_quadrant_t), (size_t)numCoarseQuads));
2544: for (i = 0; i < numFineQuads; i++) {
2545: PetscInt c = i + offset;
2546: p4est_quadrant_t *quad = &fineQuads[i];
2547: p4est_quadrant_t *quadCoarse = NULL;
2548: ssize_t disjoint = -1;
2550: while (disjoint < 0 && coarseCount < numCoarseQuads) {
2551: quadCoarse = &coarseQuads[coarseCount];
2552: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2553: if (disjoint < 0) coarseCount++;
2554: }
2555: PetscCheck(disjoint == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "did not find overlapping coarse quad");
2556: if (quadCoarse->level > quad->level || (quadCoarse->level == quad->level && !transferIdent)) { /* the "coarse" mesh is finer than the fine mesh at the point: continue */
2557: if (transferIdent) { /* find corners */
2558: PetscInt j = 0;
2560: do {
2561: if (j < P4EST_CHILDREN) {
2562: p4est_quadrant_t cornerQuad;
2563: int equal;
2565: PetscCallP4est(p4est_quadrant_corner_descendant, (quad, &cornerQuad, j, quadCoarse->level));
2566: PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (&cornerQuad, quadCoarse));
2567: if (equal) {
2568: PetscInt petscJ = P4estVertToPetscVert[j];
2569: PetscInt p = closurePointsF[numClosureIndices * c + (P4EST_INSUL - P4EST_CHILDREN) + petscJ].index;
2570: PetscSFNode q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + (P4EST_INSUL - P4EST_CHILDREN) + petscJ];
2572: roots[p - pStartF] = q;
2573: rootType[p - pStartF] = PETSC_MAX_INT;
2574: cids[p - pStartF] = -1;
2575: j++;
2576: }
2577: }
2578: coarseCount++;
2579: disjoint = 1;
2580: if (coarseCount < numCoarseQuads) {
2581: quadCoarse = &coarseQuads[coarseCount];
2582: PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad));
2583: }
2584: } while (!disjoint);
2585: }
2586: continue;
2587: }
2588: if (quadCoarse->level == quad->level) { /* same quad present in coarse and fine mesh */
2589: PetscInt j;
2590: for (j = 0; j < numClosureIndices; j++) {
2591: PetscInt p = closurePointsF[numClosureIndices * c + j].index;
2593: roots[p - pStartF] = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + j];
2594: rootType[p - pStartF] = PETSC_MAX_INT; /* unconditionally accept */
2595: cids[p - pStartF] = -1;
2596: }
2597: } else {
2598: PetscInt levelDiff = quad->level - quadCoarse->level;
2599: PetscInt proposedCids[P4EST_INSUL] = {0};
2601: if (formCids) {
2602: PetscInt cl;
2603: PetscInt *pointClosure = NULL;
2604: int cid;
2606: PetscCheck(levelDiff <= 1, PETSC_COMM_SELF, PETSC_ERR_USER, "Recursive child ids not implemented");
2607: PetscCallP4estReturn(cid, p4est_quadrant_child_id, (quad));
2608: PetscCall(DMPlexGetTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2609: for (cl = 0; cl < P4EST_INSUL; cl++) {
2610: PetscInt p = pointClosure[2 * cl];
2611: PetscInt point = childClosures[cid][2 * cl];
2612: PetscInt ornt = childClosures[cid][2 * cl + 1];
2613: PetscInt newcid = -1;
2614: DMPolytopeType ct;
2616: if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2617: PetscCall(DMPlexGetCellType(refTree, point, &ct));
2618: ornt = DMPolytopeConvertNewOrientation_Internal(ct, ornt);
2619: if (!cl) {
2620: newcid = cid + 1;
2621: } else {
2622: PetscInt rcl, parent, parentOrnt = 0;
2624: PetscCall(DMPlexGetTreeParent(refTree, point, &parent, NULL));
2625: if (parent == point) {
2626: newcid = -1;
2627: } else if (!parent) { /* in the root */
2628: newcid = point;
2629: } else {
2630: DMPolytopeType rct = DM_POLYTOPE_UNKNOWN;
2632: for (rcl = 1; rcl < P4EST_INSUL; rcl++) {
2633: if (rootClosure[2 * rcl] == parent) {
2634: PetscCall(DMPlexGetCellType(refTree, parent, &rct));
2635: parentOrnt = DMPolytopeConvertNewOrientation_Internal(rct, rootClosure[2 * rcl + 1]);
2636: break;
2637: }
2638: }
2639: PetscCheck(rcl < P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Couldn't find parent in root closure");
2640: PetscCall(DMPlexReferenceTreeGetChildSymmetry(refTree, parent, parentOrnt, ornt, point, DMPolytopeConvertNewOrientation_Internal(rct, pointClosure[2 * rcl + 1]), NULL, &newcid));
2641: }
2642: }
2643: if (newcid >= 0) {
2644: if (canonical) PetscCall(DMLabelGetValue(canonical, newcid, &newcid));
2645: proposedCids[cl] = newcid;
2646: }
2647: }
2648: PetscCall(DMPlexRestoreTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure));
2649: }
2650: p4est_qcoord_t coarseBound[2][P4EST_DIM] = {
2651: {quadCoarse->x, quadCoarse->y,
2652: #if defined(P4_TO_P8)
2653: quadCoarse->z
2654: #endif
2655: },
2656: {0}
2657: };
2658: p4est_qcoord_t fineBound[2][P4EST_DIM] = {
2659: {quad->x, quad->y,
2660: #if defined(P4_TO_P8)
2661: quad->z
2662: #endif
2663: },
2664: {0}
2665: };
2666: PetscInt j;
2667: for (j = 0; j < P4EST_DIM; j++) { /* get the coordinates of cell boundaries in each direction */
2668: coarseBound[1][j] = coarseBound[0][j] + P4EST_QUADRANT_LEN(quadCoarse->level);
2669: fineBound[1][j] = fineBound[0][j] + P4EST_QUADRANT_LEN(quad->level);
2670: }
2671: for (j = 0; j < numClosureIndices; j++) {
2672: PetscInt l, p;
2673: PetscSFNode q;
2675: p = closurePointsF[numClosureIndices * c + j].index;
2676: if (rootType[p - pStartF] == PETSC_MAX_INT) continue;
2677: if (j == 0) { /* volume: ancestor is volume */
2678: l = 0;
2679: } else if (j < 1 + P4EST_FACES) { /* facet */
2680: PetscInt face = PetscFaceToP4estFace[j - 1];
2681: PetscInt direction = face / 2;
2682: PetscInt coarseFace = -1;
2684: if (coarseBound[face % 2][direction] == fineBound[face % 2][direction]) {
2685: coarseFace = face;
2686: l = 1 + P4estFaceToPetscFace[coarseFace];
2687: } else {
2688: l = 0;
2689: }
2690: #if defined(P4_TO_P8)
2691: } else if (j < 1 + P4EST_FACES + P8EST_EDGES) {
2692: PetscInt edge = PetscEdgeToP4estEdge[j - (1 + P4EST_FACES)];
2693: PetscInt direction = edge / 4;
2694: PetscInt mod = edge % 4;
2695: PetscInt coarseEdge = -1, coarseFace = -1;
2696: PetscInt minDir = PetscMin((direction + 1) % 3, (direction + 2) % 3);
2697: PetscInt maxDir = PetscMax((direction + 1) % 3, (direction + 2) % 3);
2698: PetscBool dirTest[2];
2700: dirTest[0] = (PetscBool)(coarseBound[mod % 2][minDir] == fineBound[mod % 2][minDir]);
2701: dirTest[1] = (PetscBool)(coarseBound[mod / 2][maxDir] == fineBound[mod / 2][maxDir]);
2703: if (dirTest[0] && dirTest[1]) { /* fine edge falls on coarse edge */
2704: coarseEdge = edge;
2705: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2706: } else if (dirTest[0]) { /* fine edge falls on a coarse face in the minDir direction */
2707: coarseFace = 2 * minDir + (mod % 2);
2708: l = 1 + P4estFaceToPetscFace[coarseFace];
2709: } else if (dirTest[1]) { /* fine edge falls on a coarse face in the maxDir direction */
2710: coarseFace = 2 * maxDir + (mod / 2);
2711: l = 1 + P4estFaceToPetscFace[coarseFace];
2712: } else {
2713: l = 0;
2714: }
2715: #endif
2716: } else {
2717: PetscInt vertex = PetscVertToP4estVert[P4EST_CHILDREN - (P4EST_INSUL - j)];
2718: PetscBool dirTest[P4EST_DIM];
2719: PetscInt m;
2720: PetscInt numMatch = 0;
2721: PetscInt coarseVertex = -1, coarseFace = -1;
2722: #if defined(P4_TO_P8)
2723: PetscInt coarseEdge = -1;
2724: #endif
2726: for (m = 0; m < P4EST_DIM; m++) {
2727: dirTest[m] = (PetscBool)(coarseBound[(vertex >> m) & 1][m] == fineBound[(vertex >> m) & 1][m]);
2728: if (dirTest[m]) numMatch++;
2729: }
2730: if (numMatch == P4EST_DIM) { /* vertex on vertex */
2731: coarseVertex = vertex;
2732: l = P4EST_INSUL - (P4EST_CHILDREN - P4estVertToPetscVert[coarseVertex]);
2733: } else if (numMatch == 1) { /* vertex on face */
2734: for (m = 0; m < P4EST_DIM; m++) {
2735: if (dirTest[m]) {
2736: coarseFace = 2 * m + ((vertex >> m) & 1);
2737: break;
2738: }
2739: }
2740: l = 1 + P4estFaceToPetscFace[coarseFace];
2741: #if defined(P4_TO_P8)
2742: } else if (numMatch == 2) { /* vertex on edge */
2743: for (m = 0; m < P4EST_DIM; m++) {
2744: if (!dirTest[m]) {
2745: PetscInt otherDir1 = (m + 1) % 3;
2746: PetscInt otherDir2 = (m + 2) % 3;
2747: PetscInt minDir = PetscMin(otherDir1, otherDir2);
2748: PetscInt maxDir = PetscMax(otherDir1, otherDir2);
2750: coarseEdge = m * 4 + 2 * ((vertex >> maxDir) & 1) + ((vertex >> minDir) & 1);
2751: break;
2752: }
2753: }
2754: l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge];
2755: #endif
2756: } else { /* volume */
2757: l = 0;
2758: }
2759: }
2760: q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + l];
2761: if (l > rootType[p - pStartF]) {
2762: if (l >= P4EST_INSUL - P4EST_CHILDREN) { /* vertex on vertex: unconditional acceptance */
2763: if (transferIdent) {
2764: roots[p - pStartF] = q;
2765: rootType[p - pStartF] = PETSC_MAX_INT;
2766: if (formCids) cids[p - pStartF] = -1;
2767: }
2768: } else {
2769: PetscInt k, thisp = p, limit;
2771: roots[p - pStartF] = q;
2772: rootType[p - pStartF] = l;
2773: if (formCids) cids[p - pStartF] = proposedCids[j];
2774: limit = transferIdent ? levelDiff : (levelDiff - 1);
2775: for (k = 0; k < limit; k++) {
2776: PetscInt parent;
2778: PetscCall(DMPlexGetTreeParent(plexF, thisp, &parent, NULL));
2779: if (parent == thisp) break;
2781: roots[parent - pStartF] = q;
2782: rootType[parent - pStartF] = PETSC_MAX_INT;
2783: if (formCids) cids[parent - pStartF] = -1;
2784: thisp = parent;
2785: }
2786: }
2787: }
2788: }
2789: }
2790: }
2791: }
2793: /* 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 */
2794: if (size > 1) {
2795: PetscInt *rootTypeCopy, p;
2797: PetscCall(PetscMalloc1(pEndF - pStartF, &rootTypeCopy));
2798: PetscCall(PetscArraycpy(rootTypeCopy, rootType, pEndF - pStartF));
2799: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2800: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX));
2801: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2802: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE));
2803: for (p = pStartF; p < pEndF; p++) {
2804: 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 */
2805: roots[p - pStartF].rank = -1;
2806: roots[p - pStartF].index = -1;
2807: }
2808: if (formCids && rootTypeCopy[p - pStartF] == PETSC_MAX_INT) { cids[p - pStartF] = -1; /* we have found an antecedent that is the same: no child id */ }
2809: }
2810: PetscCall(PetscFree(rootTypeCopy));
2811: PetscCall(PetscSFReduceBegin(pointSF, nodeType, roots, roots, sfNodeReduce));
2812: PetscCall(PetscSFReduceEnd(pointSF, nodeType, roots, roots, sfNodeReduce));
2813: PetscCall(PetscSFBcastBegin(pointSF, nodeType, roots, roots, MPI_REPLACE));
2814: PetscCall(PetscSFBcastEnd(pointSF, nodeType, roots, roots, MPI_REPLACE));
2815: }
2816: PetscCall(PetscFree(rootType));
2818: {
2819: PetscInt numRoots;
2820: PetscInt numLeaves;
2821: PetscInt *leaves;
2822: PetscSFNode *iremote;
2823: /* count leaves */
2825: numRoots = pEndC - pStartC;
2827: numLeaves = 0;
2828: for (p = pStartF; p < pEndF; p++) {
2829: if (roots[p - pStartF].index >= 0) numLeaves++;
2830: }
2831: PetscCall(PetscMalloc1(numLeaves, &leaves));
2832: PetscCall(PetscMalloc1(numLeaves, &iremote));
2833: numLeaves = 0;
2834: for (p = pStartF; p < pEndF; p++) {
2835: if (roots[p - pStartF].index >= 0) {
2836: leaves[numLeaves] = p - pStartF;
2837: iremote[numLeaves] = roots[p - pStartF];
2838: numLeaves++;
2839: }
2840: }
2841: PetscCall(PetscFree(roots));
2842: PetscCall(PetscSFCreate(comm, sf));
2843: if (numLeaves == (pEndF - pStartF)) {
2844: PetscCall(PetscFree(leaves));
2845: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2846: } else {
2847: PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, leaves, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2848: }
2849: }
2850: if (formCids) {
2851: PetscSF pointSF;
2852: PetscInt child;
2854: PetscCall(DMPlexGetReferenceTree(plexF, &refTree));
2855: PetscCall(DMGetPointSF(plexF, &pointSF));
2856: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2857: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, cids, cids, MPI_MAX));
2858: if (childIds) *childIds = cids;
2859: for (child = 0; child < P4EST_CHILDREN; child++) PetscCall(DMPlexRestoreTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child]));
2860: PetscCall(DMPlexRestoreTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure));
2861: }
2862: }
2863: if (saveInCoarse) { /* cache results */
2864: PetscCall(PetscObjectReference((PetscObject)*sf));
2865: pforestC->pointSelfToAdaptSF = *sf;
2866: if (!childIds) {
2867: pforestC->pointSelfToAdaptCids = cids;
2868: } else {
2869: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestC->pointSelfToAdaptCids));
2870: PetscCall(PetscArraycpy(pforestC->pointSelfToAdaptCids, cids, pEndF - pStartF));
2871: }
2872: } else if (saveInFine) {
2873: PetscCall(PetscObjectReference((PetscObject)*sf));
2874: pforestF->pointAdaptToSelfSF = *sf;
2875: if (!childIds) {
2876: pforestF->pointAdaptToSelfCids = cids;
2877: } else {
2878: PetscCall(PetscMalloc1(pEndF - pStartF, &pforestF->pointAdaptToSelfCids));
2879: PetscCall(PetscArraycpy(pforestF->pointAdaptToSelfCids, cids, pEndF - pStartF));
2880: }
2881: }
2882: PetscCall(PetscFree2(treeQuads, treeQuadCounts));
2883: PetscCall(PetscFree(coverQuads));
2884: PetscCall(PetscFree(closurePointsC));
2885: PetscCall(PetscFree(closurePointsF));
2886: PetscCallMPI(MPI_Type_free(&nodeClosureType));
2887: PetscCallMPI(MPI_Op_free(&sfNodeReduce));
2888: PetscCallMPI(MPI_Type_free(&nodeType));
2889: PetscFunctionReturn(PETSC_SUCCESS);
2890: }
2892: /* children are sf leaves of parents */
2893: static PetscErrorCode DMPforestGetTransferSF_Internal(DM coarse, DM fine, const PetscInt dofPerDim[], PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[])
2894: {
2895: MPI_Comm comm;
2896: PetscMPIInt rank;
2897: DM_Forest_pforest *pforestC, *pforestF;
2898: DM plexC, plexF;
2899: PetscInt pStartC, pEndC, pStartF, pEndF;
2900: PetscSF pointTransferSF;
2901: PetscBool allOnes = PETSC_TRUE;
2903: PetscFunctionBegin;
2904: pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data;
2905: pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data;
2906: PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM");
2907: comm = PetscObjectComm((PetscObject)coarse);
2908: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2910: {
2911: PetscInt i;
2912: for (i = 0; i <= P4EST_DIM; i++) {
2913: if (dofPerDim[i] != 1) {
2914: allOnes = PETSC_FALSE;
2915: break;
2916: }
2917: }
2918: }
2919: PetscCall(DMPforestGetTransferSF_Point(coarse, fine, &pointTransferSF, transferIdent, childIds));
2920: if (allOnes) {
2921: *sf = pointTransferSF;
2922: PetscFunctionReturn(PETSC_SUCCESS);
2923: }
2925: PetscCall(DMPforestGetPlex(fine, &plexF));
2926: PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF));
2927: PetscCall(DMPforestGetPlex(coarse, &plexC));
2928: PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC));
2929: {
2930: PetscInt numRoots;
2931: PetscInt numLeaves;
2932: const PetscInt *leaves;
2933: const PetscSFNode *iremote;
2934: PetscInt d;
2935: PetscSection leafSection, rootSection;
2937: /* count leaves */
2938: PetscCall(PetscSFGetGraph(pointTransferSF, &numRoots, &numLeaves, &leaves, &iremote));
2939: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &rootSection));
2940: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &leafSection));
2941: PetscCall(PetscSectionSetChart(rootSection, pStartC, pEndC));
2942: PetscCall(PetscSectionSetChart(leafSection, pStartF, pEndF));
2944: for (d = 0; d <= P4EST_DIM; d++) {
2945: PetscInt startC, endC, e;
2947: PetscCall(DMPlexGetSimplexOrBoxCells(plexC, P4EST_DIM - d, &startC, &endC));
2948: for (e = startC; e < endC; e++) PetscCall(PetscSectionSetDof(rootSection, e, dofPerDim[d]));
2949: }
2951: for (d = 0; d <= P4EST_DIM; d++) {
2952: PetscInt startF, endF, e;
2954: PetscCall(DMPlexGetSimplexOrBoxCells(plexF, P4EST_DIM - d, &startF, &endF));
2955: for (e = startF; e < endF; e++) PetscCall(PetscSectionSetDof(leafSection, e, dofPerDim[d]));
2956: }
2958: PetscCall(PetscSectionSetUp(rootSection));
2959: PetscCall(PetscSectionSetUp(leafSection));
2960: {
2961: PetscInt nroots, nleaves;
2962: PetscInt *mine, i, p;
2963: PetscInt *offsets, *offsetsRoot;
2964: PetscSFNode *remote;
2966: PetscCall(PetscMalloc1(pEndF - pStartF, &offsets));
2967: PetscCall(PetscMalloc1(pEndC - pStartC, &offsetsRoot));
2968: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionGetOffset(rootSection, p, &offsetsRoot[p - pStartC]));
2969: PetscCall(PetscSFBcastBegin(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2970: PetscCall(PetscSFBcastEnd(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE));
2971: PetscCall(PetscSectionGetStorageSize(rootSection, &nroots));
2972: nleaves = 0;
2973: for (i = 0; i < numLeaves; i++) {
2974: PetscInt leaf = leaves ? leaves[i] : i;
2975: PetscInt dof;
2977: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2978: nleaves += dof;
2979: }
2980: PetscCall(PetscMalloc1(nleaves, &mine));
2981: PetscCall(PetscMalloc1(nleaves, &remote));
2982: nleaves = 0;
2983: for (i = 0; i < numLeaves; i++) {
2984: PetscInt leaf = leaves ? leaves[i] : i;
2985: PetscInt dof;
2986: PetscInt off, j;
2988: PetscCall(PetscSectionGetDof(leafSection, leaf, &dof));
2989: PetscCall(PetscSectionGetOffset(leafSection, leaf, &off));
2990: for (j = 0; j < dof; j++) {
2991: remote[nleaves].rank = iremote[i].rank;
2992: remote[nleaves].index = offsets[leaf] + j;
2993: mine[nleaves++] = off + j;
2994: }
2995: }
2996: PetscCall(PetscFree(offsetsRoot));
2997: PetscCall(PetscFree(offsets));
2998: PetscCall(PetscSFCreate(comm, sf));
2999: PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
3000: }
3001: PetscCall(PetscSectionDestroy(&leafSection));
3002: PetscCall(PetscSectionDestroy(&rootSection));
3003: PetscCall(PetscSFDestroy(&pointTransferSF));
3004: }
3005: PetscFunctionReturn(PETSC_SUCCESS);
3006: }
3008: static PetscErrorCode DMPforestGetTransferSF(DM dmA, DM dmB, const PetscInt dofPerDim[], PetscSF *sfAtoB, PetscSF *sfBtoA)
3009: {
3010: DM adaptA, adaptB;
3011: DMAdaptFlag purpose;
3013: PetscFunctionBegin;
3014: PetscCall(DMForestGetAdaptivityForest(dmA, &adaptA));
3015: PetscCall(DMForestGetAdaptivityForest(dmB, &adaptB));
3016: /* it is more efficient when the coarser mesh is the first argument: reorder if we know one is coarser than the other */
3017: if (adaptA && adaptA->data == dmB->data) { /* dmA was adapted from dmB */
3018: PetscCall(DMForestGetAdaptivityPurpose(dmA, &purpose));
3019: if (purpose == DM_ADAPT_REFINE) {
3020: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3021: PetscFunctionReturn(PETSC_SUCCESS);
3022: }
3023: } else if (adaptB && adaptB->data == dmA->data) { /* dmB was adapted from dmA */
3024: PetscCall(DMForestGetAdaptivityPurpose(dmB, &purpose));
3025: if (purpose == DM_ADAPT_COARSEN) {
3026: PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB));
3027: PetscFunctionReturn(PETSC_SUCCESS);
3028: }
3029: }
3030: if (sfAtoB) PetscCall(DMPforestGetTransferSF_Internal(dmA, dmB, dofPerDim, sfAtoB, PETSC_TRUE, NULL));
3031: if (sfBtoA) PetscCall(DMPforestGetTransferSF_Internal(dmB, dmA, dofPerDim, sfBtoA, (PetscBool)(sfAtoB == NULL), NULL));
3032: PetscFunctionReturn(PETSC_SUCCESS);
3033: }
3035: static PetscErrorCode DMPforestLabelsInitialize(DM dm, DM plex)
3036: {
3037: DM_Forest *forest = (DM_Forest *)dm->data;
3038: DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data;
3039: PetscInt cLocalStart, cLocalEnd, cStart, cEnd, fStart, fEnd, eStart, eEnd, vStart, vEnd;
3040: PetscInt cStartBase, cEndBase, fStartBase, fEndBase, vStartBase, vEndBase, eStartBase, eEndBase;
3041: PetscInt pStart, pEnd, pStartBase, pEndBase, p;
3042: DM base;
3043: PetscInt *star = NULL, starSize;
3044: DMLabelLink next = dm->labels;
3045: PetscInt guess = 0;
3046: p4est_topidx_t num_trees = pforest->topo->conn->num_trees;
3048: PetscFunctionBegin;
3049: pforest->labelsFinalized = PETSC_TRUE;
3050: cLocalStart = pforest->cLocalStart;
3051: cLocalEnd = pforest->cLocalEnd;
3052: PetscCall(DMForestGetBaseDM(dm, &base));
3053: if (!base) {
3054: if (pforest->ghostName) { /* insert a label to make the boundaries, with stratum values denoting which face of the element touches the boundary */
3055: p4est_connectivity_t *conn = pforest->topo->conn;
3056: p4est_t *p4est = pforest->forest;
3057: p4est_tree_t *trees = (p4est_tree_t *)p4est->trees->array;
3058: p4est_topidx_t t, flt = p4est->first_local_tree;
3059: p4est_topidx_t llt = pforest->forest->last_local_tree;
3060: DMLabel ghostLabel;
3061: PetscInt c;
3063: PetscCall(DMCreateLabel(plex, pforest->ghostName));
3064: PetscCall(DMGetLabel(plex, pforest->ghostName, &ghostLabel));
3065: for (c = cLocalStart, t = flt; t <= llt; t++) {
3066: p4est_tree_t *tree = &trees[t];
3067: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3068: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3069: PetscInt q;
3071: for (q = 0; q < numQuads; q++, c++) {
3072: p4est_quadrant_t *quad = &quads[q];
3073: PetscInt f;
3075: for (f = 0; f < P4EST_FACES; f++) {
3076: p4est_quadrant_t neigh;
3077: int isOutside;
3079: PetscCallP4est(p4est_quadrant_face_neighbor, (quad, f, &neigh));
3080: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&neigh));
3081: if (isOutside) {
3082: p4est_topidx_t nt;
3083: PetscInt nf;
3085: nt = conn->tree_to_tree[t * P4EST_FACES + f];
3086: nf = (PetscInt)conn->tree_to_face[t * P4EST_FACES + f];
3087: nf = nf % P4EST_FACES;
3088: if (nt == t && nf == f) {
3089: PetscInt plexF = P4estFaceToPetscFace[f];
3090: const PetscInt *cone;
3092: PetscCall(DMPlexGetCone(plex, c, &cone));
3093: PetscCall(DMLabelSetValue(ghostLabel, cone[plexF], plexF + 1));
3094: }
3095: }
3096: }
3097: }
3098: }
3099: }
3100: PetscFunctionReturn(PETSC_SUCCESS);
3101: }
3102: PetscCall(DMPlexGetSimplexOrBoxCells(base, 0, &cStartBase, &cEndBase));
3103: PetscCall(DMPlexGetSimplexOrBoxCells(base, 1, &fStartBase, &fEndBase));
3104: PetscCall(DMPlexGetSimplexOrBoxCells(base, P4EST_DIM - 1, &eStartBase, &eEndBase));
3105: PetscCall(DMPlexGetDepthStratum(base, 0, &vStartBase, &vEndBase));
3107: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
3108: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 1, &fStart, &fEnd));
3109: PetscCall(DMPlexGetSimplexOrBoxCells(plex, P4EST_DIM - 1, &eStart, &eEnd));
3110: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3112: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3113: PetscCall(DMPlexGetChart(base, &pStartBase, &pEndBase));
3114: /* go through the mesh: use star to find a quadrant that borders a point. Use the closure to determine the
3115: * orientation of the quadrant relative to that point. Use that to relate the point to the numbering in the base
3116: * mesh, and extract a label value (since the base mesh is redundantly distributed, can be found locally). */
3117: while (next) {
3118: DMLabel baseLabel;
3119: DMLabel label = next->label;
3120: PetscBool isDepth, isCellType, isGhost, isVTK, isSpmap;
3121: const char *name;
3123: PetscCall(PetscObjectGetName((PetscObject)label, &name));
3124: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3125: if (isDepth) {
3126: next = next->next;
3127: continue;
3128: }
3129: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3130: if (isCellType) {
3131: next = next->next;
3132: continue;
3133: }
3134: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3135: if (isGhost) {
3136: next = next->next;
3137: continue;
3138: }
3139: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3140: if (isVTK) {
3141: next = next->next;
3142: continue;
3143: }
3144: PetscCall(PetscStrcmp(name, "_forest_base_subpoint_map", &isSpmap));
3145: if (!isSpmap) {
3146: PetscCall(DMGetLabel(base, name, &baseLabel));
3147: if (!baseLabel) {
3148: next = next->next;
3149: continue;
3150: }
3151: PetscCall(DMLabelCreateIndex(baseLabel, pStartBase, pEndBase));
3152: } else baseLabel = NULL;
3154: for (p = pStart; p < pEnd; p++) {
3155: PetscInt s, c = -1, l;
3156: PetscInt *closure = NULL, closureSize;
3157: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3158: p4est_tree_t *trees = (p4est_tree_t *)pforest->forest->trees->array;
3159: p4est_quadrant_t *q;
3160: PetscInt t, val;
3161: PetscBool zerosupportpoint = PETSC_FALSE;
3163: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3164: for (s = 0; s < starSize; s++) {
3165: PetscInt point = star[2 * s];
3167: if (cStart <= point && point < cEnd) {
3168: PetscCall(DMPlexGetTransitiveClosure(plex, point, PETSC_TRUE, &closureSize, &closure));
3169: for (l = 0; l < closureSize; l++) {
3170: PetscInt qParent = closure[2 * l], q, pp = p, pParent = p;
3171: do { /* check parents of q */
3172: q = qParent;
3173: if (q == p) {
3174: c = point;
3175: break;
3176: }
3177: PetscCall(DMPlexGetTreeParent(plex, q, &qParent, NULL));
3178: } while (qParent != q);
3179: if (c != -1) break;
3180: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3181: q = closure[2 * l];
3182: while (pParent != pp) { /* check parents of p */
3183: pp = pParent;
3184: if (pp == q) {
3185: c = point;
3186: break;
3187: }
3188: PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL));
3189: }
3190: if (c != -1) break;
3191: }
3192: PetscCall(DMPlexRestoreTransitiveClosure(plex, point, PETSC_TRUE, NULL, &closure));
3193: if (l < closureSize) break;
3194: } else {
3195: PetscInt supportSize;
3197: PetscCall(DMPlexGetSupportSize(plex, point, &supportSize));
3198: zerosupportpoint = (PetscBool)(zerosupportpoint || !supportSize);
3199: }
3200: }
3201: if (c < 0) {
3202: const char *prefix;
3203: PetscBool print = PETSC_FALSE;
3205: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
3206: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_forest_print_label_error", &print, NULL));
3207: if (print) {
3208: PetscInt i;
3210: 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));
3211: 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]));
3212: }
3213: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3214: if (zerosupportpoint) continue;
3215: else
3216: 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");
3217: }
3218: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star));
3220: if (c < cLocalStart) {
3221: /* get from the beginning of the ghost layer */
3222: q = &ghosts[c];
3223: t = (PetscInt)q->p.which_tree;
3224: } else if (c < cLocalEnd) {
3225: PetscInt lo = 0, hi = num_trees;
3226: /* get from local quadrants: have to find the right tree */
3228: c -= cLocalStart;
3230: do {
3231: p4est_tree_t *tree;
3233: PetscCheck(guess >= lo && guess < num_trees && lo < hi, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed binary search");
3234: tree = &trees[guess];
3235: if (c < tree->quadrants_offset) {
3236: hi = guess;
3237: } else if (c < tree->quadrants_offset + (PetscInt)tree->quadrants.elem_count) {
3238: q = &((p4est_quadrant_t *)tree->quadrants.array)[c - (PetscInt)tree->quadrants_offset];
3239: t = guess;
3240: break;
3241: } else {
3242: lo = guess + 1;
3243: }
3244: guess = lo + (hi - lo) / 2;
3245: } while (1);
3246: } else {
3247: /* get from the end of the ghost layer */
3248: c -= (cLocalEnd - cLocalStart);
3250: q = &ghosts[c];
3251: t = (PetscInt)q->p.which_tree;
3252: }
3254: if (l == 0) { /* cell */
3255: if (baseLabel) {
3256: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3257: } else {
3258: val = t + cStartBase;
3259: }
3260: PetscCall(DMLabelSetValue(label, p, val));
3261: } else if (l >= 1 && l < 1 + P4EST_FACES) { /* facet */
3262: p4est_quadrant_t nq;
3263: int isInside;
3265: l = PetscFaceToP4estFace[l - 1];
3266: PetscCallP4est(p4est_quadrant_face_neighbor, (q, l, &nq));
3267: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3268: if (isInside) {
3269: /* this facet is in the interior of a tree, so it inherits the label of the tree */
3270: if (baseLabel) {
3271: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3272: } else {
3273: val = t + cStartBase;
3274: }
3275: PetscCall(DMLabelSetValue(label, p, val));
3276: } else {
3277: PetscInt f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + l];
3279: if (baseLabel) {
3280: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3281: } else {
3282: val = f + fStartBase;
3283: }
3284: PetscCall(DMLabelSetValue(label, p, val));
3285: }
3286: #if defined(P4_TO_P8)
3287: } else if (l >= 1 + P4EST_FACES && l < 1 + P4EST_FACES + P8EST_EDGES) { /* edge */
3288: p4est_quadrant_t nq;
3289: int isInside;
3291: l = PetscEdgeToP4estEdge[l - (1 + P4EST_FACES)];
3292: PetscCallP4est(p8est_quadrant_edge_neighbor, (q, l, &nq));
3293: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3294: if (isInside) {
3295: /* this edge is in the interior of a tree, so it inherits the label of the tree */
3296: if (baseLabel) {
3297: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3298: } else {
3299: val = t + cStartBase;
3300: }
3301: PetscCall(DMLabelSetValue(label, p, val));
3302: } else {
3303: int isOutsideFace;
3305: PetscCallP4estReturn(isOutsideFace, p4est_quadrant_is_outside_face, (&nq));
3306: if (isOutsideFace) {
3307: PetscInt f;
3309: if (nq.x < 0) {
3310: f = 0;
3311: } else if (nq.x >= P4EST_ROOT_LEN) {
3312: f = 1;
3313: } else if (nq.y < 0) {
3314: f = 2;
3315: } else if (nq.y >= P4EST_ROOT_LEN) {
3316: f = 3;
3317: } else if (nq.z < 0) {
3318: f = 4;
3319: } else {
3320: f = 5;
3321: }
3322: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3323: if (baseLabel) {
3324: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3325: } else {
3326: val = f + fStartBase;
3327: }
3328: PetscCall(DMLabelSetValue(label, p, val));
3329: } else { /* the quadrant edge corresponds to the tree edge */
3330: PetscInt e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + l];
3332: if (baseLabel) {
3333: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3334: } else {
3335: val = e + eStartBase;
3336: }
3337: PetscCall(DMLabelSetValue(label, p, val));
3338: }
3339: }
3340: #endif
3341: } else { /* vertex */
3342: p4est_quadrant_t nq;
3343: int isInside;
3345: #if defined(P4_TO_P8)
3346: l = PetscVertToP4estVert[l - (1 + P4EST_FACES + P8EST_EDGES)];
3347: #else
3348: l = PetscVertToP4estVert[l - (1 + P4EST_FACES)];
3349: #endif
3350: PetscCallP4est(p4est_quadrant_corner_neighbor, (q, l, &nq));
3351: PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq));
3352: if (isInside) {
3353: if (baseLabel) {
3354: PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val));
3355: } else {
3356: val = t + cStartBase;
3357: }
3358: PetscCall(DMLabelSetValue(label, p, val));
3359: } else {
3360: int isOutside;
3362: PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&nq));
3363: if (isOutside) {
3364: PetscInt f = -1;
3366: if (nq.x < 0) {
3367: f = 0;
3368: } else if (nq.x >= P4EST_ROOT_LEN) {
3369: f = 1;
3370: } else if (nq.y < 0) {
3371: f = 2;
3372: } else if (nq.y >= P4EST_ROOT_LEN) {
3373: f = 3;
3374: #if defined(P4_TO_P8)
3375: } else if (nq.z < 0) {
3376: f = 4;
3377: } else {
3378: f = 5;
3379: #endif
3380: }
3381: f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f];
3382: if (baseLabel) {
3383: PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val));
3384: } else {
3385: val = f + fStartBase;
3386: }
3387: PetscCall(DMLabelSetValue(label, p, val));
3388: continue;
3389: }
3390: #if defined(P4_TO_P8)
3391: PetscCallP4estReturn(isOutside, p8est_quadrant_is_outside_edge, (&nq));
3392: if (isOutside) {
3393: /* outside edge */
3394: PetscInt e = -1;
3396: if (nq.x >= 0 && nq.x < P4EST_ROOT_LEN) {
3397: if (nq.z < 0) {
3398: if (nq.y < 0) {
3399: e = 0;
3400: } else {
3401: e = 1;
3402: }
3403: } else {
3404: if (nq.y < 0) {
3405: e = 2;
3406: } else {
3407: e = 3;
3408: }
3409: }
3410: } else if (nq.y >= 0 && nq.y < P4EST_ROOT_LEN) {
3411: if (nq.z < 0) {
3412: if (nq.x < 0) {
3413: e = 4;
3414: } else {
3415: e = 5;
3416: }
3417: } else {
3418: if (nq.x < 0) {
3419: e = 6;
3420: } else {
3421: e = 7;
3422: }
3423: }
3424: } else {
3425: if (nq.y < 0) {
3426: if (nq.x < 0) {
3427: e = 8;
3428: } else {
3429: e = 9;
3430: }
3431: } else {
3432: if (nq.x < 0) {
3433: e = 10;
3434: } else {
3435: e = 11;
3436: }
3437: }
3438: }
3440: e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + e];
3441: if (baseLabel) {
3442: PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val));
3443: } else {
3444: val = e + eStartBase;
3445: }
3446: PetscCall(DMLabelSetValue(label, p, val));
3447: continue;
3448: }
3449: #endif
3450: {
3451: /* outside vertex: same corner as quadrant corner */
3452: PetscInt v = pforest->topo->conn->tree_to_corner[P4EST_CHILDREN * t + l];
3454: if (baseLabel) {
3455: PetscCall(DMLabelGetValue(baseLabel, v + vStartBase, &val));
3456: } else {
3457: val = v + vStartBase;
3458: }
3459: PetscCall(DMLabelSetValue(label, p, val));
3460: }
3461: }
3462: }
3463: }
3464: next = next->next;
3465: }
3466: PetscFunctionReturn(PETSC_SUCCESS);
3467: }
3469: static PetscErrorCode DMPforestLabelsFinalize(DM dm, DM plex)
3470: {
3471: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
3472: DM adapt;
3474: PetscFunctionBegin;
3475: if (pforest->labelsFinalized) PetscFunctionReturn(PETSC_SUCCESS);
3476: pforest->labelsFinalized = PETSC_TRUE;
3477: PetscCall(DMForestGetAdaptivityForest(dm, &adapt));
3478: if (!adapt) {
3479: /* Initialize labels from the base dm */
3480: PetscCall(DMPforestLabelsInitialize(dm, plex));
3481: } else {
3482: PetscInt dofPerDim[4] = {1, 1, 1, 1};
3483: PetscSF transferForward, transferBackward, pointSF;
3484: PetscInt pStart, pEnd, pStartA, pEndA;
3485: PetscInt *values, *adaptValues;
3486: DMLabelLink next = adapt->labels;
3487: DMLabel adaptLabel;
3488: DM adaptPlex;
3490: PetscCall(DMForestGetAdaptivityLabel(dm, &adaptLabel));
3491: PetscCall(DMPforestGetPlex(adapt, &adaptPlex));
3492: PetscCall(DMPforestGetTransferSF(adapt, dm, dofPerDim, &transferForward, &transferBackward));
3493: PetscCall(DMPlexGetChart(plex, &pStart, &pEnd));
3494: PetscCall(DMPlexGetChart(adaptPlex, &pStartA, &pEndA));
3495: PetscCall(PetscMalloc2(pEnd - pStart, &values, pEndA - pStartA, &adaptValues));
3496: PetscCall(DMGetPointSF(plex, &pointSF));
3497: if (PetscDefined(USE_DEBUG)) {
3498: PetscInt p;
3499: for (p = pStartA; p < pEndA; p++) adaptValues[p - pStartA] = -1;
3500: for (p = pStart; p < pEnd; p++) values[p - pStart] = -2;
3501: if (transferForward) {
3502: PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3503: PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3504: }
3505: if (transferBackward) {
3506: PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3507: PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3508: }
3509: for (p = pStart; p < pEnd; p++) {
3510: PetscInt q = p, parent;
3512: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3513: while (parent != q) {
3514: if (values[parent] == -2) values[parent] = values[q];
3515: q = parent;
3516: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3517: }
3518: }
3519: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3520: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3521: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3522: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3523: for (p = pStart; p < pEnd; p++) PetscCheck(values[p - pStart] != -2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "uncovered point %" PetscInt_FMT, p);
3524: }
3525: while (next) {
3526: DMLabel nextLabel = next->label;
3527: const char *name;
3528: PetscBool isDepth, isCellType, isGhost, isVTK;
3529: DMLabel label;
3530: PetscInt p;
3532: PetscCall(PetscObjectGetName((PetscObject)nextLabel, &name));
3533: PetscCall(PetscStrcmp(name, "depth", &isDepth));
3534: if (isDepth) {
3535: next = next->next;
3536: continue;
3537: }
3538: PetscCall(PetscStrcmp(name, "celltype", &isCellType));
3539: if (isCellType) {
3540: next = next->next;
3541: continue;
3542: }
3543: PetscCall(PetscStrcmp(name, "ghost", &isGhost));
3544: if (isGhost) {
3545: next = next->next;
3546: continue;
3547: }
3548: PetscCall(PetscStrcmp(name, "vtk", &isVTK));
3549: if (isVTK) {
3550: next = next->next;
3551: continue;
3552: }
3553: if (nextLabel == adaptLabel) {
3554: next = next->next;
3555: continue;
3556: }
3557: /* label was created earlier */
3558: PetscCall(DMGetLabel(dm, name, &label));
3559: for (p = pStartA; p < pEndA; p++) PetscCall(DMLabelGetValue(nextLabel, p, &adaptValues[p]));
3560: for (p = pStart; p < pEnd; p++) values[p] = PETSC_MIN_INT;
3562: if (transferForward) PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3563: if (transferBackward) PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3564: if (transferForward) PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE));
3565: if (transferBackward) PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX));
3566: for (p = pStart; p < pEnd; p++) {
3567: PetscInt q = p, parent;
3569: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3570: while (parent != q) {
3571: if (values[parent] == PETSC_MIN_INT) values[parent] = values[q];
3572: q = parent;
3573: PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL));
3574: }
3575: }
3576: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX));
3577: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX));
3578: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3579: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE));
3581: for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(label, p, values[p]));
3582: next = next->next;
3583: }
3584: PetscCall(PetscFree2(values, adaptValues));
3585: PetscCall(PetscSFDestroy(&transferForward));
3586: PetscCall(PetscSFDestroy(&transferBackward));
3587: pforest->labelsFinalized = PETSC_TRUE;
3588: }
3589: PetscFunctionReturn(PETSC_SUCCESS);
3590: }
3592: 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)
3593: {
3594: PetscInt closureSize, c, coordStart, coordEnd, coordDim;
3595: PetscInt *closure = NULL;
3596: PetscSection coordSec;
3598: PetscFunctionBegin;
3599: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3600: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3601: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3602: PetscCall(DMPlexGetTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3603: for (c = 0; c < closureSize; c++) {
3604: PetscInt point = closure[2 * c];
3606: if (point >= coordStart && point < coordEnd) {
3607: PetscInt dof, off;
3608: PetscInt nCoords, i;
3609: PetscCall(PetscSectionGetDof(coordSec, point, &dof));
3610: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3611: nCoords = dof / coordDim;
3612: PetscCall(PetscSectionGetOffset(coordSec, point, &off));
3613: for (i = 0; i < nCoords; i++) {
3614: PetscScalar *coord = &coords[off + i * coordDim];
3615: double coordP4est[3] = {0.};
3616: double coordP4estMapped[3] = {0.};
3617: PetscInt j;
3618: PetscReal treeCoords[P4EST_CHILDREN][3] = {{0.}};
3619: PetscReal eta[3] = {0.};
3620: PetscInt numRounds = 10;
3621: PetscReal coordGuess[3] = {0.};
3623: eta[0] = (PetscReal)q->x / (PetscReal)P4EST_ROOT_LEN;
3624: eta[1] = (PetscReal)q->y / (PetscReal)P4EST_ROOT_LEN;
3625: #if defined(P4_TO_P8)
3626: eta[2] = (PetscReal)q->z / (PetscReal)P4EST_ROOT_LEN;
3627: #endif
3629: for (j = 0; j < P4EST_CHILDREN; j++) {
3630: PetscInt k;
3632: for (k = 0; k < 3; k++) treeCoords[j][k] = conn->vertices[3 * conn->tree_to_vertex[P4EST_CHILDREN * t + j] + k];
3633: }
3635: for (j = 0; j < P4EST_CHILDREN; j++) {
3636: PetscInt k;
3637: PetscReal prod = 1.;
3639: for (k = 0; k < P4EST_DIM; k++) prod *= (j & (1 << k)) ? eta[k] : (1. - eta[k]);
3640: for (k = 0; k < 3; k++) coordGuess[k] += prod * treeCoords[j][k];
3641: }
3643: for (j = 0; j < numRounds; j++) {
3644: PetscInt dir;
3646: for (dir = 0; dir < P4EST_DIM; dir++) {
3647: PetscInt k;
3648: PetscReal diff[3];
3649: PetscReal dXdeta[3] = {0.};
3650: PetscReal rhs, scale, update;
3652: for (k = 0; k < 3; k++) diff[k] = coordP4est[k] - coordGuess[k];
3653: for (k = 0; k < P4EST_CHILDREN; k++) {
3654: PetscInt l;
3655: PetscReal prod = 1.;
3657: for (l = 0; l < P4EST_DIM; l++) {
3658: if (l == dir) {
3659: prod *= (k & (1 << l)) ? 1. : -1.;
3660: } else {
3661: prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3662: }
3663: }
3664: for (l = 0; l < 3; l++) dXdeta[l] += prod * treeCoords[k][l];
3665: }
3666: rhs = 0.;
3667: scale = 0;
3668: for (k = 0; k < 3; k++) {
3669: rhs += diff[k] * dXdeta[k];
3670: scale += dXdeta[k] * dXdeta[k];
3671: }
3672: update = rhs / scale;
3673: eta[dir] += update;
3674: eta[dir] = PetscMin(eta[dir], 1.);
3675: eta[dir] = PetscMax(eta[dir], 0.);
3677: coordGuess[0] = coordGuess[1] = coordGuess[2] = 0.;
3678: for (k = 0; k < P4EST_CHILDREN; k++) {
3679: PetscInt l;
3680: PetscReal prod = 1.;
3682: for (l = 0; l < P4EST_DIM; l++) prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]);
3683: for (l = 0; l < 3; l++) coordGuess[l] += prod * treeCoords[k][l];
3684: }
3685: }
3686: }
3687: for (j = 0; j < 3; j++) coordP4est[j] = (double)eta[j];
3689: if (geom) {
3690: (geom->X)(geom, t, coordP4est, coordP4estMapped);
3691: for (j = 0; j < coordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3692: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded");
3693: }
3694: }
3695: }
3696: PetscCall(DMPlexRestoreTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure));
3697: PetscFunctionReturn(PETSC_SUCCESS);
3698: }
3700: static PetscErrorCode DMPforestMapCoordinates(DM dm, DM plex)
3701: {
3702: DM_Forest *forest;
3703: DM_Forest_pforest *pforest;
3704: p4est_geometry_t *geom;
3705: PetscInt cLocalStart, cLocalEnd;
3706: Vec coordLocalVec;
3707: PetscScalar *coords;
3708: p4est_topidx_t flt, llt, t;
3709: p4est_tree_t *trees;
3710: PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *);
3711: void *mapCtx;
3713: PetscFunctionBegin;
3714: forest = (DM_Forest *)dm->data;
3715: pforest = (DM_Forest_pforest *)forest->data;
3716: geom = pforest->topo->geom;
3717: PetscCall(DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx));
3718: if (!geom && !map) PetscFunctionReturn(PETSC_SUCCESS);
3719: PetscCall(DMGetCoordinatesLocal(plex, &coordLocalVec));
3720: PetscCall(VecGetArray(coordLocalVec, &coords));
3721: cLocalStart = pforest->cLocalStart;
3722: cLocalEnd = pforest->cLocalEnd;
3723: flt = pforest->forest->first_local_tree;
3724: llt = pforest->forest->last_local_tree;
3725: trees = (p4est_tree_t *)pforest->forest->trees->array;
3726: if (map) { /* apply the map directly to the existing coordinates */
3727: PetscSection coordSec;
3728: PetscInt coordStart, coordEnd, p, coordDim, p4estCoordDim, cStart, cEnd, cEndInterior;
3729: DM base;
3731: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3732: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3733: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3734: PetscCall(DMForestGetBaseDM(dm, &base));
3735: PetscCall(DMGetCoordinateSection(plex, &coordSec));
3736: PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd));
3737: PetscCall(DMGetCoordinateDim(plex, &coordDim));
3738: p4estCoordDim = PetscMin(coordDim, 3);
3739: for (p = coordStart; p < coordEnd; p++) {
3740: PetscInt *star = NULL, starSize;
3741: PetscInt dof, off, cell = -1, coarsePoint = -1;
3742: PetscInt nCoords, i;
3743: PetscCall(PetscSectionGetDof(coordSec, p, &dof));
3744: PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout");
3745: nCoords = dof / coordDim;
3746: PetscCall(PetscSectionGetOffset(coordSec, p, &off));
3747: PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3748: for (i = 0; i < starSize; i++) {
3749: PetscInt point = star[2 * i];
3751: if (cStart <= point && point < cEnd) {
3752: cell = point;
3753: break;
3754: }
3755: }
3756: PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star));
3757: if (cell >= 0) {
3758: if (cell < cLocalStart) {
3759: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3761: coarsePoint = ghosts[cell].p.which_tree;
3762: } else if (cell < cLocalEnd) {
3763: cell -= cLocalStart;
3764: for (t = flt; t <= llt; t++) {
3765: p4est_tree_t *tree = &trees[t];
3767: if (cell >= tree->quadrants_offset && (size_t)cell < tree->quadrants_offset + tree->quadrants.elem_count) {
3768: coarsePoint = t;
3769: break;
3770: }
3771: }
3772: } else {
3773: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3775: coarsePoint = ghosts[cell - cLocalEnd].p.which_tree;
3776: }
3777: }
3778: for (i = 0; i < nCoords; i++) {
3779: PetscScalar *coord = &coords[off + i * coordDim];
3780: PetscReal coordP4est[3] = {0.};
3781: PetscReal coordP4estMapped[3] = {0.};
3782: PetscInt j;
3784: for (j = 0; j < p4estCoordDim; j++) coordP4est[j] = PetscRealPart(coord[j]);
3785: PetscCall((map)(base, coarsePoint, p4estCoordDim, coordP4est, coordP4estMapped, mapCtx));
3786: for (j = 0; j < p4estCoordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j];
3787: }
3788: }
3789: } else { /* we have to transform coordinates back to the unit cube (where geom is defined), and then apply geom */
3790: PetscInt cStart, cEnd, cEndInterior;
3792: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3793: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3794: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3795: if (cLocalStart > 0) {
3796: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3797: PetscInt count;
3799: for (count = 0; count < cLocalStart; count++) {
3800: p4est_quadrant_t *quad = &ghosts[count];
3801: p4est_topidx_t t = quad->p.which_tree;
3803: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, quad, t, pforest->topo->conn, coords));
3804: }
3805: }
3806: for (t = flt; t <= llt; t++) {
3807: p4est_tree_t *tree = &trees[t];
3808: PetscInt offset = cLocalStart + tree->quadrants_offset, i;
3809: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
3810: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
3812: for (i = 0; i < numQuads; i++) {
3813: PetscInt count = i + offset;
3815: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, &quads[i], t, pforest->topo->conn, coords));
3816: }
3817: }
3818: if (cLocalEnd - cLocalStart < cEnd - cStart) {
3819: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3820: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
3821: PetscInt count;
3823: for (count = 0; count < numGhosts - cLocalStart; count++) {
3824: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
3825: p4est_topidx_t t = quad->p.which_tree;
3827: PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count + cLocalEnd, quad, t, pforest->topo->conn, coords));
3828: }
3829: }
3830: }
3831: PetscCall(VecRestoreArray(coordLocalVec, &coords));
3832: PetscFunctionReturn(PETSC_SUCCESS);
3833: }
3835: static PetscErrorCode PforestQuadrantIsInterior(p4est_quadrant_t *quad, PetscBool *is_interior)
3836: {
3837: PetscFunctionBegin;
3838: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3839: if ((quad->x > 0) && (quad->x + h < P4EST_ROOT_LEN)
3840: #ifdef P4_TO_P8
3841: && (quad->z > 0) && (quad->z + h < P4EST_ROOT_LEN)
3842: #endif
3843: && (quad->y > 0) && (quad->y + h < P4EST_ROOT_LEN)) {
3844: *is_interior = PETSC_TRUE;
3845: } else {
3846: *is_interior = PETSC_FALSE;
3847: }
3848: PetscFunctionReturn(PETSC_SUCCESS);
3849: }
3851: /* We always use DG coordinates with p4est: if they do not match the vertex
3852: coordinates, add space for them in the section */
3853: static PetscErrorCode PforestCheckLocalizeCell(DM plex, PetscInt cDim, Vec cVecOld, DM_Forest_pforest *pforest, PetscSection oldSection, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad)
3854: {
3855: PetscBool is_interior;
3857: PetscFunctionBegin;
3858: PetscCall(PforestQuadrantIsInterior(quad, &is_interior));
3859: if (is_interior) { // quads in the interior of a coarse cell can't touch periodic interfaces
3860: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3861: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3862: } else {
3863: PetscInt cSize;
3864: PetscScalar *values = NULL;
3865: PetscBool same_coords = PETSC_TRUE;
3867: PetscCall(DMPlexVecGetClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3868: PetscAssert(cSize == cDim * P4EST_CHILDREN, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected closure size");
3869: for (int c = 0; c < P4EST_CHILDREN; c++) {
3870: p4est_qcoord_t quad_coords[3];
3871: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3872: double corner_coords[3];
3873: double vert_coords[3];
3874: PetscInt corner = PetscVertToP4estVert[c];
3876: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) vert_coords[d] = PetscRealPart(values[c * cDim + d]);
3878: quad_coords[0] = quad->x;
3879: quad_coords[1] = quad->y;
3880: #ifdef P4_TO_P8
3881: quad_coords[2] = quad->z;
3882: #endif
3883: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3884: #ifndef P4_TO_P8
3885: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3886: #else
3887: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3888: #endif
3889: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) {
3890: if (fabs(vert_coords[d] - corner_coords[d]) > PETSC_SMALL) {
3891: same_coords = PETSC_FALSE;
3892: break;
3893: }
3894: }
3895: }
3896: if (same_coords) {
3897: PetscCall(PetscSectionSetDof(newSection, cell, 0));
3898: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0));
3899: } else {
3900: PetscCall(PetscSectionSetDof(newSection, cell, cSize));
3901: PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, cSize));
3902: }
3903: PetscCall(DMPlexVecRestoreClosure(plex, oldSection, cVecOld, cell, &cSize, &values));
3904: }
3905: PetscFunctionReturn(PETSC_SUCCESS);
3906: }
3908: static PetscErrorCode PforestLocalizeCell(DM plex, PetscInt cDim, DM_Forest_pforest *pforest, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad, PetscScalar coords[])
3909: {
3910: PetscInt cdof, off;
3912: PetscFunctionBegin;
3913: PetscCall(PetscSectionGetDof(newSection, cell, &cdof));
3914: if (!cdof) PetscFunctionReturn(PETSC_SUCCESS);
3916: PetscCall(PetscSectionGetOffset(newSection, cell, &off));
3917: for (PetscInt c = 0, pos = off; c < P4EST_CHILDREN; c++) {
3918: p4est_qcoord_t quad_coords[3];
3919: p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level);
3920: double corner_coords[3];
3921: PetscInt corner = PetscVertToP4estVert[c];
3923: quad_coords[0] = quad->x;
3924: quad_coords[1] = quad->y;
3925: #ifdef P4_TO_P8
3926: quad_coords[2] = quad->z;
3927: #endif
3928: for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0;
3929: #ifndef P4_TO_P8
3930: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords));
3931: #else
3932: PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords));
3933: #endif
3934: for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) coords[pos++] = corner_coords[d];
3935: for (PetscInt d = PetscMin(cDim, 3); d < cDim; d++) coords[pos++] = 0.;
3936: }
3937: PetscFunctionReturn(PETSC_SUCCESS);
3938: }
3940: static PetscErrorCode DMPforestLocalizeCoordinates(DM dm, DM plex)
3941: {
3942: DM_Forest *forest;
3943: DM_Forest_pforest *pforest;
3944: DM base, cdm, cdmCell;
3945: Vec cVec, cVecOld;
3946: PetscSection oldSection, newSection;
3947: PetscScalar *coords2;
3948: const PetscReal *L;
3949: PetscInt cLocalStart, cLocalEnd, coarsePoint;
3950: PetscInt cDim, newStart, newEnd;
3951: PetscInt v, vStart, vEnd, cp, cStart, cEnd, cEndInterior;
3952: p4est_topidx_t flt, llt, t;
3953: p4est_tree_t *trees;
3954: PetscBool baseLocalized = PETSC_FALSE;
3956: PetscFunctionBegin;
3957: PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
3958: /* we localize on all cells if we don't have a base DM or the base DM coordinates have not been localized */
3959: PetscCall(DMGetCoordinateDim(dm, &cDim));
3960: PetscCall(DMForestGetBaseDM(dm, &base));
3961: if (base) PetscCall(DMGetCoordinatesLocalized(base, &baseLocalized));
3962: if (!baseLocalized) base = NULL;
3963: if (!baseLocalized && !L) PetscFunctionReturn(PETSC_SUCCESS);
3964: PetscCall(DMPlexGetChart(plex, &newStart, &newEnd));
3966: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newSection));
3967: PetscCall(PetscSectionSetNumFields(newSection, 1));
3968: PetscCall(PetscSectionSetFieldComponents(newSection, 0, cDim));
3969: PetscCall(PetscSectionSetChart(newSection, newStart, newEnd));
3971: PetscCall(DMGetCoordinateSection(plex, &oldSection));
3972: PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
3973: PetscCall(DMGetCoordinatesLocal(plex, &cVecOld));
3975: forest = (DM_Forest *)dm->data;
3976: pforest = (DM_Forest_pforest *)forest->data;
3977: cLocalStart = pforest->cLocalStart;
3978: cLocalEnd = pforest->cLocalEnd;
3979: flt = pforest->forest->first_local_tree;
3980: llt = pforest->forest->last_local_tree;
3981: trees = (p4est_tree_t *)pforest->forest->trees->array;
3983: PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd));
3984: PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3985: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
3986: cp = 0;
3987: if (cLocalStart > 0) {
3988: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
3989: PetscInt cell;
3991: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
3992: p4est_quadrant_t *quad = &ghosts[cell];
3994: coarsePoint = quad->p.which_tree;
3995: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
3996: }
3997: }
3998: for (t = flt; t <= llt; t++) {
3999: p4est_tree_t *tree = &trees[t];
4000: PetscInt offset = cLocalStart + tree->quadrants_offset;
4001: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
4002: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4003: PetscInt i;
4005: if (!numQuads) continue;
4006: coarsePoint = t;
4007: for (i = 0; i < numQuads; i++, cp++) {
4008: PetscInt cell = i + offset;
4009: p4est_quadrant_t *quad = &quads[i];
4011: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4012: }
4013: }
4014: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4015: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4016: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4017: PetscInt count;
4019: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4020: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4021: coarsePoint = quad->p.which_tree;
4022: PetscInt cell = count + cLocalEnd;
4024: PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad));
4025: }
4026: }
4027: PetscAssert(cp == cEnd - cStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of fine cells %" PetscInt_FMT " != %" PetscInt_FMT, cp, cEnd - cStart);
4029: PetscCall(PetscSectionSetUp(newSection));
4030: PetscCall(DMGetCoordinateDM(plex, &cdm));
4031: PetscCall(DMClone(cdm, &cdmCell));
4032: PetscCall(DMSetCellCoordinateDM(plex, cdmCell));
4033: PetscCall(DMDestroy(&cdmCell));
4034: PetscCall(DMSetCellCoordinateSection(plex, cDim, newSection));
4035: PetscCall(PetscSectionGetStorageSize(newSection, &v));
4036: PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
4037: PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
4038: PetscCall(VecSetBlockSize(cVec, cDim));
4039: PetscCall(VecSetSizes(cVec, v, PETSC_DETERMINE));
4040: PetscCall(VecSetType(cVec, VECSTANDARD));
4041: PetscCall(VecSet(cVec, PETSC_MIN_REAL));
4043: /* Localize coordinates on cells if needed */
4044: PetscCall(VecGetArray(cVec, &coords2));
4045: cp = 0;
4046: if (cLocalStart > 0) {
4047: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4048: PetscInt cell;
4050: for (cell = 0; cell < cLocalStart; ++cell, cp++) {
4051: p4est_quadrant_t *quad = &ghosts[cell];
4053: coarsePoint = quad->p.which_tree;
4054: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4055: }
4056: }
4057: for (t = flt; t <= llt; t++) {
4058: p4est_tree_t *tree = &trees[t];
4059: PetscInt offset = cLocalStart + tree->quadrants_offset;
4060: PetscInt numQuads = (PetscInt)tree->quadrants.elem_count;
4061: p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array;
4062: PetscInt i;
4064: if (!numQuads) continue;
4065: coarsePoint = t;
4066: for (i = 0; i < numQuads; i++, cp++) {
4067: PetscInt cell = i + offset;
4068: p4est_quadrant_t *quad = &quads[i];
4070: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4071: }
4072: }
4073: if (cLocalEnd - cLocalStart < cEnd - cStart) {
4074: p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array;
4075: PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count;
4076: PetscInt count;
4078: for (count = 0; count < numGhosts - cLocalStart; count++, cp++) {
4079: p4est_quadrant_t *quad = &ghosts[count + cLocalStart];
4080: coarsePoint = quad->p.which_tree;
4081: PetscInt cell = count + cLocalEnd;
4083: PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2));
4084: }
4085: }
4086: PetscCall(VecRestoreArray(cVec, &coords2));
4087: PetscCall(DMSetCellCoordinatesLocal(plex, cVec));
4088: PetscCall(VecDestroy(&cVec));
4089: PetscCall(PetscSectionDestroy(&newSection));
4090: PetscFunctionReturn(PETSC_SUCCESS);
4091: }
4093: #define DMForestClearAdaptivityForest_pforest _append_pforest(DMForestClearAdaptivityForest)
4094: static PetscErrorCode DMForestClearAdaptivityForest_pforest(DM dm)
4095: {
4096: DM_Forest *forest;
4097: DM_Forest_pforest *pforest;
4099: PetscFunctionBegin;
4100: forest = (DM_Forest *)dm->data;
4101: pforest = (DM_Forest_pforest *)forest->data;
4102: PetscCall(PetscSFDestroy(&pforest->pointAdaptToSelfSF));
4103: PetscCall(PetscSFDestroy(&pforest->pointSelfToAdaptSF));
4104: PetscCall(PetscFree(pforest->pointAdaptToSelfCids));
4105: PetscCall(PetscFree(pforest->pointSelfToAdaptCids));
4106: PetscFunctionReturn(PETSC_SUCCESS);
4107: }
4109: static PetscErrorCode DMConvert_pforest_plex(DM dm, DMType newtype, DM *plex)
4110: {
4111: DM_Forest *forest;
4112: DM_Forest_pforest *pforest;
4113: DM refTree, newPlex, base;
4114: PetscInt adjDim, adjCodim, coordDim;
4115: MPI_Comm comm;
4116: PetscBool isPforest;
4117: PetscInt dim;
4118: PetscInt overlap;
4119: p4est_connect_type_t ctype;
4120: p4est_locidx_t first_local_quad = -1;
4121: sc_array_t *points_per_dim, *cone_sizes, *cones, *cone_orientations, *coords, *children, *parents, *childids, *leaves, *remotes;
4122: PetscSection parentSection;
4123: PetscSF pointSF;
4124: size_t zz, count;
4125: PetscInt pStart, pEnd;
4126: DMLabel ghostLabelBase = NULL;
4128: PetscFunctionBegin;
4130: comm = PetscObjectComm((PetscObject)dm);
4131: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPFOREST, &isPforest));
4132: PetscCheck(isPforest, comm, PETSC_ERR_ARG_WRONG, "Expected DM type %s, got %s", DMPFOREST, ((PetscObject)dm)->type_name);
4133: PetscCall(DMSetUp(dm));
4134: PetscCall(DMGetDimension(dm, &dim));
4135: PetscCheck(dim == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Expected DM dimension %d, got %" PetscInt_FMT, P4EST_DIM, dim);
4136: forest = (DM_Forest *)dm->data;
4137: pforest = (DM_Forest_pforest *)forest->data;
4138: PetscCall(DMForestGetBaseDM(dm, &base));
4139: if (base) PetscCall(DMGetLabel(base, "ghost", &ghostLabelBase));
4140: if (!pforest->plex) {
4141: PetscMPIInt size;
4142: const char *name;
4144: PetscCallMPI(MPI_Comm_size(comm, &size));
4145: PetscCall(DMCreate(comm, &newPlex));
4146: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
4147: PetscCall(PetscObjectSetName((PetscObject)newPlex, name));
4148: PetscCall(DMSetType(newPlex, DMPLEX));
4149: PetscCall(DMSetMatType(newPlex, dm->mattype));
4150: /* share labels */
4151: PetscCall(DMCopyLabels(dm, newPlex, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4152: PetscCall(DMForestGetAdjacencyDimension(dm, &adjDim));
4153: PetscCall(DMForestGetAdjacencyCodimension(dm, &adjCodim));
4154: PetscCall(DMGetCoordinateDim(dm, &coordDim));
4155: if (adjDim == 0) {
4156: ctype = P4EST_CONNECT_FULL;
4157: } else if (adjCodim == 1) {
4158: ctype = P4EST_CONNECT_FACE;
4159: #if defined(P4_TO_P8)
4160: } else if (adjDim == 1) {
4161: ctype = P8EST_CONNECT_EDGE;
4162: #endif
4163: } else {
4164: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid adjacency dimension %" PetscInt_FMT, adjDim);
4165: }
4166: PetscCheck(ctype == P4EST_CONNECT_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Adjacency dimension %" PetscInt_FMT " / codimension %" PetscInt_FMT " not supported yet", adjDim, adjCodim);
4167: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
4168: PetscCall(DMPlexSetOverlap_Plex(newPlex, NULL, overlap));
4170: points_per_dim = sc_array_new(sizeof(p4est_locidx_t));
4171: cone_sizes = sc_array_new(sizeof(p4est_locidx_t));
4172: cones = sc_array_new(sizeof(p4est_locidx_t));
4173: cone_orientations = sc_array_new(sizeof(p4est_locidx_t));
4174: coords = sc_array_new(3 * sizeof(double));
4175: children = sc_array_new(sizeof(p4est_locidx_t));
4176: parents = sc_array_new(sizeof(p4est_locidx_t));
4177: childids = sc_array_new(sizeof(p4est_locidx_t));
4178: leaves = sc_array_new(sizeof(p4est_locidx_t));
4179: remotes = sc_array_new(2 * sizeof(p4est_locidx_t));
4181: PetscCallP4est(p4est_get_plex_data_ext, (pforest->forest, &pforest->ghost, &pforest->lnodes, ctype, (int)((size > 1) ? overlap : 0), &first_local_quad, points_per_dim, cone_sizes, cones, cone_orientations, coords, children, parents, childids, leaves, remotes, 1));
4183: pforest->cLocalStart = (PetscInt)first_local_quad;
4184: pforest->cLocalEnd = pforest->cLocalStart + (PetscInt)pforest->forest->local_num_quadrants;
4185: PetscCall(locidx_to_PetscInt(points_per_dim));
4186: PetscCall(locidx_to_PetscInt(cone_sizes));
4187: PetscCall(locidx_to_PetscInt(cones));
4188: PetscCall(locidx_to_PetscInt(cone_orientations));
4189: PetscCall(coords_double_to_PetscScalar(coords, coordDim));
4190: PetscCall(locidx_to_PetscInt(children));
4191: PetscCall(locidx_to_PetscInt(parents));
4192: PetscCall(locidx_to_PetscInt(childids));
4193: PetscCall(locidx_to_PetscInt(leaves));
4194: PetscCall(locidx_pair_to_PetscSFNode(remotes));
4196: PetscCall(DMSetDimension(newPlex, P4EST_DIM));
4197: PetscCall(DMSetCoordinateDim(newPlex, coordDim));
4198: PetscCall(DMPlexSetMaxProjectionHeight(newPlex, P4EST_DIM - 1));
4199: PetscCall(DMPlexCreateFromDAG(newPlex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array));
4200: PetscCall(DMPlexConvertOldOrientations_Internal(newPlex));
4201: PetscCall(DMCreateReferenceTree_pforest(comm, &refTree));
4202: PetscCall(DMPlexSetReferenceTree(newPlex, refTree));
4203: PetscCall(PetscSectionCreate(comm, &parentSection));
4204: PetscCall(DMPlexGetChart(newPlex, &pStart, &pEnd));
4205: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
4206: count = children->elem_count;
4207: for (zz = 0; zz < count; zz++) {
4208: PetscInt child = *((PetscInt *)sc_array_index(children, zz));
4210: PetscCall(PetscSectionSetDof(parentSection, child, 1));
4211: }
4212: PetscCall(PetscSectionSetUp(parentSection));
4213: PetscCall(DMPlexSetTree(newPlex, parentSection, (PetscInt *)parents->array, (PetscInt *)childids->array));
4214: PetscCall(PetscSectionDestroy(&parentSection));
4215: PetscCall(PetscSFCreate(comm, &pointSF));
4216: /*
4217: These arrays defining the sf are from the p4est library, but the code there shows the leaves being populated in increasing order.
4218: https://gitlab.com/petsc/petsc/merge_requests/2248#note_240186391
4219: */
4220: PetscCall(PetscSFSetGraph(pointSF, pEnd - pStart, (PetscInt)leaves->elem_count, (PetscInt *)leaves->array, PETSC_COPY_VALUES, (PetscSFNode *)remotes->array, PETSC_COPY_VALUES));
4221: PetscCall(DMSetPointSF(newPlex, pointSF));
4222: PetscCall(DMSetPointSF(dm, pointSF));
4223: {
4224: DM coordDM;
4226: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4227: PetscCall(DMSetPointSF(coordDM, pointSF));
4228: }
4229: PetscCall(PetscSFDestroy(&pointSF));
4230: sc_array_destroy(points_per_dim);
4231: sc_array_destroy(cone_sizes);
4232: sc_array_destroy(cones);
4233: sc_array_destroy(cone_orientations);
4234: sc_array_destroy(coords);
4235: sc_array_destroy(children);
4236: sc_array_destroy(parents);
4237: sc_array_destroy(childids);
4238: sc_array_destroy(leaves);
4239: sc_array_destroy(remotes);
4241: {
4242: const PetscReal *maxCell, *Lstart, *L;
4244: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
4245: PetscCall(DMSetPeriodicity(newPlex, maxCell, Lstart, L));
4246: PetscCall(DMPforestLocalizeCoordinates(dm, newPlex));
4247: }
4249: if (overlap > 0) { /* the p4est routine can't set all of the coordinates in its routine if there is overlap */
4250: Vec coordsGlobal, coordsLocal;
4251: const PetscScalar *globalArray;
4252: PetscScalar *localArray;
4253: PetscSF coordSF;
4254: DM coordDM;
4256: PetscCall(DMGetCoordinateDM(newPlex, &coordDM));
4257: PetscCall(DMGetSectionSF(coordDM, &coordSF));
4258: PetscCall(DMGetCoordinates(newPlex, &coordsGlobal));
4259: PetscCall(DMGetCoordinatesLocal(newPlex, &coordsLocal));
4260: PetscCall(VecGetArrayRead(coordsGlobal, &globalArray));
4261: PetscCall(VecGetArray(coordsLocal, &localArray));
4262: PetscCall(PetscSFBcastBegin(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4263: PetscCall(PetscSFBcastEnd(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE));
4264: PetscCall(VecRestoreArray(coordsLocal, &localArray));
4265: PetscCall(VecRestoreArrayRead(coordsGlobal, &globalArray));
4266: PetscCall(DMSetCoordinatesLocal(newPlex, coordsLocal));
4267: }
4268: PetscCall(DMPforestMapCoordinates(dm, newPlex));
4270: pforest->plex = newPlex;
4272: /* copy labels */
4273: PetscCall(DMPforestLabelsFinalize(dm, newPlex));
4275: if (ghostLabelBase || pforest->ghostName) { /* we have to do this after copying labels because the labels drive the construction of ghost cells */
4276: PetscInt numAdded;
4277: DM newPlexGhosted;
4278: void *ctx;
4280: PetscCall(DMPlexConstructGhostCells(newPlex, pforest->ghostName, &numAdded, &newPlexGhosted));
4281: PetscCall(DMGetApplicationContext(newPlex, &ctx));
4282: PetscCall(DMSetApplicationContext(newPlexGhosted, ctx));
4283: /* we want the sf for the ghost dm to be the one for the p4est dm as well */
4284: PetscCall(DMGetPointSF(newPlexGhosted, &pointSF));
4285: PetscCall(DMSetPointSF(dm, pointSF));
4286: PetscCall(DMDestroy(&newPlex));
4287: PetscCall(DMPlexSetReferenceTree(newPlexGhosted, refTree));
4288: PetscCall(DMForestClearAdaptivityForest_pforest(dm));
4289: newPlex = newPlexGhosted;
4291: /* share the labels back */
4292: PetscCall(DMDestroyLabelLinkList_Internal(dm));
4293: PetscCall(DMCopyLabels(newPlex, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL));
4294: pforest->plex = newPlex;
4295: }
4296: PetscCall(DMDestroy(&refTree));
4297: if (dm->setfromoptionscalled) {
4298: PetscObjectOptionsBegin((PetscObject)newPlex);
4299: PetscCall(DMSetFromOptions_NonRefinement_Plex(newPlex, PetscOptionsObject));
4300: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)newPlex, PetscOptionsObject));
4301: PetscOptionsEnd();
4302: }
4303: PetscCall(DMViewFromOptions(newPlex, NULL, "-dm_p4est_plex_view"));
4304: {
4305: DM cdm;
4306: PetscSection coordsSec;
4307: Vec coords;
4308: PetscInt cDim;
4310: PetscCall(DMGetCoordinateDim(newPlex, &cDim));
4311: PetscCall(DMGetCoordinateSection(newPlex, &coordsSec));
4312: PetscCall(DMSetCoordinateSection(dm, cDim, coordsSec));
4313: PetscCall(DMGetCoordinatesLocal(newPlex, &coords));
4314: PetscCall(DMSetCoordinatesLocal(dm, coords));
4315: PetscCall(DMGetCoordinateDM(newPlex, &cdm));
4316: if (cdm) {
4317: PetscFE fe;
4318: #if !defined(P4_TO_P8)
4319: DMPolytopeType celltype = DM_POLYTOPE_QUADRILATERAL;
4320: #else
4321: DMPolytopeType celltype = DM_POLYTOPE_HEXAHEDRON;
4322: #endif
4324: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dim, celltype, 1, PETSC_DEFAULT, &fe));
4325: PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe));
4326: PetscCall(PetscFEDestroy(&fe));
4327: PetscCall(DMCreateDS(cdm));
4328: }
4329: PetscCall(DMGetCellCoordinateDM(newPlex, &cdm));
4330: if (cdm) PetscCall(DMSetCellCoordinateDM(dm, cdm));
4331: PetscCall(DMGetCellCoordinateSection(newPlex, &coordsSec));
4332: if (coordsSec) PetscCall(DMSetCellCoordinateSection(dm, cDim, coordsSec));
4333: PetscCall(DMGetCellCoordinatesLocal(newPlex, &coords));
4334: if (coords) PetscCall(DMSetCellCoordinatesLocal(dm, coords));
4335: }
4336: } else {
4337: PetscCall(DMCopyLabels(dm, pforest->plex, PETSC_OWN_POINTER, PETSC_FALSE, DM_COPY_LABELS_REPLACE));
4338: }
4339: newPlex = pforest->plex;
4340: if (plex) {
4341: PetscCall(DMClone(newPlex, plex));
4342: #if 0
4343: PetscCall(DMGetCoordinateDM(newPlex,&coordDM));
4344: PetscCall(DMSetCoordinateDM(*plex,coordDM));
4345: PetscCall(DMGetCellCoordinateDM(newPlex,&coordDM));
4346: PetscCall(DMSetCellCoordinateDM(*plex,coordDM));
4347: #endif
4348: PetscCall(DMShareDiscretization(dm, *plex));
4349: }
4350: PetscFunctionReturn(PETSC_SUCCESS);
4351: }
4353: static PetscErrorCode DMSetFromOptions_pforest(DM dm, PetscOptionItems *PetscOptionsObject)
4354: {
4355: DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4356: char stringBuffer[256];
4357: PetscBool flg;
4359: PetscFunctionBegin;
4360: PetscCall(DMSetFromOptions_Forest(dm, PetscOptionsObject));
4361: PetscOptionsHeadBegin(PetscOptionsObject, "DM" P4EST_STRING " options");
4362: PetscCall(PetscOptionsBool("-dm_p4est_partition_for_coarsening", "partition forest to allow for coarsening", "DMP4estSetPartitionForCoarsening", pforest->partition_for_coarsening, &pforest->partition_for_coarsening, NULL));
4363: PetscCall(PetscOptionsString("-dm_p4est_ghost_label_name", "the name of the ghost label when converting from a DMPlex", NULL, NULL, stringBuffer, sizeof(stringBuffer), &flg));
4364: PetscOptionsHeadEnd();
4365: if (flg) {
4366: PetscCall(PetscFree(pforest->ghostName));
4367: PetscCall(PetscStrallocpy(stringBuffer, &pforest->ghostName));
4368: }
4369: PetscFunctionReturn(PETSC_SUCCESS);
4370: }
4372: #if !defined(P4_TO_P8)
4373: #define DMPforestGetPartitionForCoarsening DMP4estGetPartitionForCoarsening
4374: #define DMPforestSetPartitionForCoarsening DMP4estSetPartitionForCoarsening
4375: #else
4376: #define DMPforestGetPartitionForCoarsening DMP8estGetPartitionForCoarsening
4377: #define DMPforestSetPartitionForCoarsening DMP8estSetPartitionForCoarsening
4378: #endif
4380: PETSC_EXTERN PetscErrorCode DMPforestGetPartitionForCoarsening(DM dm, PetscBool *flg)
4381: {
4382: DM_Forest_pforest *pforest;
4384: PetscFunctionBegin;
4386: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4387: *flg = pforest->partition_for_coarsening;
4388: PetscFunctionReturn(PETSC_SUCCESS);
4389: }
4391: PETSC_EXTERN PetscErrorCode DMPforestSetPartitionForCoarsening(DM dm, PetscBool flg)
4392: {
4393: DM_Forest_pforest *pforest;
4395: PetscFunctionBegin;
4397: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4398: pforest->partition_for_coarsening = flg;
4399: PetscFunctionReturn(PETSC_SUCCESS);
4400: }
4402: static PetscErrorCode DMPforestGetPlex(DM dm, DM *plex)
4403: {
4404: DM_Forest_pforest *pforest;
4406: PetscFunctionBegin;
4407: if (plex) *plex = NULL;
4408: PetscCall(DMSetUp(dm));
4409: pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data;
4410: if (!pforest->plex) PetscCall(DMConvert_pforest_plex(dm, DMPLEX, NULL));
4411: PetscCall(DMShareDiscretization(dm, pforest->plex));
4412: if (plex) *plex = pforest->plex;
4413: PetscFunctionReturn(PETSC_SUCCESS);
4414: }
4416: #define DMCreateInterpolation_pforest _append_pforest(DMCreateInterpolation)
4417: static PetscErrorCode DMCreateInterpolation_pforest(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
4418: {
4419: PetscSection gsc, gsf;
4420: PetscInt m, n;
4421: DM cdm;
4423: PetscFunctionBegin;
4424: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4425: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
4426: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4427: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &n));
4429: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), interpolation));
4430: PetscCall(MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4431: PetscCall(MatSetType(*interpolation, MATAIJ));
4433: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4434: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only interpolation from coarse DM for now");
4436: {
4437: DM plexF, plexC;
4438: PetscSF sf;
4439: PetscInt *cids;
4440: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4442: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4443: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4444: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4445: PetscCall(PetscSFSetUp(sf));
4446: PetscCall(DMPlexComputeInterpolatorTree(plexC, plexF, sf, cids, *interpolation));
4447: PetscCall(PetscSFDestroy(&sf));
4448: PetscCall(PetscFree(cids));
4449: }
4450: PetscCall(MatViewFromOptions(*interpolation, NULL, "-interp_mat_view"));
4451: /* Use naive scaling */
4452: PetscCall(DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling));
4453: PetscFunctionReturn(PETSC_SUCCESS);
4454: }
4456: #define DMCreateInjection_pforest _append_pforest(DMCreateInjection)
4457: static PetscErrorCode DMCreateInjection_pforest(DM dmCoarse, DM dmFine, Mat *injection)
4458: {
4459: PetscSection gsc, gsf;
4460: PetscInt m, n;
4461: DM cdm;
4463: PetscFunctionBegin;
4464: PetscCall(DMGetGlobalSection(dmFine, &gsf));
4465: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &n));
4466: PetscCall(DMGetGlobalSection(dmCoarse, &gsc));
4467: PetscCall(PetscSectionGetConstrainedStorageSize(gsc, &m));
4469: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmFine), injection));
4470: PetscCall(MatSetSizes(*injection, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
4471: PetscCall(MatSetType(*injection, MATAIJ));
4473: PetscCall(DMGetCoarseDM(dmFine, &cdm));
4474: PetscCheck(cdm == dmCoarse, PetscObjectComm((PetscObject)dmFine), PETSC_ERR_SUP, "Only injection to coarse DM for now");
4476: {
4477: DM plexF, plexC;
4478: PetscSF sf;
4479: PetscInt *cids;
4480: PetscInt dofPerDim[4] = {1, 1, 1, 1};
4482: PetscCall(DMPforestGetPlex(dmCoarse, &plexC));
4483: PetscCall(DMPforestGetPlex(dmFine, &plexF));
4484: PetscCall(DMPforestGetTransferSF_Internal(dmCoarse, dmFine, dofPerDim, &sf, PETSC_TRUE, &cids));
4485: PetscCall(PetscSFSetUp(sf));
4486: PetscCall(DMPlexComputeInjectorTree(plexC, plexF, sf, cids, *injection));
4487: PetscCall(PetscSFDestroy(&sf));
4488: PetscCall(PetscFree(cids));
4489: }
4490: PetscCall(MatViewFromOptions(*injection, NULL, "-inject_mat_view"));
4491: /* Use naive scaling */
4492: PetscFunctionReturn(PETSC_SUCCESS);
4493: }
4495: #define DMForestTransferVecFromBase_pforest _append_pforest(DMForestTransferVecFromBase)
4496: static PetscErrorCode DMForestTransferVecFromBase_pforest(DM dm, Vec vecIn, Vec vecOut)
4497: {
4498: DM dmIn, dmVecIn, base, basec, plex, coarseDM;
4499: DM *hierarchy;
4500: PetscSF sfRed = NULL;
4501: PetscDS ds;
4502: Vec vecInLocal, vecOutLocal;
4503: DMLabel subpointMap;
4504: PetscInt minLevel, mh, n_hi, i;
4505: PetscBool hiforest, *hierarchy_forest;
4507: PetscFunctionBegin;
4508: PetscCall(VecGetDM(vecIn, &dmVecIn));
4509: PetscCall(DMGetDS(dmVecIn, &ds));
4510: PetscCheck(ds, PetscObjectComm((PetscObject)dmVecIn), PETSC_ERR_SUP, "Cannot transfer without a PetscDS object");
4511: { /* we cannot stick user contexts into function callbacks for DMProjectFieldLocal! */
4512: PetscSection section;
4513: PetscInt Nf;
4515: PetscCall(DMGetLocalSection(dmVecIn, §ion));
4516: PetscCall(PetscSectionGetNumFields(section, &Nf));
4517: PetscCheck(Nf <= 3, PetscObjectComm((PetscObject)dmVecIn), PETSC_ERR_SUP, "Number of fields %" PetscInt_FMT " are currently not supported! Send an email at petsc-dev@mcs.anl.gov", Nf);
4518: }
4519: PetscCall(DMForestGetMinimumRefinement(dm, &minLevel));
4520: PetscCheck(!minLevel, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot transfer with minimum refinement set to %" PetscInt_FMT ". Rerun with DMForestSetMinimumRefinement(dm,0)", minLevel);
4521: PetscCall(DMForestGetBaseDM(dm, &base));
4522: PetscCheck(base, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing base DM");
4524: PetscCall(VecSet(vecOut, 0.0));
4525: if (dmVecIn == base) { /* sequential runs */
4526: PetscCall(PetscObjectReference((PetscObject)vecIn));
4527: } else {
4528: PetscSection secIn, secInRed;
4529: Vec vecInRed, vecInLocal;
4531: PetscCall(PetscObjectQuery((PetscObject)base, "_base_migration_sf", (PetscObject *)&sfRed));
4532: PetscCheck(sfRed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not the DM set with DMForestSetBaseDM()");
4533: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmVecIn), &secInRed));
4534: PetscCall(VecCreate(PETSC_COMM_SELF, &vecInRed));
4535: PetscCall(DMGetLocalSection(dmVecIn, &secIn));
4536: PetscCall(DMGetLocalVector(dmVecIn, &vecInLocal));
4537: PetscCall(DMGlobalToLocalBegin(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4538: PetscCall(DMGlobalToLocalEnd(dmVecIn, vecIn, INSERT_VALUES, vecInLocal));
4539: PetscCall(DMPlexDistributeField(dmVecIn, sfRed, secIn, vecInLocal, secInRed, vecInRed));
4540: PetscCall(DMRestoreLocalVector(dmVecIn, &vecInLocal));
4541: PetscCall(PetscSectionDestroy(&secInRed));
4542: vecIn = vecInRed;
4543: }
4545: /* we first search through the AdaptivityForest hierarchy
4546: once we found the first disconnected forest, we upsweep the DM hierarchy */
4547: hiforest = PETSC_TRUE;
4549: /* upsweep to the coarsest DM */
4550: n_hi = 0;
4551: coarseDM = dm;
4552: do {
4553: PetscBool isforest;
4555: dmIn = coarseDM;
4556: /* need to call DMSetUp to have the hierarchy recursively setup */
4557: PetscCall(DMSetUp(dmIn));
4558: PetscCall(DMIsForest(dmIn, &isforest));
4559: PetscCheck(isforest, PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Cannot currently transfer through a mixed hierarchy! Found DM type %s", ((PetscObject)dmIn)->type_name);
4560: coarseDM = NULL;
4561: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4562: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4563: hiforest = PETSC_FALSE;
4564: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4565: }
4566: n_hi++;
4567: } while (coarseDM);
4569: PetscCall(PetscMalloc2(n_hi, &hierarchy, n_hi, &hierarchy_forest));
4571: i = 0;
4572: hiforest = PETSC_TRUE;
4573: coarseDM = dm;
4574: do {
4575: dmIn = coarseDM;
4576: coarseDM = NULL;
4577: if (hiforest) PetscCall(DMForestGetAdaptivityForest(dmIn, &coarseDM));
4578: if (!coarseDM) { /* DMForest hierarchy ended, we keep upsweeping through the DM hierarchy */
4579: hiforest = PETSC_FALSE;
4580: PetscCall(DMGetCoarseDM(dmIn, &coarseDM));
4581: }
4582: i++;
4583: hierarchy[n_hi - i] = dmIn;
4584: } while (coarseDM);
4586: /* project base vector on the coarsest forest (minimum refinement = 0) */
4587: PetscCall(DMPforestGetPlex(dmIn, &plex));
4589: /* Check this plex is compatible with the base */
4590: {
4591: IS gnum[2];
4592: PetscInt ncells[2], gncells[2];
4594: PetscCall(DMPlexGetCellNumbering(base, &gnum[0]));
4595: PetscCall(DMPlexGetCellNumbering(plex, &gnum[1]));
4596: PetscCall(ISGetMinMax(gnum[0], NULL, &ncells[0]));
4597: PetscCall(ISGetMinMax(gnum[1], NULL, &ncells[1]));
4598: PetscCall(MPIU_Allreduce(ncells, gncells, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
4599: PetscCheck(gncells[0] == gncells[1], PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Invalid number of base cells! Expected %" PetscInt_FMT ", found %" PetscInt_FMT, gncells[0] + 1, gncells[1] + 1);
4600: }
4602: PetscCall(DMGetLabel(dmIn, "_forest_base_subpoint_map", &subpointMap));
4603: PetscCheck(subpointMap, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing _forest_base_subpoint_map label");
4605: PetscCall(DMPlexGetMaxProjectionHeight(base, &mh));
4606: PetscCall(DMPlexSetMaxProjectionHeight(plex, mh));
4608: PetscCall(DMClone(base, &basec));
4609: PetscCall(DMCopyDisc(dmVecIn, basec));
4610: if (sfRed) {
4611: PetscCall(PetscObjectReference((PetscObject)vecIn));
4612: vecInLocal = vecIn;
4613: } else {
4614: PetscCall(DMCreateLocalVector(basec, &vecInLocal));
4615: PetscCall(DMGlobalToLocalBegin(basec, vecIn, INSERT_VALUES, vecInLocal));
4616: PetscCall(DMGlobalToLocalEnd(basec, vecIn, INSERT_VALUES, vecInLocal));
4617: }
4619: PetscCall(DMGetLocalVector(dmIn, &vecOutLocal));
4620: { /* get degrees of freedom ordered onto dmIn */
4621: PetscSF basetocoarse;
4622: PetscInt bStart, bEnd, nroots;
4623: PetscInt iStart, iEnd, nleaves, leaf;
4624: PetscMPIInt rank;
4625: PetscSFNode *remotes;
4626: PetscSection secIn, secOut;
4627: PetscInt *remoteOffsets;
4628: PetscSF transferSF;
4629: const PetscScalar *inArray;
4630: PetscScalar *outArray;
4632: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)basec), &rank));
4633: PetscCall(DMPlexGetChart(basec, &bStart, &bEnd));
4634: nroots = PetscMax(bEnd - bStart, 0);
4635: PetscCall(DMPlexGetChart(plex, &iStart, &iEnd));
4636: nleaves = PetscMax(iEnd - iStart, 0);
4638: PetscCall(PetscMalloc1(nleaves, &remotes));
4639: for (leaf = iStart; leaf < iEnd; leaf++) {
4640: PetscInt index;
4642: remotes[leaf - iStart].rank = rank;
4643: PetscCall(DMLabelGetValue(subpointMap, leaf, &index));
4644: remotes[leaf - iStart].index = index;
4645: }
4647: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)basec), &basetocoarse));
4648: PetscCall(PetscSFSetGraph(basetocoarse, nroots, nleaves, NULL, PETSC_OWN_POINTER, remotes, PETSC_OWN_POINTER));
4649: PetscCall(PetscSFSetUp(basetocoarse));
4650: PetscCall(DMGetLocalSection(basec, &secIn));
4651: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmIn), &secOut));
4652: PetscCall(PetscSFDistributeSection(basetocoarse, secIn, &remoteOffsets, secOut));
4653: PetscCall(PetscSFCreateSectionSF(basetocoarse, secIn, remoteOffsets, secOut, &transferSF));
4654: PetscCall(PetscFree(remoteOffsets));
4655: PetscCall(VecGetArrayWrite(vecOutLocal, &outArray));
4656: PetscCall(VecGetArrayRead(vecInLocal, &inArray));
4657: PetscCall(PetscSFBcastBegin(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4658: PetscCall(PetscSFBcastEnd(transferSF, MPIU_SCALAR, inArray, outArray, MPI_REPLACE));
4659: PetscCall(VecRestoreArrayRead(vecInLocal, &inArray));
4660: PetscCall(VecRestoreArrayWrite(vecOutLocal, &outArray));
4661: PetscCall(PetscSFDestroy(&transferSF));
4662: PetscCall(PetscSectionDestroy(&secOut));
4663: PetscCall(PetscSFDestroy(&basetocoarse));
4664: }
4665: PetscCall(VecDestroy(&vecInLocal));
4666: PetscCall(DMDestroy(&basec));
4667: PetscCall(VecDestroy(&vecIn));
4669: /* output */
4670: if (n_hi > 1) { /* downsweep the stored hierarchy */
4671: Vec vecOut1, vecOut2;
4672: DM fineDM;
4674: PetscCall(DMGetGlobalVector(dmIn, &vecOut1));
4675: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut1));
4676: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4677: for (i = 1; i < n_hi - 1; i++) {
4678: fineDM = hierarchy[i];
4679: PetscCall(DMGetGlobalVector(fineDM, &vecOut2));
4680: PetscCall(DMForestTransferVec(dmIn, vecOut1, fineDM, vecOut2, PETSC_TRUE, 0.0));
4681: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4682: vecOut1 = vecOut2;
4683: dmIn = fineDM;
4684: }
4685: PetscCall(DMForestTransferVec(dmIn, vecOut1, dm, vecOut, PETSC_TRUE, 0.0));
4686: PetscCall(DMRestoreGlobalVector(dmIn, &vecOut1));
4687: } else {
4688: PetscCall(DMLocalToGlobal(dmIn, vecOutLocal, INSERT_VALUES, vecOut));
4689: PetscCall(DMRestoreLocalVector(dmIn, &vecOutLocal));
4690: }
4691: PetscCall(PetscFree2(hierarchy, hierarchy_forest));
4692: PetscFunctionReturn(PETSC_SUCCESS);
4693: }
4695: #define DMForestTransferVec_pforest _append_pforest(DMForestTransferVec)
4696: static PetscErrorCode DMForestTransferVec_pforest(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscBool useBCs, PetscReal time)
4697: {
4698: DM adaptIn, adaptOut, plexIn, plexOut;
4699: DM_Forest *forestIn, *forestOut, *forestAdaptIn, *forestAdaptOut;
4700: PetscInt dofPerDim[] = {1, 1, 1, 1};
4701: PetscSF inSF = NULL, outSF = NULL;
4702: PetscInt *inCids = NULL, *outCids = NULL;
4703: DMAdaptFlag purposeIn, purposeOut;
4705: PetscFunctionBegin;
4706: forestOut = (DM_Forest *)dmOut->data;
4707: forestIn = (DM_Forest *)dmIn->data;
4709: PetscCall(DMForestGetAdaptivityForest(dmOut, &adaptOut));
4710: PetscCall(DMForestGetAdaptivityPurpose(dmOut, &purposeOut));
4711: forestAdaptOut = adaptOut ? (DM_Forest *)adaptOut->data : NULL;
4713: PetscCall(DMForestGetAdaptivityForest(dmIn, &adaptIn));
4714: PetscCall(DMForestGetAdaptivityPurpose(dmIn, &purposeIn));
4715: forestAdaptIn = adaptIn ? (DM_Forest *)adaptIn->data : NULL;
4717: if (forestAdaptOut == forestIn) {
4718: switch (purposeOut) {
4719: case DM_ADAPT_REFINE:
4720: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4721: PetscCall(PetscSFSetUp(inSF));
4722: break;
4723: case DM_ADAPT_COARSEN:
4724: case DM_ADAPT_COARSEN_LAST:
4725: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &outCids));
4726: PetscCall(PetscSFSetUp(outSF));
4727: break;
4728: default:
4729: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4730: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4731: PetscCall(PetscSFSetUp(inSF));
4732: PetscCall(PetscSFSetUp(outSF));
4733: }
4734: } else if (forestAdaptIn == forestOut) {
4735: switch (purposeIn) {
4736: case DM_ADAPT_REFINE:
4737: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_TRUE, &inCids));
4738: PetscCall(PetscSFSetUp(outSF));
4739: break;
4740: case DM_ADAPT_COARSEN:
4741: case DM_ADAPT_COARSEN_LAST:
4742: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4743: PetscCall(PetscSFSetUp(inSF));
4744: break;
4745: default:
4746: PetscCall(DMPforestGetTransferSF_Internal(dmIn, dmOut, dofPerDim, &inSF, PETSC_TRUE, &inCids));
4747: PetscCall(DMPforestGetTransferSF_Internal(dmOut, dmIn, dofPerDim, &outSF, PETSC_FALSE, &outCids));
4748: PetscCall(PetscSFSetUp(inSF));
4749: PetscCall(PetscSFSetUp(outSF));
4750: }
4751: } else SETERRQ(PetscObjectComm((PetscObject)dmIn), PETSC_ERR_SUP, "Only support transfer from pre-adaptivity to post-adaptivity right now");
4752: PetscCall(DMPforestGetPlex(dmIn, &plexIn));
4753: PetscCall(DMPforestGetPlex(dmOut, &plexOut));
4755: PetscCall(DMPlexTransferVecTree(plexIn, vecIn, plexOut, vecOut, inSF, outSF, inCids, outCids, useBCs, time));
4756: PetscCall(PetscFree(inCids));
4757: PetscCall(PetscFree(outCids));
4758: PetscCall(PetscSFDestroy(&inSF));
4759: PetscCall(PetscSFDestroy(&outSF));
4760: PetscCall(PetscFree(inCids));
4761: PetscCall(PetscFree(outCids));
4762: PetscFunctionReturn(PETSC_SUCCESS);
4763: }
4765: #define DMCreateCoordinateDM_pforest _append_pforest(DMCreateCoordinateDM)
4766: static PetscErrorCode DMCreateCoordinateDM_pforest(DM dm, DM *cdm)
4767: {
4768: DM plex;
4770: PetscFunctionBegin;
4772: PetscCall(DMPforestGetPlex(dm, &plex));
4773: PetscCall(DMGetCoordinateDM(plex, cdm));
4774: PetscCall(PetscObjectReference((PetscObject)*cdm));
4775: PetscFunctionReturn(PETSC_SUCCESS);
4776: }
4778: #define VecViewLocal_pforest _append_pforest(VecViewLocal)
4779: static PetscErrorCode VecViewLocal_pforest(Vec vec, PetscViewer viewer)
4780: {
4781: DM dm, plex;
4783: PetscFunctionBegin;
4784: PetscCall(VecGetDM(vec, &dm));
4785: PetscCall(PetscObjectReference((PetscObject)dm));
4786: PetscCall(DMPforestGetPlex(dm, &plex));
4787: PetscCall(VecSetDM(vec, plex));
4788: PetscCall(VecView_Plex_Local(vec, viewer));
4789: PetscCall(VecSetDM(vec, dm));
4790: PetscCall(DMDestroy(&dm));
4791: PetscFunctionReturn(PETSC_SUCCESS);
4792: }
4794: #define VecView_pforest _append_pforest(VecView)
4795: static PetscErrorCode VecView_pforest(Vec vec, PetscViewer viewer)
4796: {
4797: DM dm, plex;
4799: PetscFunctionBegin;
4800: PetscCall(VecGetDM(vec, &dm));
4801: PetscCall(PetscObjectReference((PetscObject)dm));
4802: PetscCall(DMPforestGetPlex(dm, &plex));
4803: PetscCall(VecSetDM(vec, plex));
4804: PetscCall(VecView_Plex(vec, viewer));
4805: PetscCall(VecSetDM(vec, dm));
4806: PetscCall(DMDestroy(&dm));
4807: PetscFunctionReturn(PETSC_SUCCESS);
4808: }
4810: #define VecView_pforest_Native _infix_pforest(VecView, _Native)
4811: static PetscErrorCode VecView_pforest_Native(Vec vec, PetscViewer viewer)
4812: {
4813: DM dm, plex;
4815: PetscFunctionBegin;
4816: PetscCall(VecGetDM(vec, &dm));
4817: PetscCall(PetscObjectReference((PetscObject)dm));
4818: PetscCall(DMPforestGetPlex(dm, &plex));
4819: PetscCall(VecSetDM(vec, plex));
4820: PetscCall(VecView_Plex_Native(vec, viewer));
4821: PetscCall(VecSetDM(vec, dm));
4822: PetscCall(DMDestroy(&dm));
4823: PetscFunctionReturn(PETSC_SUCCESS);
4824: }
4826: #define VecLoad_pforest _append_pforest(VecLoad)
4827: static PetscErrorCode VecLoad_pforest(Vec vec, PetscViewer viewer)
4828: {
4829: DM dm, plex;
4831: PetscFunctionBegin;
4832: PetscCall(VecGetDM(vec, &dm));
4833: PetscCall(PetscObjectReference((PetscObject)dm));
4834: PetscCall(DMPforestGetPlex(dm, &plex));
4835: PetscCall(VecSetDM(vec, plex));
4836: PetscCall(VecLoad_Plex(vec, viewer));
4837: PetscCall(VecSetDM(vec, dm));
4838: PetscCall(DMDestroy(&dm));
4839: PetscFunctionReturn(PETSC_SUCCESS);
4840: }
4842: #define VecLoad_pforest_Native _infix_pforest(VecLoad, _Native)
4843: static PetscErrorCode VecLoad_pforest_Native(Vec vec, PetscViewer viewer)
4844: {
4845: DM dm, plex;
4847: PetscFunctionBegin;
4848: PetscCall(VecGetDM(vec, &dm));
4849: PetscCall(PetscObjectReference((PetscObject)dm));
4850: PetscCall(DMPforestGetPlex(dm, &plex));
4851: PetscCall(VecSetDM(vec, plex));
4852: PetscCall(VecLoad_Plex_Native(vec, viewer));
4853: PetscCall(VecSetDM(vec, dm));
4854: PetscCall(DMDestroy(&dm));
4855: PetscFunctionReturn(PETSC_SUCCESS);
4856: }
4858: #define DMCreateGlobalVector_pforest _append_pforest(DMCreateGlobalVector)
4859: static PetscErrorCode DMCreateGlobalVector_pforest(DM dm, Vec *vec)
4860: {
4861: PetscFunctionBegin;
4862: PetscCall(DMCreateGlobalVector_Section_Private(dm, vec));
4863: /* PetscCall(VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM)); */
4864: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_pforest));
4865: PetscCall(VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void))VecView_pforest_Native));
4866: PetscCall(VecSetOperation(*vec, VECOP_LOAD, (void (*)(void))VecLoad_pforest));
4867: PetscCall(VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void))VecLoad_pforest_Native));
4868: PetscFunctionReturn(PETSC_SUCCESS);
4869: }
4871: #define DMCreateLocalVector_pforest _append_pforest(DMCreateLocalVector)
4872: static PetscErrorCode DMCreateLocalVector_pforest(DM dm, Vec *vec)
4873: {
4874: PetscFunctionBegin;
4875: PetscCall(DMCreateLocalVector_Section_Private(dm, vec));
4876: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecViewLocal_pforest));
4877: PetscFunctionReturn(PETSC_SUCCESS);
4878: }
4880: #define DMCreateMatrix_pforest _append_pforest(DMCreateMatrix)
4881: static PetscErrorCode DMCreateMatrix_pforest(DM dm, Mat *mat)
4882: {
4883: DM plex;
4885: PetscFunctionBegin;
4887: PetscCall(DMPforestGetPlex(dm, &plex));
4888: if (plex->prealloc_only != dm->prealloc_only) plex->prealloc_only = dm->prealloc_only; /* maybe this should go into forest->plex */
4889: PetscCall(DMSetMatType(plex, dm->mattype));
4890: PetscCall(DMCreateMatrix(plex, mat));
4891: PetscCall(MatSetDM(*mat, dm));
4892: PetscFunctionReturn(PETSC_SUCCESS);
4893: }
4895: #define DMProjectFunctionLocal_pforest _append_pforest(DMProjectFunctionLocal)
4896: static PetscErrorCode DMProjectFunctionLocal_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4897: {
4898: DM plex;
4900: PetscFunctionBegin;
4902: PetscCall(DMPforestGetPlex(dm, &plex));
4903: PetscCall(DMProjectFunctionLocal(plex, time, funcs, ctxs, mode, localX));
4904: PetscFunctionReturn(PETSC_SUCCESS);
4905: }
4907: #define DMProjectFunctionLabelLocal_pforest _append_pforest(DMProjectFunctionLabelLocal)
4908: static PetscErrorCode DMProjectFunctionLabelLocal_pforest(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
4909: {
4910: DM plex;
4912: PetscFunctionBegin;
4914: PetscCall(DMPforestGetPlex(dm, &plex));
4915: PetscCall(DMProjectFunctionLabelLocal(plex, time, label, numIds, ids, Ncc, comps, funcs, ctxs, mode, localX));
4916: PetscFunctionReturn(PETSC_SUCCESS);
4917: }
4919: #define DMProjectFieldLocal_pforest _append_pforest(DMProjectFieldLocal)
4920: PetscErrorCode DMProjectFieldLocal_pforest(DM dm, PetscReal time, Vec localU, void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode mode, Vec localX)
4921: {
4922: DM plex;
4924: PetscFunctionBegin;
4926: PetscCall(DMPforestGetPlex(dm, &plex));
4927: PetscCall(DMProjectFieldLocal(plex, time, localU, funcs, mode, localX));
4928: PetscFunctionReturn(PETSC_SUCCESS);
4929: }
4931: #define DMComputeL2Diff_pforest _append_pforest(DMComputeL2Diff)
4932: PetscErrorCode DMComputeL2Diff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
4933: {
4934: DM plex;
4936: PetscFunctionBegin;
4938: PetscCall(DMPforestGetPlex(dm, &plex));
4939: PetscCall(DMComputeL2Diff(plex, time, funcs, ctxs, X, diff));
4940: PetscFunctionReturn(PETSC_SUCCESS);
4941: }
4943: #define DMComputeL2FieldDiff_pforest _append_pforest(DMComputeL2FieldDiff)
4944: PetscErrorCode DMComputeL2FieldDiff_pforest(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
4945: {
4946: DM plex;
4948: PetscFunctionBegin;
4950: PetscCall(DMPforestGetPlex(dm, &plex));
4951: PetscCall(DMComputeL2FieldDiff(plex, time, funcs, ctxs, X, diff));
4952: PetscFunctionReturn(PETSC_SUCCESS);
4953: }
4955: #define DMCreatelocalsection_pforest _append_pforest(DMCreatelocalsection)
4956: static PetscErrorCode DMCreatelocalsection_pforest(DM dm)
4957: {
4958: DM plex;
4959: PetscSection section;
4961: PetscFunctionBegin;
4963: PetscCall(DMPforestGetPlex(dm, &plex));
4964: PetscCall(DMGetLocalSection(plex, §ion));
4965: PetscCall(DMSetLocalSection(dm, section));
4966: PetscFunctionReturn(PETSC_SUCCESS);
4967: }
4969: #define DMCreateDefaultConstraints_pforest _append_pforest(DMCreateDefaultConstraints)
4970: static PetscErrorCode DMCreateDefaultConstraints_pforest(DM dm)
4971: {
4972: DM plex;
4973: Mat mat;
4974: Vec bias;
4975: PetscSection section;
4977: PetscFunctionBegin;
4979: PetscCall(DMPforestGetPlex(dm, &plex));
4980: PetscCall(DMGetDefaultConstraints(plex, §ion, &mat, &bias));
4981: PetscCall(DMSetDefaultConstraints(dm, section, mat, bias));
4982: PetscFunctionReturn(PETSC_SUCCESS);
4983: }
4985: #define DMGetDimPoints_pforest _append_pforest(DMGetDimPoints)
4986: static PetscErrorCode DMGetDimPoints_pforest(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
4987: {
4988: DM plex;
4990: PetscFunctionBegin;
4992: PetscCall(DMPforestGetPlex(dm, &plex));
4993: PetscCall(DMGetDimPoints(plex, dim, cStart, cEnd));
4994: PetscFunctionReturn(PETSC_SUCCESS);
4995: }
4997: /* Need to forward declare */
4998: #define DMInitialize_pforest _append_pforest(DMInitialize)
4999: static PetscErrorCode DMInitialize_pforest(DM dm);
5001: #define DMClone_pforest _append_pforest(DMClone)
5002: static PetscErrorCode DMClone_pforest(DM dm, DM *newdm)
5003: {
5004: PetscFunctionBegin;
5005: PetscCall(DMClone_Forest(dm, newdm));
5006: PetscCall(DMInitialize_pforest(*newdm));
5007: PetscFunctionReturn(PETSC_SUCCESS);
5008: }
5010: #define DMForestCreateCellChart_pforest _append_pforest(DMForestCreateCellChart)
5011: static PetscErrorCode DMForestCreateCellChart_pforest(DM dm, PetscInt *cStart, PetscInt *cEnd)
5012: {
5013: DM_Forest *forest;
5014: DM_Forest_pforest *pforest;
5015: PetscInt overlap;
5017: PetscFunctionBegin;
5018: PetscCall(DMSetUp(dm));
5019: forest = (DM_Forest *)dm->data;
5020: pforest = (DM_Forest_pforest *)forest->data;
5021: *cStart = 0;
5022: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
5023: if (overlap && pforest->ghost) {
5024: *cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[pforest->forest->mpisize];
5025: } else {
5026: *cEnd = pforest->forest->local_num_quadrants;
5027: }
5028: PetscFunctionReturn(PETSC_SUCCESS);
5029: }
5031: #define DMForestCreateCellSF_pforest _append_pforest(DMForestCreateCellSF)
5032: static PetscErrorCode DMForestCreateCellSF_pforest(DM dm, PetscSF *cellSF)
5033: {
5034: DM_Forest *forest;
5035: DM_Forest_pforest *pforest;
5036: PetscMPIInt rank;
5037: PetscInt overlap;
5038: PetscInt cStart, cEnd, cLocalStart, cLocalEnd;
5039: PetscInt nRoots, nLeaves, *mine = NULL;
5040: PetscSFNode *remote = NULL;
5041: PetscSF sf;
5043: PetscFunctionBegin;
5044: PetscCall(DMForestGetCellChart(dm, &cStart, &cEnd));
5045: forest = (DM_Forest *)dm->data;
5046: pforest = (DM_Forest_pforest *)forest->data;
5047: nRoots = cEnd - cStart;
5048: cLocalStart = pforest->cLocalStart;
5049: cLocalEnd = pforest->cLocalEnd;
5050: nLeaves = 0;
5051: PetscCall(DMForestGetPartitionOverlap(dm, &overlap));
5052: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
5053: if (overlap && pforest->ghost) {
5054: PetscSFNode *mirror;
5055: p4est_quadrant_t *mirror_array;
5056: PetscInt nMirror, nGhostPre, nSelf, q;
5057: void **mirrorPtrs;
5059: nMirror = (PetscInt)pforest->ghost->mirrors.elem_count;
5060: nSelf = cLocalEnd - cLocalStart;
5061: nLeaves = nRoots - nSelf;
5062: nGhostPre = (PetscInt)pforest->ghost->proc_offsets[rank];
5063: PetscCall(PetscMalloc1(nLeaves, &mine));
5064: PetscCall(PetscMalloc1(nLeaves, &remote));
5065: PetscCall(PetscMalloc2(nMirror, &mirror, nMirror, &mirrorPtrs));
5066: mirror_array = (p4est_quadrant_t *)pforest->ghost->mirrors.array;
5067: for (q = 0; q < nMirror; q++) {
5068: p4est_quadrant_t *mir = &mirror_array[q];
5070: mirror[q].rank = rank;
5071: mirror[q].index = (PetscInt)mir->p.piggy3.local_num + cLocalStart;
5072: mirrorPtrs[q] = (void *)&mirror[q];
5073: }
5074: PetscCallP4est(p4est_ghost_exchange_custom, (pforest->forest, pforest->ghost, sizeof(PetscSFNode), mirrorPtrs, remote));
5075: PetscCall(PetscFree2(mirror, mirrorPtrs));
5076: for (q = 0; q < nGhostPre; q++) mine[q] = q;
5077: for (; q < nLeaves; q++) mine[q] = (q - nGhostPre) + cLocalEnd;
5078: }
5079: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf));
5080: PetscCall(PetscSFSetGraph(sf, nRoots, nLeaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
5081: *cellSF = sf;
5082: PetscFunctionReturn(PETSC_SUCCESS);
5083: }
5085: static PetscErrorCode DMCreateNeumannOverlap_pforest(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx)
5086: {
5087: DM plex;
5089: PetscFunctionBegin;
5090: PetscCall(DMPforestGetPlex(dm, &plex));
5091: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5092: PetscCall(DMCreateNeumannOverlap_Plex(plex, ovl, J, setup, setup_ctx));
5093: PetscCall(DMClearAuxiliaryVec(plex));
5094: if (!*setup) {
5095: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup));
5096: if (*setup) PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm));
5097: }
5098: PetscFunctionReturn(PETSC_SUCCESS);
5099: }
5101: #define DMCreateDomainDecomposition_pforest _append_pforest(DMCreateDomainDecomposition)
5102: static PetscErrorCode DMCreateDomainDecomposition_pforest(DM dm, PetscInt *nsub, char ***names, IS **innerises, IS **outerises, DM **dms)
5103: {
5104: DM plex;
5106: PetscFunctionBegin;
5107: PetscCall(DMPforestGetPlex(dm, &plex));
5108: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5109: PetscCall(DMCreateDomainDecomposition(plex, nsub, names, innerises, outerises, dms));
5110: PetscCall(DMClearAuxiliaryVec(plex));
5111: PetscFunctionReturn(PETSC_SUCCESS);
5112: }
5114: #define DMCreateDomainDecompositionScatters_pforest _append_pforest(DMCreateDomainDecompositionScatters)
5115: static PetscErrorCode DMCreateDomainDecompositionScatters_pforest(DM dm, PetscInt n, DM *subdms, VecScatter **iscat, VecScatter **oscat, VecScatter **lscat)
5116: {
5117: DM plex;
5119: PetscFunctionBegin;
5120: PetscCall(DMPforestGetPlex(dm, &plex));
5121: PetscCall(DMCopyAuxiliaryVec(dm, plex));
5122: PetscCall(DMCreateDomainDecompositionScatters(plex, n, subdms, iscat, oscat, lscat));
5123: PetscFunctionReturn(PETSC_SUCCESS);
5124: }
5126: static PetscErrorCode DMInitialize_pforest(DM dm)
5127: {
5128: PetscFunctionBegin;
5129: dm->ops->setup = DMSetUp_pforest;
5130: dm->ops->view = DMView_pforest;
5131: dm->ops->clone = DMClone_pforest;
5132: dm->ops->createinterpolation = DMCreateInterpolation_pforest;
5133: dm->ops->createinjection = DMCreateInjection_pforest;
5134: dm->ops->setfromoptions = DMSetFromOptions_pforest;
5135: dm->ops->createcoordinatedm = DMCreateCoordinateDM_pforest;
5136: dm->ops->createglobalvector = DMCreateGlobalVector_pforest;
5137: dm->ops->createlocalvector = DMCreateLocalVector_pforest;
5138: dm->ops->creatematrix = DMCreateMatrix_pforest;
5139: dm->ops->projectfunctionlocal = DMProjectFunctionLocal_pforest;
5140: dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_pforest;
5141: dm->ops->projectfieldlocal = DMProjectFieldLocal_pforest;
5142: dm->ops->createlocalsection = DMCreatelocalsection_pforest;
5143: dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_pforest;
5144: dm->ops->computel2diff = DMComputeL2Diff_pforest;
5145: dm->ops->computel2fielddiff = DMComputeL2FieldDiff_pforest;
5146: dm->ops->getdimpoints = DMGetDimPoints_pforest;
5147: dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_pforest;
5148: dm->ops->createddscatters = DMCreateDomainDecompositionScatters_pforest;
5150: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", DMConvert_plex_pforest));
5151: PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", DMConvert_pforest_plex));
5152: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_pforest));
5153: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMForestGetPartitionOverlap));
5154: PetscFunctionReturn(PETSC_SUCCESS);
5155: }
5157: #define DMCreate_pforest _append_pforest(DMCreate)
5158: PETSC_EXTERN PetscErrorCode DMCreate_pforest(DM dm)
5159: {
5160: DM_Forest *forest;
5161: DM_Forest_pforest *pforest;
5163: PetscFunctionBegin;
5164: PetscCall(PetscP4estInitialize());
5165: PetscCall(DMCreate_Forest(dm));
5166: PetscCall(DMInitialize_pforest(dm));
5167: PetscCall(DMSetDimension(dm, P4EST_DIM));
5169: /* set forest defaults */
5170: PetscCall(DMForestSetTopology(dm, "unit"));
5171: PetscCall(DMForestSetMinimumRefinement(dm, 0));
5172: PetscCall(DMForestSetInitialRefinement(dm, 0));
5173: PetscCall(DMForestSetMaximumRefinement(dm, P4EST_QMAXLEVEL));
5174: PetscCall(DMForestSetGradeFactor(dm, 2));
5175: PetscCall(DMForestSetAdjacencyDimension(dm, 0));
5176: PetscCall(DMForestSetPartitionOverlap(dm, 0));
5178: /* create p4est data */
5179: PetscCall(PetscNew(&pforest));
5181: forest = (DM_Forest *)dm->data;
5182: forest->data = pforest;
5183: forest->destroy = DMForestDestroy_pforest;
5184: forest->ftemplate = DMForestTemplate_pforest;
5185: forest->transfervec = DMForestTransferVec_pforest;
5186: forest->transfervecfrombase = DMForestTransferVecFromBase_pforest;
5187: forest->createcellchart = DMForestCreateCellChart_pforest;
5188: forest->createcellsf = DMForestCreateCellSF_pforest;
5189: forest->clearadaptivityforest = DMForestClearAdaptivityForest_pforest;
5190: forest->getadaptivitysuccess = DMForestGetAdaptivitySuccess_pforest;
5191: pforest->topo = NULL;
5192: pforest->forest = NULL;
5193: pforest->ghost = NULL;
5194: pforest->lnodes = NULL;
5195: pforest->partition_for_coarsening = PETSC_TRUE;
5196: pforest->coarsen_hierarchy = PETSC_FALSE;
5197: pforest->cLocalStart = -1;
5198: pforest->cLocalEnd = -1;
5199: pforest->labelsFinalized = PETSC_FALSE;
5200: pforest->ghostName = NULL;
5201: PetscFunctionReturn(PETSC_SUCCESS);
5202: }
5204: #endif /* defined(PETSC_HAVE_P4EST) */