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