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