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