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 PetscCall(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_MAX_INT; 714 ctx.maxLevel = 0; 715 ctx.currLevel = 0; 716 ctx.anyChange = PETSC_FALSE; 717 /* sanity check */ 718 PetscCall(DMForestGetAdaptivityForest(dm, &adaptFrom)); 719 PetscCall(DMForestGetBaseDM(dm, &base)); 720 PetscCall(DMForestGetTopology(dm, &topoName)); 721 PetscCheck(adaptFrom || base || topoName, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "A forest needs a topology, a base DM, or a DM to adapt from"); 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 PetscCall(MPIU_Allreduce(&numValues, &numValuesGlobal, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)adaptFrom))); 874 PetscCall(DMLabelGetDefaultValue(adaptLabel, &defaultValue)); 875 if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN_LAST) { /* uniform coarsen of the last level only (equivalent to DM_ADAPT_COARSEN for conforming grids) */ 876 PetscCall(DMForestGetMinimumRefinement(dm, &ctx.minLevel)); 877 PetscCall(DMPforestGetRefinementLevel(dm, &ctx.currLevel)); 878 pforest->forest->user_pointer = (void *)&ctx; 879 PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_currlevel, NULL)); 880 pforest->forest->user_pointer = (void *)dm; 881 PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL)); 882 /* we will have to change the offset after we compute the overlap */ 883 if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL)); 884 } else if (!numValuesGlobal && defaultValue == DM_ADAPT_COARSEN) { /* uniform coarsen */ 885 PetscCall(DMForestGetMinimumRefinement(dm, &ctx.minLevel)); 886 pforest->forest->user_pointer = (void *)&ctx; 887 PetscCallP4est(p4est_coarsen, (pforest->forest, 0, pforest_coarsen_uniform, NULL)); 888 pforest->forest->user_pointer = (void *)dm; 889 PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL)); 890 /* we will have to change the offset after we compute the overlap */ 891 if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), pforest->forest, 0, apforest->forest, apforest->cLocalStart, &coarseToPreFine, NULL)); 892 } else if (!numValuesGlobal && defaultValue == DM_ADAPT_REFINE) { /* uniform refine */ 893 PetscCall(DMForestGetMaximumRefinement(dm, &ctx.maxLevel)); 894 pforest->forest->user_pointer = (void *)&ctx; 895 PetscCallP4est(p4est_refine, (pforest->forest, 0, pforest_refine_uniform, NULL)); 896 pforest->forest->user_pointer = (void *)dm; 897 PetscCallP4est(p4est_balance, (pforest->forest, P4EST_CONNECT_FULL, NULL)); 898 /* we will have to change the offset after we compute the overlap */ 899 if (computeAdaptSF) PetscCall(DMPforestComputeLocalCellTransferSF(PetscObjectComm((PetscObject)dm), apforest->forest, apforest->cLocalStart, pforest->forest, 0, &preCoarseToFine, NULL)); 900 } else if (numValuesGlobal) { 901 p4est_t *p4est = pforest->forest; 902 PetscInt *cellFlags; 903 DMForestAdaptivityStrategy strategy; 904 PetscSF cellSF; 905 PetscInt c, cStart, cEnd; 906 PetscBool adaptAny; 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 if (!forest->cellWeights && forest->weightCapacity == 1. && forest->weightsFactor == 1.) { 1117 PetscCallP4estReturn(shipped, p4est_partition_ext, (pforest->forest, (int)pforest->partition_for_coarsening, NULL)); 1118 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Non-uniform partition cases not implemented yet"); 1119 if (shipped) ctx.anyChange = PETSC_TRUE; 1120 if (forest_copy) { 1121 if (preCoarseToFine || coarseToPreFine) { 1122 PetscSF repartSF; /* repartSF has roots in the old partition */ 1123 PetscInt pStart = -1, pEnd = -1, p; 1124 PetscInt numRoots, numLeaves; 1125 PetscSFNode *repartRoots; 1126 p4est_gloidx_t postStart = pforest->forest->global_first_quadrant[rank]; 1127 p4est_gloidx_t postEnd = pforest->forest->global_first_quadrant[rank + 1]; 1128 p4est_gloidx_t partOffset = postStart; 1129 1130 numRoots = (PetscInt)(forest_copy->global_first_quadrant[rank + 1] - forest_copy->global_first_quadrant[rank]); 1131 numLeaves = (PetscInt)(postEnd - postStart); 1132 PetscCall(DMPforestComputeOverlappingRanks(size, rank, pforest->forest, forest_copy, &pStart, &pEnd)); 1133 PetscCall(PetscMalloc1((PetscInt)pforest->forest->local_num_quadrants, &repartRoots)); 1134 for (p = pStart; p < pEnd; p++) { 1135 p4est_gloidx_t preStart = forest_copy->global_first_quadrant[p]; 1136 p4est_gloidx_t preEnd = forest_copy->global_first_quadrant[p + 1]; 1137 PetscInt q; 1138 1139 if (preEnd == preStart) continue; 1140 PetscCheck(preStart <= postStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad partition overlap computation"); 1141 preEnd = preEnd > postEnd ? postEnd : preEnd; 1142 for (q = partOffset; q < preEnd; q++) { 1143 repartRoots[q - postStart].rank = p; 1144 repartRoots[q - postStart].index = partOffset - preStart; 1145 } 1146 partOffset = preEnd; 1147 } 1148 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &repartSF)); 1149 PetscCall(PetscSFSetGraph(repartSF, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, repartRoots, PETSC_OWN_POINTER)); 1150 PetscCall(PetscSFSetUp(repartSF)); 1151 if (preCoarseToFine) { 1152 PetscSF repartSFembed, preCoarseToFineNew; 1153 PetscInt nleaves; 1154 const PetscInt *leaves; 1155 1156 PetscCall(PetscSFSetUp(preCoarseToFine)); 1157 PetscCall(PetscSFGetGraph(preCoarseToFine, NULL, &nleaves, &leaves, NULL)); 1158 if (leaves) { 1159 PetscCall(PetscSFCreateEmbeddedRootSF(repartSF, nleaves, leaves, &repartSFembed)); 1160 } else { 1161 repartSFembed = repartSF; 1162 PetscCall(PetscObjectReference((PetscObject)repartSFembed)); 1163 } 1164 PetscCall(PetscSFCompose(preCoarseToFine, repartSFembed, &preCoarseToFineNew)); 1165 PetscCall(PetscSFDestroy(&preCoarseToFine)); 1166 PetscCall(PetscSFDestroy(&repartSFembed)); 1167 preCoarseToFine = preCoarseToFineNew; 1168 } 1169 if (coarseToPreFine) { 1170 PetscSF repartSFinv, coarseToPreFineNew; 1171 1172 PetscCall(PetscSFCreateInverseSF(repartSF, &repartSFinv)); 1173 PetscCall(PetscSFCompose(repartSFinv, coarseToPreFine, &coarseToPreFineNew)); 1174 PetscCall(PetscSFDestroy(&coarseToPreFine)); 1175 PetscCall(PetscSFDestroy(&repartSFinv)); 1176 coarseToPreFine = coarseToPreFineNew; 1177 } 1178 PetscCall(PetscSFDestroy(&repartSF)); 1179 } 1180 PetscCallP4est(p4est_destroy, (forest_copy)); 1181 } 1182 } 1183 if (size > 1) { 1184 PetscInt overlap; 1185 1186 PetscCall(DMForestGetPartitionOverlap(dm, &overlap)); 1187 1188 if (adaptFrom) { 1189 PetscInt aoverlap; 1190 1191 PetscCall(DMForestGetPartitionOverlap(adaptFrom, &aoverlap)); 1192 if (aoverlap != overlap) ctx.anyChange = PETSC_TRUE; 1193 } 1194 1195 if (overlap > 0) { 1196 PetscInt i, cLocalStart; 1197 PetscInt cEnd; 1198 PetscSF preCellSF = NULL, cellSF = NULL; 1199 1200 PetscCallP4estReturn(pforest->ghost, p4est_ghost_new, (pforest->forest, P4EST_CONNECT_FULL)); 1201 PetscCallP4estReturn(pforest->lnodes, p4est_lnodes_new, (pforest->forest, pforest->ghost, -P4EST_DIM)); 1202 PetscCallP4est(p4est_ghost_support_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost)); 1203 for (i = 1; i < overlap; i++) PetscCallP4est(p4est_ghost_expand_by_lnodes, (pforest->forest, pforest->lnodes, pforest->ghost)); 1204 1205 cLocalStart = pforest->cLocalStart = pforest->ghost->proc_offsets[rank]; 1206 cEnd = pforest->forest->local_num_quadrants + pforest->ghost->proc_offsets[size]; 1207 1208 /* shift sfs by cLocalStart, expand by cell SFs */ 1209 if (preCoarseToFine || coarseToPreFine) { 1210 if (adaptFrom) PetscCall(DMForestGetCellSF(adaptFrom, &preCellSF)); 1211 dm->setupcalled = PETSC_TRUE; 1212 PetscCall(DMForestGetCellSF(dm, &cellSF)); 1213 } 1214 if (preCoarseToFine) { 1215 PetscSF preCoarseToFineNew; 1216 PetscInt nleaves, nroots, *leavesNew, i, nleavesNew; 1217 const PetscInt *leaves; 1218 const PetscSFNode *remotes; 1219 PetscSFNode *remotesAll; 1220 1221 PetscCall(PetscSFSetUp(preCoarseToFine)); 1222 PetscCall(PetscSFGetGraph(preCoarseToFine, &nroots, &nleaves, &leaves, &remotes)); 1223 PetscCall(PetscMalloc1(cEnd, &remotesAll)); 1224 for (i = 0; i < cEnd; i++) { 1225 remotesAll[i].rank = -1; 1226 remotesAll[i].index = -1; 1227 } 1228 for (i = 0; i < nleaves; i++) remotesAll[(leaves ? leaves[i] : i) + cLocalStart] = remotes[i]; 1229 PetscCall(PetscSFSetUp(cellSF)); 1230 PetscCall(PetscSFBcastBegin(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE)); 1231 PetscCall(PetscSFBcastEnd(cellSF, MPIU_2INT, remotesAll, remotesAll, MPI_REPLACE)); 1232 nleavesNew = 0; 1233 for (i = 0; i < nleaves; i++) { 1234 if (remotesAll[i].rank >= 0) nleavesNew++; 1235 } 1236 PetscCall(PetscMalloc1(nleavesNew, &leavesNew)); 1237 nleavesNew = 0; 1238 for (i = 0; i < nleaves; i++) { 1239 if (remotesAll[i].rank >= 0) { 1240 leavesNew[nleavesNew] = i; 1241 if (i > nleavesNew) remotesAll[nleavesNew] = remotesAll[i]; 1242 nleavesNew++; 1243 } 1244 } 1245 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &preCoarseToFineNew)); 1246 if (nleavesNew < cEnd) { 1247 PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, leavesNew, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES)); 1248 } else { /* all cells are leaves */ 1249 PetscCall(PetscFree(leavesNew)); 1250 PetscCall(PetscSFSetGraph(preCoarseToFineNew, nroots, nleavesNew, NULL, PETSC_OWN_POINTER, remotesAll, PETSC_COPY_VALUES)); 1251 } 1252 PetscCall(PetscFree(remotesAll)); 1253 PetscCall(PetscSFDestroy(&preCoarseToFine)); 1254 preCoarseToFine = preCoarseToFineNew; 1255 preCoarseToFine = preCoarseToFineNew; 1256 } 1257 if (coarseToPreFine) { 1258 PetscSF coarseToPreFineNew; 1259 PetscInt nleaves, nroots, i, nleavesCellSF, nleavesExpanded, *leavesNew; 1260 const PetscInt *leaves; 1261 const PetscSFNode *remotes; 1262 PetscSFNode *remotesNew, *remotesNewRoot, *remotesExpanded; 1263 1264 PetscCall(PetscSFSetUp(coarseToPreFine)); 1265 PetscCall(PetscSFGetGraph(coarseToPreFine, &nroots, &nleaves, &leaves, &remotes)); 1266 PetscCall(PetscSFGetGraph(preCellSF, NULL, &nleavesCellSF, NULL, NULL)); 1267 PetscCall(PetscMalloc1(nroots, &remotesNewRoot)); 1268 PetscCall(PetscMalloc1(nleaves, &remotesNew)); 1269 for (i = 0; i < nroots; i++) { 1270 remotesNewRoot[i].rank = rank; 1271 remotesNewRoot[i].index = i + cLocalStart; 1272 } 1273 PetscCall(PetscSFBcastBegin(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE)); 1274 PetscCall(PetscSFBcastEnd(coarseToPreFine, MPIU_2INT, remotesNewRoot, remotesNew, MPI_REPLACE)); 1275 PetscCall(PetscFree(remotesNewRoot)); 1276 PetscCall(PetscMalloc1(nleavesCellSF, &remotesExpanded)); 1277 for (i = 0; i < nleavesCellSF; i++) { 1278 remotesExpanded[i].rank = -1; 1279 remotesExpanded[i].index = -1; 1280 } 1281 for (i = 0; i < nleaves; i++) remotesExpanded[leaves ? leaves[i] : i] = remotesNew[i]; 1282 PetscCall(PetscFree(remotesNew)); 1283 PetscCall(PetscSFBcastBegin(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE)); 1284 PetscCall(PetscSFBcastEnd(preCellSF, MPIU_2INT, remotesExpanded, remotesExpanded, MPI_REPLACE)); 1285 1286 nleavesExpanded = 0; 1287 for (i = 0; i < nleavesCellSF; i++) { 1288 if (remotesExpanded[i].rank >= 0) nleavesExpanded++; 1289 } 1290 PetscCall(PetscMalloc1(nleavesExpanded, &leavesNew)); 1291 nleavesExpanded = 0; 1292 for (i = 0; i < nleavesCellSF; i++) { 1293 if (remotesExpanded[i].rank >= 0) { 1294 leavesNew[nleavesExpanded] = i; 1295 if (i > nleavesExpanded) remotesExpanded[nleavesExpanded] = remotes[i]; 1296 nleavesExpanded++; 1297 } 1298 } 1299 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), &coarseToPreFineNew)); 1300 if (nleavesExpanded < nleavesCellSF) { 1301 PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, leavesNew, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES)); 1302 } else { 1303 PetscCall(PetscFree(leavesNew)); 1304 PetscCall(PetscSFSetGraph(coarseToPreFineNew, cEnd, nleavesExpanded, NULL, PETSC_OWN_POINTER, remotesExpanded, PETSC_COPY_VALUES)); 1305 } 1306 PetscCall(PetscFree(remotesExpanded)); 1307 PetscCall(PetscSFDestroy(&coarseToPreFine)); 1308 coarseToPreFine = coarseToPreFineNew; 1309 } 1310 } 1311 } 1312 } 1313 forest->preCoarseToFine = preCoarseToFine; 1314 forest->coarseToPreFine = coarseToPreFine; 1315 dm->setupcalled = PETSC_TRUE; 1316 PetscCall(MPIU_Allreduce(&ctx.anyChange, &pforest->adaptivitySuccess, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm))); 1317 PetscCall(DMPforestGetPlex(dm, NULL)); 1318 PetscFunctionReturn(PETSC_SUCCESS); 1319 } 1320 1321 #define DMForestGetAdaptivitySuccess_pforest _append_pforest(DMForestGetAdaptivitySuccess) 1322 static PetscErrorCode DMForestGetAdaptivitySuccess_pforest(DM dm, PetscBool *success) 1323 { 1324 DM_Forest *forest; 1325 DM_Forest_pforest *pforest; 1326 1327 PetscFunctionBegin; 1328 forest = (DM_Forest *)dm->data; 1329 pforest = (DM_Forest_pforest *)forest->data; 1330 *success = pforest->adaptivitySuccess; 1331 PetscFunctionReturn(PETSC_SUCCESS); 1332 } 1333 1334 #define DMView_ASCII_pforest _append_pforest(DMView_ASCII) 1335 static PetscErrorCode DMView_ASCII_pforest(PetscObject odm, PetscViewer viewer) 1336 { 1337 DM dm = (DM)odm; 1338 1339 PetscFunctionBegin; 1340 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1341 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1342 PetscCall(DMSetUp(dm)); 1343 switch (viewer->format) { 1344 case PETSC_VIEWER_DEFAULT: 1345 case PETSC_VIEWER_ASCII_INFO: { 1346 PetscInt dim; 1347 const char *name; 1348 1349 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 1350 PetscCall(DMGetDimension(dm, &dim)); 1351 if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "Forest %s in %" PetscInt_FMT " dimensions:\n", name, dim)); 1352 else PetscCall(PetscViewerASCIIPrintf(viewer, "Forest in %" PetscInt_FMT " dimensions:\n", dim)); 1353 } /* fall through */ 1354 case PETSC_VIEWER_ASCII_INFO_DETAIL: 1355 case PETSC_VIEWER_LOAD_BALANCE: { 1356 DM plex; 1357 1358 PetscCall(DMPforestGetPlex(dm, &plex)); 1359 PetscCall(DMView(plex, viewer)); 1360 } break; 1361 default: 1362 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]); 1363 } 1364 PetscFunctionReturn(PETSC_SUCCESS); 1365 } 1366 1367 #define DMView_VTK_pforest _append_pforest(DMView_VTK) 1368 static PetscErrorCode DMView_VTK_pforest(PetscObject odm, PetscViewer viewer) 1369 { 1370 DM dm = (DM)odm; 1371 DM_Forest *forest = (DM_Forest *)dm->data; 1372 DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data; 1373 PetscBool isvtk; 1374 PetscReal vtkScale = 1. - PETSC_MACHINE_EPSILON; 1375 PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data; 1376 const char *name; 1377 char *filenameStrip = NULL; 1378 PetscBool hasExt; 1379 size_t len; 1380 p4est_geometry_t *geom; 1381 1382 PetscFunctionBegin; 1383 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1384 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1385 PetscCall(DMSetUp(dm)); 1386 geom = pforest->topo->geom; 1387 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk)); 1388 PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name); 1389 switch (viewer->format) { 1390 case PETSC_VIEWER_VTK_VTU: 1391 PetscCheck(pforest->forest, PetscObjectComm(odm), PETSC_ERR_ARG_WRONG, "DM has not been setup with a valid forest"); 1392 name = vtk->filename; 1393 PetscCall(PetscStrlen(name, &len)); 1394 PetscCall(PetscStrcasecmp(name + len - 4, ".vtu", &hasExt)); 1395 if (hasExt) { 1396 PetscCall(PetscStrallocpy(name, &filenameStrip)); 1397 filenameStrip[len - 4] = '\0'; 1398 name = filenameStrip; 1399 } 1400 if (!pforest->topo->geom) PetscCallP4estReturn(geom, p4est_geometry_new_connectivity, (pforest->topo->conn)); 1401 { 1402 p4est_vtk_context_t *pvtk; 1403 int footerr; 1404 1405 PetscCallP4estReturn(pvtk, p4est_vtk_context_new, (pforest->forest, name)); 1406 PetscCallP4est(p4est_vtk_context_set_geom, (pvtk, geom)); 1407 PetscCallP4est(p4est_vtk_context_set_scale, (pvtk, (double)vtkScale)); 1408 PetscCallP4estReturn(pvtk, p4est_vtk_write_header, (pvtk)); 1409 PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_header() failed"); 1410 PetscCallP4estReturn(pvtk, p4est_vtk_write_cell_dataf, 1411 (pvtk, 1, /* write tree */ 1412 1, /* write level */ 1413 1, /* write rank */ 1414 0, /* do not wrap rank */ 1415 0, /* no scalar fields */ 1416 0, /* no vector fields */ 1417 pvtk)); 1418 PetscCheck(pvtk, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_cell_dataf() failed"); 1419 PetscCallP4estReturn(footerr, p4est_vtk_write_footer, (pvtk)); 1420 PetscCheck(!footerr, PetscObjectComm((PetscObject)odm), PETSC_ERR_LIB, P4EST_STRING "_vtk_write_footer() failed"); 1421 } 1422 if (!pforest->topo->geom) PetscCallP4est(p4est_geometry_destroy, (geom)); 1423 PetscCall(PetscFree(filenameStrip)); 1424 break; 1425 default: 1426 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]); 1427 } 1428 PetscFunctionReturn(PETSC_SUCCESS); 1429 } 1430 1431 #define DMView_HDF5_pforest _append_pforest(DMView_HDF5) 1432 static PetscErrorCode DMView_HDF5_pforest(DM dm, PetscViewer viewer) 1433 { 1434 DM plex; 1435 1436 PetscFunctionBegin; 1437 PetscCall(DMSetUp(dm)); 1438 PetscCall(DMPforestGetPlex(dm, &plex)); 1439 PetscCall(DMView(plex, viewer)); 1440 PetscFunctionReturn(PETSC_SUCCESS); 1441 } 1442 1443 #define DMView_GLVis_pforest _append_pforest(DMView_GLVis) 1444 static PetscErrorCode DMView_GLVis_pforest(DM dm, PetscViewer viewer) 1445 { 1446 DM plex; 1447 1448 PetscFunctionBegin; 1449 PetscCall(DMSetUp(dm)); 1450 PetscCall(DMPforestGetPlex(dm, &plex)); 1451 PetscCall(DMView(plex, viewer)); 1452 PetscFunctionReturn(PETSC_SUCCESS); 1453 } 1454 1455 #define DMView_pforest _append_pforest(DMView) 1456 static PetscErrorCode DMView_pforest(DM dm, PetscViewer viewer) 1457 { 1458 PetscBool isascii, isvtk, ishdf5, isglvis; 1459 1460 PetscFunctionBegin; 1461 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1462 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1463 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 1464 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk)); 1465 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 1466 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis)); 1467 if (isascii) { 1468 PetscCall(DMView_ASCII_pforest((PetscObject)dm, viewer)); 1469 } else if (isvtk) { 1470 PetscCall(DMView_VTK_pforest((PetscObject)dm, viewer)); 1471 } else if (ishdf5) { 1472 PetscCall(DMView_HDF5_pforest(dm, viewer)); 1473 } else if (isglvis) { 1474 PetscCall(DMView_GLVis_pforest(dm, viewer)); 1475 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Viewer not supported (not VTK, HDF5, or GLVis)"); 1476 PetscFunctionReturn(PETSC_SUCCESS); 1477 } 1478 1479 static PetscErrorCode PforestConnectivityEnumerateFacets(p4est_connectivity_t *conn, PetscInt **tree_face_to_uniq) 1480 { 1481 PetscInt *ttf, f, t, g, count; 1482 PetscInt numFacets; 1483 1484 PetscFunctionBegin; 1485 numFacets = conn->num_trees * P4EST_FACES; 1486 PetscCall(PetscMalloc1(numFacets, &ttf)); 1487 for (f = 0; f < numFacets; f++) ttf[f] = -1; 1488 for (g = 0, count = 0, t = 0; t < conn->num_trees; t++) { 1489 for (f = 0; f < P4EST_FACES; f++, g++) { 1490 if (ttf[g] == -1) { 1491 PetscInt ng; 1492 1493 ttf[g] = count++; 1494 ng = conn->tree_to_tree[g] * P4EST_FACES + (conn->tree_to_face[g] % P4EST_FACES); 1495 ttf[ng] = ttf[g]; 1496 } 1497 } 1498 } 1499 *tree_face_to_uniq = ttf; 1500 PetscFunctionReturn(PETSC_SUCCESS); 1501 } 1502 1503 static PetscErrorCode DMPlexCreateConnectivity_pforest(DM dm, p4est_connectivity_t **connOut, PetscInt **tree_face_to_uniq) 1504 { 1505 p4est_topidx_t numTrees, numVerts, numCorns, numCtt; 1506 PetscSection ctt; 1507 #if defined(P4_TO_P8) 1508 p4est_topidx_t numEdges, numEtt; 1509 PetscSection ett; 1510 PetscInt eStart, eEnd, e, ettSize; 1511 PetscInt vertOff = 1 + P4EST_FACES + P8EST_EDGES; 1512 PetscInt edgeOff = 1 + P4EST_FACES; 1513 #else 1514 PetscInt vertOff = 1 + P4EST_FACES; 1515 #endif 1516 p4est_connectivity_t *conn; 1517 PetscInt cStart, cEnd, c, vStart, vEnd, v, fStart, fEnd, f; 1518 PetscInt *star = NULL, *closure = NULL, closureSize, starSize, cttSize; 1519 PetscInt *ttf; 1520 1521 PetscFunctionBegin; 1522 /* 1: count objects, allocate */ 1523 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 1524 PetscCall(P4estTopidxCast(cEnd - cStart, &numTrees)); 1525 numVerts = P4EST_CHILDREN * numTrees; 1526 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1527 PetscCall(P4estTopidxCast(vEnd - vStart, &numCorns)); 1528 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ctt)); 1529 PetscCall(PetscSectionSetChart(ctt, vStart, vEnd)); 1530 for (v = vStart; v < vEnd; v++) { 1531 PetscInt s; 1532 1533 PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star)); 1534 for (s = 0; s < starSize; s++) { 1535 PetscInt p = star[2 * s]; 1536 1537 if (p >= cStart && p < cEnd) { 1538 /* we want to count every time cell p references v, so we see how many times it comes up in the closure. This 1539 * only protects against periodicity problems */ 1540 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1541 PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " with wrong closure size %" PetscInt_FMT " != %d", p, closureSize, P4EST_INSUL); 1542 for (c = 0; c < P4EST_CHILDREN; c++) { 1543 PetscInt cellVert = closure[2 * (c + vertOff)]; 1544 1545 PetscCheck(cellVert >= vStart && cellVert < vEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: vertices"); 1546 if (cellVert == v) PetscCall(PetscSectionAddDof(ctt, v, 1)); 1547 } 1548 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1549 } 1550 } 1551 PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star)); 1552 } 1553 PetscCall(PetscSectionSetUp(ctt)); 1554 PetscCall(PetscSectionGetStorageSize(ctt, &cttSize)); 1555 PetscCall(P4estTopidxCast(cttSize, &numCtt)); 1556 #if defined(P4_TO_P8) 1557 PetscCall(DMPlexGetSimplexOrBoxCells(dm, P4EST_DIM - 1, &eStart, &eEnd)); 1558 PetscCall(P4estTopidxCast(eEnd - eStart, &numEdges)); 1559 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &ett)); 1560 PetscCall(PetscSectionSetChart(ett, eStart, eEnd)); 1561 for (e = eStart; e < eEnd; e++) { 1562 PetscInt s; 1563 1564 PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star)); 1565 for (s = 0; s < starSize; s++) { 1566 PetscInt p = star[2 * s]; 1567 1568 if (p >= cStart && p < cEnd) { 1569 /* we want to count every time cell p references e, so we see how many times it comes up in the closure. This 1570 * only protects against periodicity problems */ 1571 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1572 PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell with wrong closure size"); 1573 for (c = 0; c < P8EST_EDGES; c++) { 1574 PetscInt cellEdge = closure[2 * (c + edgeOff)]; 1575 1576 PetscCheck(cellEdge >= eStart && cellEdge < eEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure: edges"); 1577 if (cellEdge == e) PetscCall(PetscSectionAddDof(ett, e, 1)); 1578 } 1579 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1580 } 1581 } 1582 PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star)); 1583 } 1584 PetscCall(PetscSectionSetUp(ett)); 1585 PetscCall(PetscSectionGetStorageSize(ett, &ettSize)); 1586 PetscCall(P4estTopidxCast(ettSize, &numEtt)); 1587 1588 /* This routine allocates space for the arrays, which we fill below */ 1589 PetscCallP4estReturn(conn, p8est_connectivity_new, (numVerts, numTrees, numEdges, numEtt, numCorns, numCtt)); 1590 #else 1591 PetscCallP4estReturn(conn, p4est_connectivity_new, (numVerts, numTrees, numCorns, numCtt)); 1592 #endif 1593 1594 /* 2: visit every face, determine neighboring cells(trees) */ 1595 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 1, &fStart, &fEnd)); 1596 PetscCall(PetscMalloc1((cEnd - cStart) * P4EST_FACES, &ttf)); 1597 for (f = fStart; f < fEnd; f++) { 1598 PetscInt numSupp, s; 1599 PetscInt myFace[2] = {-1, -1}; 1600 PetscInt myOrnt[2] = {PETSC_MIN_INT, PETSC_MIN_INT}; 1601 const PetscInt *supp; 1602 1603 PetscCall(DMPlexGetSupportSize(dm, f, &numSupp)); 1604 PetscCheck(numSupp == 1 || numSupp == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "point %" PetscInt_FMT " has facet with %" PetscInt_FMT " sides: must be 1 or 2 (boundary or conformal)", f, numSupp); 1605 PetscCall(DMPlexGetSupport(dm, f, &supp)); 1606 1607 for (s = 0; s < numSupp; s++) { 1608 PetscInt p = supp[s]; 1609 1610 if (p >= cEnd) { 1611 numSupp--; 1612 if (s) supp = &supp[1 - s]; 1613 break; 1614 } 1615 } 1616 for (s = 0; s < numSupp; s++) { 1617 PetscInt p = supp[s], i; 1618 PetscInt numCone; 1619 DMPolytopeType ct; 1620 const PetscInt *cone; 1621 const PetscInt *ornt; 1622 PetscInt orient = PETSC_MIN_INT; 1623 1624 PetscCall(DMPlexGetConeSize(dm, p, &numCone)); 1625 PetscCheck(numCone == P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " has %" PetscInt_FMT " facets, expect %d", p, numCone, P4EST_FACES); 1626 PetscCall(DMPlexGetCone(dm, p, &cone)); 1627 PetscCall(DMPlexGetCellType(dm, cone[0], &ct)); 1628 PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); 1629 for (i = 0; i < P4EST_FACES; i++) { 1630 if (cone[i] == f) { 1631 orient = DMPolytopeConvertNewOrientation_Internal(ct, ornt[i]); 1632 break; 1633 } 1634 } 1635 PetscCheck(i < P4EST_FACES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " faced %" PetscInt_FMT " mismatch", p, f); 1636 if (p < cStart || p >= cEnd) { 1637 DMPolytopeType ct; 1638 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1639 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "cell %" PetscInt_FMT " (%s) should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], cStart, cEnd); 1640 } 1641 ttf[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = f - fStart; 1642 if (numSupp == 1) { 1643 /* boundary faces indicated by self reference */ 1644 conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = p - cStart; 1645 conn->tree_to_face[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = (int8_t)PetscFaceToP4estFace[i]; 1646 } else { 1647 const PetscInt N = P4EST_CHILDREN / 2; 1648 1649 conn->tree_to_tree[P4EST_FACES * (p - cStart) + PetscFaceToP4estFace[i]] = supp[1 - s] - cStart; 1650 myFace[s] = PetscFaceToP4estFace[i]; 1651 /* get the orientation of cell p in p4est-type closure to facet f, by composing the p4est-closure to 1652 * petsc-closure permutation and the petsc-closure to facet orientation */ 1653 myOrnt[s] = DihedralCompose(N, orient, DMPolytopeConvertNewOrientation_Internal(ct, P4estFaceToPetscOrnt[myFace[s]])); 1654 } 1655 } 1656 if (numSupp == 2) { 1657 for (s = 0; s < numSupp; s++) { 1658 PetscInt p = supp[s]; 1659 PetscInt orntAtoB; 1660 PetscInt p4estOrient; 1661 const PetscInt N = P4EST_CHILDREN / 2; 1662 1663 /* composing the forward permutation with the other cell's inverse permutation gives the self-to-neighbor 1664 * permutation of this cell-facet's cone */ 1665 orntAtoB = DihedralCompose(N, DihedralInvert(N, myOrnt[1 - s]), myOrnt[s]); 1666 1667 /* convert cone-description permutation (i.e., edges around facet) to cap-description permutation (i.e., 1668 * vertices around facet) */ 1669 #if !defined(P4_TO_P8) 1670 p4estOrient = orntAtoB < 0 ? -(orntAtoB + 1) : orntAtoB; 1671 #else 1672 { 1673 PetscInt firstVert = orntAtoB < 0 ? ((-orntAtoB) % N) : orntAtoB; 1674 PetscInt p4estFirstVert = firstVert < 2 ? firstVert : (firstVert ^ 1); 1675 1676 /* swap bits */ 1677 p4estOrient = ((myFace[s] <= myFace[1 - s]) || (orntAtoB < 0)) ? p4estFirstVert : ((p4estFirstVert >> 1) | ((p4estFirstVert & 1) << 1)); 1678 } 1679 #endif 1680 /* encode neighbor face and orientation in tree_to_face per p4est_connectivity standard (see 1681 * p4est_connectivity.h, p8est_connectivity.h) */ 1682 conn->tree_to_face[P4EST_FACES * (p - cStart) + myFace[s]] = (int8_t)myFace[1 - s] + p4estOrient * P4EST_FACES; 1683 } 1684 } 1685 } 1686 1687 #if defined(P4_TO_P8) 1688 /* 3: visit every edge */ 1689 conn->ett_offset[0] = 0; 1690 for (e = eStart; e < eEnd; e++) { 1691 PetscInt off, s; 1692 1693 PetscCall(PetscSectionGetOffset(ett, e, &off)); 1694 conn->ett_offset[e - eStart] = (p4est_topidx_t)off; 1695 PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star)); 1696 for (s = 0; s < starSize; s++) { 1697 PetscInt p = star[2 * s]; 1698 1699 if (p >= cStart && p < cEnd) { 1700 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1701 PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure"); 1702 for (c = 0; c < P8EST_EDGES; c++) { 1703 PetscInt cellEdge = closure[2 * (c + edgeOff)]; 1704 PetscInt cellOrnt = closure[2 * (c + edgeOff) + 1]; 1705 DMPolytopeType ct; 1706 1707 PetscCall(DMPlexGetCellType(dm, cellEdge, &ct)); 1708 cellOrnt = DMPolytopeConvertNewOrientation_Internal(ct, cellOrnt); 1709 if (cellEdge == e) { 1710 PetscInt p4estEdge = PetscEdgeToP4estEdge[c]; 1711 PetscInt totalOrient; 1712 1713 /* compose p4est-closure to petsc-closure permutation and petsc-closure to edge orientation */ 1714 totalOrient = DihedralCompose(2, cellOrnt, DMPolytopeConvertNewOrientation_Internal(DM_POLYTOPE_SEGMENT, P4estEdgeToPetscOrnt[p4estEdge])); 1715 /* p4est orientations are positive: -2 => 1, -1 => 0 */ 1716 totalOrient = (totalOrient < 0) ? -(totalOrient + 1) : totalOrient; 1717 conn->edge_to_tree[off] = (p4est_locidx_t)(p - cStart); 1718 /* encode cell-edge and orientation in edge_to_edge per p8est_connectivity standard (see 1719 * p8est_connectivity.h) */ 1720 conn->edge_to_edge[off++] = (int8_t)p4estEdge + P8EST_EDGES * totalOrient; 1721 conn->tree_to_edge[P8EST_EDGES * (p - cStart) + p4estEdge] = e - eStart; 1722 } 1723 } 1724 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1725 } 1726 } 1727 PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &starSize, &star)); 1728 } 1729 PetscCall(PetscSectionDestroy(&ett)); 1730 #endif 1731 1732 /* 4: visit every vertex */ 1733 conn->ctt_offset[0] = 0; 1734 for (v = vStart; v < vEnd; v++) { 1735 PetscInt off, s; 1736 1737 PetscCall(PetscSectionGetOffset(ctt, v, &off)); 1738 conn->ctt_offset[v - vStart] = (p4est_topidx_t)off; 1739 PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star)); 1740 for (s = 0; s < starSize; s++) { 1741 PetscInt p = star[2 * s]; 1742 1743 if (p >= cStart && p < cEnd) { 1744 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1745 PetscCheck(closureSize == P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Non-standard closure"); 1746 for (c = 0; c < P4EST_CHILDREN; c++) { 1747 PetscInt cellVert = closure[2 * (c + vertOff)]; 1748 1749 if (cellVert == v) { 1750 PetscInt p4estVert = PetscVertToP4estVert[c]; 1751 1752 conn->corner_to_tree[off] = (p4est_locidx_t)(p - cStart); 1753 conn->corner_to_corner[off++] = (int8_t)p4estVert; 1754 conn->tree_to_corner[P4EST_CHILDREN * (p - cStart) + p4estVert] = v - vStart; 1755 } 1756 } 1757 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure)); 1758 } 1759 } 1760 PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star)); 1761 } 1762 PetscCall(PetscSectionDestroy(&ctt)); 1763 1764 /* 5: Compute the coordinates */ 1765 { 1766 PetscInt coordDim; 1767 1768 PetscCall(DMGetCoordinateDim(dm, &coordDim)); 1769 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 1770 for (c = cStart; c < cEnd; c++) { 1771 PetscInt dof; 1772 PetscBool isDG; 1773 PetscScalar *cellCoords = NULL; 1774 const PetscScalar *array; 1775 1776 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords)); 1777 PetscCheck(dof == P4EST_CHILDREN * coordDim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need coordinates at the corners: (dof) %" PetscInt_FMT " != %d * %" PetscInt_FMT " (sdim)", dof, P4EST_CHILDREN, coordDim); 1778 for (v = 0; v < P4EST_CHILDREN; v++) { 1779 PetscInt i, lim = PetscMin(3, coordDim); 1780 PetscInt p4estVert = PetscVertToP4estVert[v]; 1781 1782 conn->tree_to_vertex[P4EST_CHILDREN * (c - cStart) + v] = P4EST_CHILDREN * (c - cStart) + v; 1783 /* p4est vertices are always embedded in R^3 */ 1784 for (i = 0; i < 3; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = 0.; 1785 for (i = 0; i < lim; i++) conn->vertices[3 * (P4EST_CHILDREN * (c - cStart) + p4estVert) + i] = PetscRealPart(cellCoords[v * coordDim + i]); 1786 } 1787 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &dof, &array, &cellCoords)); 1788 } 1789 } 1790 1791 #if defined(P4EST_ENABLE_DEBUG) 1792 PetscCheck(p4est_connectivity_is_valid(conn), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Plex to p4est conversion failed"); 1793 #endif 1794 1795 *connOut = conn; 1796 1797 *tree_face_to_uniq = ttf; 1798 PetscFunctionReturn(PETSC_SUCCESS); 1799 } 1800 1801 static PetscErrorCode locidx_to_PetscInt(sc_array_t *array) 1802 { 1803 sc_array_t *newarray; 1804 size_t zz, count = array->elem_count; 1805 1806 PetscFunctionBegin; 1807 PetscCheck(array->elem_size == sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size"); 1808 1809 if (sizeof(p4est_locidx_t) == sizeof(PetscInt)) PetscFunctionReturn(PETSC_SUCCESS); 1810 1811 newarray = sc_array_new_size(sizeof(PetscInt), array->elem_count); 1812 for (zz = 0; zz < count; zz++) { 1813 p4est_locidx_t il = *((p4est_locidx_t *)sc_array_index(array, zz)); 1814 PetscInt *ip = (PetscInt *)sc_array_index(newarray, zz); 1815 1816 *ip = (PetscInt)il; 1817 } 1818 1819 sc_array_reset(array); 1820 sc_array_init_size(array, sizeof(PetscInt), count); 1821 sc_array_copy(array, newarray); 1822 sc_array_destroy(newarray); 1823 PetscFunctionReturn(PETSC_SUCCESS); 1824 } 1825 1826 static PetscErrorCode coords_double_to_PetscScalar(sc_array_t *array, PetscInt dim) 1827 { 1828 sc_array_t *newarray; 1829 size_t zz, count = array->elem_count; 1830 1831 PetscFunctionBegin; 1832 PetscCheck(array->elem_size == 3 * sizeof(double), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong coordinate size"); 1833 #if !defined(PETSC_USE_COMPLEX) 1834 if (sizeof(double) == sizeof(PetscScalar) && dim == 3) PetscFunctionReturn(PETSC_SUCCESS); 1835 #endif 1836 1837 newarray = sc_array_new_size(dim * sizeof(PetscScalar), array->elem_count); 1838 for (zz = 0; zz < count; zz++) { 1839 int i; 1840 double *id = (double *)sc_array_index(array, zz); 1841 PetscScalar *ip = (PetscScalar *)sc_array_index(newarray, zz); 1842 1843 for (i = 0; i < dim; i++) ip[i] = 0.; 1844 for (i = 0; i < PetscMin(dim, 3); i++) ip[i] = (PetscScalar)id[i]; 1845 } 1846 1847 sc_array_reset(array); 1848 sc_array_init_size(array, dim * sizeof(PetscScalar), count); 1849 sc_array_copy(array, newarray); 1850 sc_array_destroy(newarray); 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 static PetscErrorCode locidx_pair_to_PetscSFNode(sc_array_t *array) 1855 { 1856 sc_array_t *newarray; 1857 size_t zz, count = array->elem_count; 1858 1859 PetscFunctionBegin; 1860 PetscCheck(array->elem_size == 2 * sizeof(p4est_locidx_t), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong locidx size"); 1861 1862 newarray = sc_array_new_size(sizeof(PetscSFNode), array->elem_count); 1863 for (zz = 0; zz < count; zz++) { 1864 p4est_locidx_t *il = (p4est_locidx_t *)sc_array_index(array, zz); 1865 PetscSFNode *ip = (PetscSFNode *)sc_array_index(newarray, zz); 1866 1867 ip->rank = (PetscInt)il[0]; 1868 ip->index = (PetscInt)il[1]; 1869 } 1870 1871 sc_array_reset(array); 1872 sc_array_init_size(array, sizeof(PetscSFNode), count); 1873 sc_array_copy(array, newarray); 1874 sc_array_destroy(newarray); 1875 PetscFunctionReturn(PETSC_SUCCESS); 1876 } 1877 1878 static PetscErrorCode P4estToPlex_Local(p4est_t *p4est, DM *plex) 1879 { 1880 PetscFunctionBegin; 1881 { 1882 sc_array_t *points_per_dim = sc_array_new(sizeof(p4est_locidx_t)); 1883 sc_array_t *cone_sizes = sc_array_new(sizeof(p4est_locidx_t)); 1884 sc_array_t *cones = sc_array_new(sizeof(p4est_locidx_t)); 1885 sc_array_t *cone_orientations = sc_array_new(sizeof(p4est_locidx_t)); 1886 sc_array_t *coords = sc_array_new(3 * sizeof(double)); 1887 sc_array_t *children = sc_array_new(sizeof(p4est_locidx_t)); 1888 sc_array_t *parents = sc_array_new(sizeof(p4est_locidx_t)); 1889 sc_array_t *childids = sc_array_new(sizeof(p4est_locidx_t)); 1890 sc_array_t *leaves = sc_array_new(sizeof(p4est_locidx_t)); 1891 sc_array_t *remotes = sc_array_new(2 * sizeof(p4est_locidx_t)); 1892 p4est_locidx_t first_local_quad; 1893 1894 PetscCallP4est(p4est_get_plex_data, (p4est, P4EST_CONNECT_FULL, 0, &first_local_quad, points_per_dim, cone_sizes, cones, cone_orientations, coords, children, parents, childids, leaves, remotes)); 1895 1896 PetscCall(locidx_to_PetscInt(points_per_dim)); 1897 PetscCall(locidx_to_PetscInt(cone_sizes)); 1898 PetscCall(locidx_to_PetscInt(cones)); 1899 PetscCall(locidx_to_PetscInt(cone_orientations)); 1900 PetscCall(coords_double_to_PetscScalar(coords, P4EST_DIM)); 1901 1902 PetscCall(DMPlexCreate(PETSC_COMM_SELF, plex)); 1903 PetscCall(DMSetDimension(*plex, P4EST_DIM)); 1904 PetscCall(DMPlexCreateFromDAG(*plex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array)); 1905 PetscCall(DMPlexConvertOldOrientations_Internal(*plex)); 1906 sc_array_destroy(points_per_dim); 1907 sc_array_destroy(cone_sizes); 1908 sc_array_destroy(cones); 1909 sc_array_destroy(cone_orientations); 1910 sc_array_destroy(coords); 1911 sc_array_destroy(children); 1912 sc_array_destroy(parents); 1913 sc_array_destroy(childids); 1914 sc_array_destroy(leaves); 1915 sc_array_destroy(remotes); 1916 } 1917 PetscFunctionReturn(PETSC_SUCCESS); 1918 } 1919 1920 #define DMReferenceTreeGetChildSymmetry_pforest _append_pforest(DMReferenceTreeGetChildSymmetry) 1921 static PetscErrorCode DMReferenceTreeGetChildSymmetry_pforest(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 1922 { 1923 PetscInt coneSize, dStart, dEnd, vStart, vEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 1924 1925 PetscFunctionBegin; 1926 if (parentOrientA == parentOrientB) { 1927 if (childOrientB) *childOrientB = childOrientA; 1928 if (childB) *childB = childA; 1929 PetscFunctionReturn(PETSC_SUCCESS); 1930 } 1931 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1932 if (childA >= vStart && childA < vEnd) { /* vertices (always in the middle) are invariant under rotation */ 1933 if (childOrientB) *childOrientB = 0; 1934 if (childB) *childB = childA; 1935 PetscFunctionReturn(PETSC_SUCCESS); 1936 } 1937 for (dim = 0; dim < 3; dim++) { 1938 PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd)); 1939 if (parent >= dStart && parent <= dEnd) break; 1940 } 1941 PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim); 1942 PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children"); 1943 if (childA < dStart || childA >= dEnd) { /* a 1-cell in a 2-cell */ 1944 /* this is a lower-dimensional child: bootstrap */ 1945 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 1946 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 1947 1948 PetscCall(DMPlexGetSupportSize(dm, childA, &size)); 1949 PetscCall(DMPlexGetSupport(dm, childA, &supp)); 1950 1951 /* find a point sA in supp(childA) that has the same parent */ 1952 for (i = 0; i < size; i++) { 1953 PetscInt sParent; 1954 1955 sA = supp[i]; 1956 if (sA == parent) continue; 1957 PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL)); 1958 if (sParent == parent) break; 1959 } 1960 PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children"); 1961 /* find out which point sB is in an equivalent position to sA under 1962 * parentOrientB */ 1963 PetscCall(DMReferenceTreeGetChildSymmetry_pforest(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB)); 1964 PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize)); 1965 PetscCall(DMPlexGetCone(dm, sA, &coneA)); 1966 PetscCall(DMPlexGetCone(dm, sB, &coneB)); 1967 PetscCall(DMPlexGetConeOrientation(dm, sA, &oA)); 1968 PetscCall(DMPlexGetConeOrientation(dm, sB, &oB)); 1969 /* step through the cone of sA in natural order */ 1970 for (i = 0; i < sConeSize; i++) { 1971 if (coneA[i] == childA) { 1972 /* if childA is at position i in coneA, 1973 * then we want the point that is at sOrientB*i in coneB */ 1974 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize); 1975 if (childB) *childB = coneB[j]; 1976 if (childOrientB) { 1977 DMPolytopeType ct; 1978 PetscInt oBtrue; 1979 1980 PetscCall(DMPlexGetConeSize(dm, childA, &coneSize)); 1981 /* compose sOrientB and oB[j] */ 1982 PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge"); 1983 ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT; 1984 /* we may have to flip an edge */ 1985 oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientationInv(ct, -1, oB[j]); 1986 oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue); 1987 ABswap = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue); 1988 *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap); 1989 } 1990 break; 1991 } 1992 } 1993 PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch"); 1994 PetscFunctionReturn(PETSC_SUCCESS); 1995 } 1996 /* get the cone size and symmetry swap */ 1997 PetscCall(DMPlexGetConeSize(dm, parent, &coneSize)); 1998 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 1999 if (dim == 2) { 2000 /* orientations refer to cones: we want them to refer to vertices: 2001 * if it's a rotation, they are the same, but if the order is reversed, a 2002 * permutation that puts side i first does *not* put vertex i first */ 2003 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 2004 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 2005 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 2006 } else { 2007 oAvert = parentOrientA; 2008 oBvert = parentOrientB; 2009 ABswapVert = ABswap; 2010 } 2011 if (childB) { 2012 /* assume that each child corresponds to a vertex, in the same order */ 2013 PetscInt p, posA = -1, numChildren, i; 2014 const PetscInt *children; 2015 2016 /* count which position the child is in */ 2017 PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children)); 2018 for (i = 0; i < numChildren; i++) { 2019 p = children[i]; 2020 if (p == childA) { 2021 if (dim == 1) { 2022 posA = i; 2023 } else { /* 2D Morton to rotation */ 2024 posA = (i & 2) ? (i ^ 1) : i; 2025 } 2026 break; 2027 } 2028 } 2029 if (posA >= coneSize) { 2030 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find childA in children of parent"); 2031 } else { 2032 /* figure out position B by applying ABswapVert */ 2033 PetscInt posB, childIdB; 2034 2035 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize); 2036 if (dim == 1) { 2037 childIdB = posB; 2038 } else { /* 2D rotation to Morton */ 2039 childIdB = (posB & 2) ? (posB ^ 1) : posB; 2040 } 2041 if (childB) *childB = children[childIdB]; 2042 } 2043 } 2044 if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap); 2045 PetscFunctionReturn(PETSC_SUCCESS); 2046 } 2047 2048 #define DMCreateReferenceTree_pforest _append_pforest(DMCreateReferenceTree) 2049 static PetscErrorCode DMCreateReferenceTree_pforest(MPI_Comm comm, DM *dm) 2050 { 2051 p4est_connectivity_t *refcube; 2052 p4est_t *root, *refined; 2053 DM dmRoot, dmRefined; 2054 DM_Plex *mesh; 2055 PetscMPIInt rank; 2056 #if defined(PETSC_HAVE_MPIUNI) 2057 sc_MPI_Comm comm_self = sc_MPI_COMM_SELF; 2058 #else 2059 MPI_Comm comm_self = PETSC_COMM_SELF; 2060 #endif 2061 2062 PetscFunctionBegin; 2063 PetscCallP4estReturn(refcube, p4est_connectivity_new_byname, ("unit")); 2064 { /* [-1,1]^d geometry */ 2065 PetscInt i, j; 2066 2067 for (i = 0; i < P4EST_CHILDREN; i++) { 2068 for (j = 0; j < 3; j++) { 2069 refcube->vertices[3 * i + j] *= 2.; 2070 refcube->vertices[3 * i + j] -= 1.; 2071 } 2072 } 2073 } 2074 PetscCallP4estReturn(root, p4est_new, (comm_self, refcube, 0, NULL, NULL)); 2075 PetscCallP4estReturn(refined, p4est_new_ext, (comm_self, refcube, 0, 1, 1, 0, NULL, NULL)); 2076 PetscCall(P4estToPlex_Local(root, &dmRoot)); 2077 PetscCall(P4estToPlex_Local(refined, &dmRefined)); 2078 { 2079 #if !defined(P4_TO_P8) 2080 PetscInt nPoints = 25; 2081 PetscInt perm[25] = {0, 1, 2, 3, 4, 12, 8, 14, 6, 9, 15, 5, 13, 10, 7, 11, 16, 22, 20, 24, 17, 21, 18, 23, 19}; 2082 PetscInt ident[25] = {0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0, 0, 0, 5, 6, 7, 8, 1, 2, 3, 4, 0}; 2083 #else 2084 PetscInt nPoints = 125; 2085 PetscInt perm[125] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 32, 16, 36, 24, 40, 12, 17, 37, 25, 41, 9, 33, 20, 26, 42, 13, 21, 27, 43, 10, 34, 18, 38, 28, 14, 19, 39, 29, 11, 35, 22, 30, 15, 2086 23, 31, 44, 84, 76, 92, 52, 86, 68, 94, 60, 78, 70, 96, 45, 85, 77, 93, 54, 72, 62, 74, 46, 80, 53, 87, 69, 95, 64, 82, 47, 81, 55, 73, 66, 48, 88, 56, 90, 61, 79, 71, 2087 97, 49, 89, 58, 63, 75, 50, 57, 91, 65, 83, 51, 59, 67, 98, 106, 110, 122, 114, 120, 118, 124, 99, 111, 115, 119, 100, 107, 116, 121, 101, 117, 102, 108, 112, 123, 103, 113, 104, 109, 105}; 2088 PetscInt ident[125] = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 2089 16, 17, 17, 18, 18, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 0, 0, 0, 0, 0, 0, 19, 20, 21, 22, 23, 24, 25, 26, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1, 2, 3, 4, 5, 6, 0}; 2090 2091 #endif 2092 IS permIS; 2093 DM dmPerm; 2094 2095 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nPoints, perm, PETSC_USE_POINTER, &permIS)); 2096 PetscCall(DMPlexPermute(dmRefined, permIS, &dmPerm)); 2097 if (dmPerm) { 2098 PetscCall(DMDestroy(&dmRefined)); 2099 dmRefined = dmPerm; 2100 } 2101 PetscCall(ISDestroy(&permIS)); 2102 { 2103 PetscInt p; 2104 PetscCall(DMCreateLabel(dmRoot, "identity")); 2105 PetscCall(DMCreateLabel(dmRefined, "identity")); 2106 for (p = 0; p < P4EST_INSUL; p++) PetscCall(DMSetLabelValue(dmRoot, "identity", p, p)); 2107 for (p = 0; p < nPoints; p++) PetscCall(DMSetLabelValue(dmRefined, "identity", p, ident[p])); 2108 } 2109 } 2110 PetscCall(DMPlexCreateReferenceTree_Union(dmRoot, dmRefined, "identity", dm)); 2111 mesh = (DM_Plex *)(*dm)->data; 2112 mesh->getchildsymmetry = DMReferenceTreeGetChildSymmetry_pforest; 2113 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2114 if (rank == 0) { 2115 PetscCall(DMViewFromOptions(dmRoot, NULL, "-dm_p4est_ref_root_view")); 2116 PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_refined_view")); 2117 PetscCall(DMViewFromOptions(dmRefined, NULL, "-dm_p4est_ref_tree_view")); 2118 } 2119 PetscCall(DMDestroy(&dmRefined)); 2120 PetscCall(DMDestroy(&dmRoot)); 2121 PetscCallP4est(p4est_destroy, (refined)); 2122 PetscCallP4est(p4est_destroy, (root)); 2123 PetscCallP4est(p4est_connectivity_destroy, (refcube)); 2124 PetscFunctionReturn(PETSC_SUCCESS); 2125 } 2126 2127 static PetscErrorCode DMShareDiscretization(DM dmA, DM dmB) 2128 { 2129 void *ctx; 2130 PetscInt num; 2131 PetscReal val; 2132 2133 PetscFunctionBegin; 2134 PetscCall(DMGetApplicationContext(dmA, &ctx)); 2135 PetscCall(DMSetApplicationContext(dmB, ctx)); 2136 PetscCall(DMCopyDisc(dmA, dmB)); 2137 PetscCall(DMGetOutputSequenceNumber(dmA, &num, &val)); 2138 PetscCall(DMSetOutputSequenceNumber(dmB, num, val)); 2139 if (dmB->localSection != dmA->localSection || dmB->globalSection != dmA->globalSection) { 2140 PetscCall(DMClearLocalVectors(dmB)); 2141 PetscCall(PetscObjectReference((PetscObject)dmA->localSection)); 2142 PetscCall(PetscSectionDestroy(&dmB->localSection)); 2143 dmB->localSection = dmA->localSection; 2144 PetscCall(DMClearGlobalVectors(dmB)); 2145 PetscCall(PetscObjectReference((PetscObject)dmA->globalSection)); 2146 PetscCall(PetscSectionDestroy(&dmB->globalSection)); 2147 dmB->globalSection = dmA->globalSection; 2148 PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.section)); 2149 PetscCall(PetscSectionDestroy(&dmB->defaultConstraint.section)); 2150 dmB->defaultConstraint.section = dmA->defaultConstraint.section; 2151 PetscCall(PetscObjectReference((PetscObject)dmA->defaultConstraint.mat)); 2152 PetscCall(MatDestroy(&dmB->defaultConstraint.mat)); 2153 dmB->defaultConstraint.mat = dmA->defaultConstraint.mat; 2154 if (dmA->map) PetscCall(PetscLayoutReference(dmA->map, &dmB->map)); 2155 } 2156 if (dmB->sectionSF != dmA->sectionSF) { 2157 PetscCall(PetscObjectReference((PetscObject)dmA->sectionSF)); 2158 PetscCall(PetscSFDestroy(&dmB->sectionSF)); 2159 dmB->sectionSF = dmA->sectionSF; 2160 } 2161 PetscFunctionReturn(PETSC_SUCCESS); 2162 } 2163 2164 /* Get an SF that broadcasts a coarse-cell covering of the local fine cells */ 2165 static PetscErrorCode DMPforestGetCellCoveringSF(MPI_Comm comm, p4est_t *p4estC, p4est_t *p4estF, PetscInt cStart, PetscInt cEnd, PetscSF *coveringSF) 2166 { 2167 PetscInt startF, endF, startC, endC, p, nLeaves; 2168 PetscSFNode *leaves; 2169 PetscSF sf; 2170 PetscInt *recv, *send; 2171 PetscMPIInt tag; 2172 MPI_Request *recvReqs, *sendReqs; 2173 PetscSection section; 2174 2175 PetscFunctionBegin; 2176 PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estF, p4estC, &startC, &endC)); 2177 PetscCall(PetscMalloc2(2 * (endC - startC), &recv, endC - startC, &recvReqs)); 2178 PetscCall(PetscCommGetNewTag(comm, &tag)); 2179 for (p = startC; p < endC; p++) { 2180 recvReqs[p - startC] = MPI_REQUEST_NULL; /* just in case we don't initiate a receive */ 2181 if (p4estC->global_first_quadrant[p] == p4estC->global_first_quadrant[p + 1]) { /* empty coarse partition */ 2182 recv[2 * (p - startC)] = 0; 2183 recv[2 * (p - startC) + 1] = 0; 2184 continue; 2185 } 2186 2187 PetscCallMPI(MPI_Irecv(&recv[2 * (p - startC)], 2, MPIU_INT, p, tag, comm, &recvReqs[p - startC])); 2188 } 2189 PetscCall(DMPforestComputeOverlappingRanks(p4estC->mpisize, p4estC->mpirank, p4estC, p4estF, &startF, &endF)); 2190 PetscCall(PetscMalloc2(2 * (endF - startF), &send, endF - startF, &sendReqs)); 2191 /* count the quadrants rank will send to each of [startF,endF) */ 2192 for (p = startF; p < endF; p++) { 2193 p4est_quadrant_t *myFineStart = &p4estF->global_first_position[p]; 2194 p4est_quadrant_t *myFineEnd = &p4estF->global_first_position[p + 1]; 2195 PetscInt tStart = (PetscInt)myFineStart->p.which_tree; 2196 PetscInt tEnd = (PetscInt)myFineEnd->p.which_tree; 2197 PetscInt firstCell = -1, lastCell = -1; 2198 p4est_tree_t *treeStart = &(((p4est_tree_t *)p4estC->trees->array)[tStart]); 2199 p4est_tree_t *treeEnd = (size_t)tEnd < p4estC->trees->elem_count ? &(((p4est_tree_t *)p4estC->trees->array)[tEnd]) : NULL; 2200 ssize_t overlapIndex; 2201 2202 sendReqs[p - startF] = MPI_REQUEST_NULL; /* just in case we don't initiate a send */ 2203 if (p4estF->global_first_quadrant[p] == p4estF->global_first_quadrant[p + 1]) continue; 2204 2205 /* locate myFineStart in (or before) a cell */ 2206 if (treeStart->quadrants.elem_count) { 2207 PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeStart->quadrants, myFineStart, p4est_quadrant_disjoint)); 2208 if (overlapIndex < 0) { 2209 firstCell = 0; 2210 } else { 2211 firstCell = treeStart->quadrants_offset + overlapIndex; 2212 } 2213 } else { 2214 firstCell = 0; 2215 } 2216 if (treeEnd && treeEnd->quadrants.elem_count) { 2217 PetscCallP4estReturn(overlapIndex, sc_array_bsearch, (&treeEnd->quadrants, myFineEnd, p4est_quadrant_disjoint)); 2218 if (overlapIndex < 0) { /* all of this local section is overlapped */ 2219 lastCell = p4estC->local_num_quadrants; 2220 } else { 2221 p4est_quadrant_t *container = &(((p4est_quadrant_t *)treeEnd->quadrants.array)[overlapIndex]); 2222 p4est_quadrant_t first_desc; 2223 int equal; 2224 2225 PetscCallP4est(p4est_quadrant_first_descendant, (container, &first_desc, P4EST_QMAXLEVEL)); 2226 PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (myFineEnd, &first_desc)); 2227 if (equal) { 2228 lastCell = treeEnd->quadrants_offset + overlapIndex; 2229 } else { 2230 lastCell = treeEnd->quadrants_offset + overlapIndex + 1; 2231 } 2232 } 2233 } else { 2234 lastCell = p4estC->local_num_quadrants; 2235 } 2236 send[2 * (p - startF)] = firstCell; 2237 send[2 * (p - startF) + 1] = lastCell - firstCell; 2238 PetscCallMPI(MPI_Isend(&send[2 * (p - startF)], 2, MPIU_INT, p, tag, comm, &sendReqs[p - startF])); 2239 } 2240 PetscCallMPI(MPI_Waitall((PetscMPIInt)(endC - startC), recvReqs, MPI_STATUSES_IGNORE)); 2241 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion)); 2242 PetscCall(PetscSectionSetChart(section, startC, endC)); 2243 for (p = startC; p < endC; p++) { 2244 PetscInt numCells = recv[2 * (p - startC) + 1]; 2245 PetscCall(PetscSectionSetDof(section, p, numCells)); 2246 } 2247 PetscCall(PetscSectionSetUp(section)); 2248 PetscCall(PetscSectionGetStorageSize(section, &nLeaves)); 2249 PetscCall(PetscMalloc1(nLeaves, &leaves)); 2250 for (p = startC; p < endC; p++) { 2251 PetscInt firstCell = recv[2 * (p - startC)]; 2252 PetscInt numCells = recv[2 * (p - startC) + 1]; 2253 PetscInt off, i; 2254 2255 PetscCall(PetscSectionGetOffset(section, p, &off)); 2256 for (i = 0; i < numCells; i++) { 2257 leaves[off + i].rank = p; 2258 leaves[off + i].index = firstCell + i; 2259 } 2260 } 2261 PetscCall(PetscSFCreate(comm, &sf)); 2262 PetscCall(PetscSFSetGraph(sf, cEnd - cStart, nLeaves, NULL, PETSC_OWN_POINTER, leaves, PETSC_OWN_POINTER)); 2263 PetscCall(PetscSectionDestroy(§ion)); 2264 PetscCallMPI(MPI_Waitall((PetscMPIInt)(endF - startF), sendReqs, MPI_STATUSES_IGNORE)); 2265 PetscCall(PetscFree2(send, sendReqs)); 2266 PetscCall(PetscFree2(recv, recvReqs)); 2267 *coveringSF = sf; 2268 PetscFunctionReturn(PETSC_SUCCESS); 2269 } 2270 2271 /* closure points for locally-owned cells */ 2272 static PetscErrorCode DMPforestGetCellSFNodes(DM dm, PetscInt numClosureIndices, PetscInt *numClosurePoints, PetscSFNode **closurePoints, PetscBool redirect) 2273 { 2274 PetscInt cStart, cEnd; 2275 PetscInt count, c; 2276 PetscMPIInt rank; 2277 PetscInt closureSize = -1; 2278 PetscInt *closure = NULL; 2279 PetscSF pointSF; 2280 PetscInt nleaves, nroots; 2281 const PetscInt *ilocal; 2282 const PetscSFNode *iremote; 2283 DM plex; 2284 DM_Forest *forest; 2285 DM_Forest_pforest *pforest; 2286 2287 PetscFunctionBegin; 2288 forest = (DM_Forest *)dm->data; 2289 pforest = (DM_Forest_pforest *)forest->data; 2290 cStart = pforest->cLocalStart; 2291 cEnd = pforest->cLocalEnd; 2292 PetscCall(DMPforestGetPlex(dm, &plex)); 2293 PetscCall(DMGetPointSF(dm, &pointSF)); 2294 PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote)); 2295 nleaves = PetscMax(0, nleaves); 2296 nroots = PetscMax(0, nroots); 2297 *numClosurePoints = numClosureIndices * (cEnd - cStart); 2298 PetscCall(PetscMalloc1(*numClosurePoints, closurePoints)); 2299 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2300 for (c = cStart, count = 0; c < cEnd; c++) { 2301 PetscInt i; 2302 PetscCall(DMPlexGetTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure)); 2303 2304 for (i = 0; i < numClosureIndices; i++, count++) { 2305 PetscInt p = closure[2 * i]; 2306 PetscInt loc = -1; 2307 2308 PetscCall(PetscFindInt(p, nleaves, ilocal, &loc)); 2309 if (redirect && loc >= 0) { 2310 (*closurePoints)[count].rank = iremote[loc].rank; 2311 (*closurePoints)[count].index = iremote[loc].index; 2312 } else { 2313 (*closurePoints)[count].rank = rank; 2314 (*closurePoints)[count].index = p; 2315 } 2316 } 2317 PetscCall(DMPlexRestoreTransitiveClosure(plex, c, PETSC_TRUE, &closureSize, &closure)); 2318 } 2319 PetscFunctionReturn(PETSC_SUCCESS); 2320 } 2321 2322 static void MPIAPI DMPforestMaxSFNode(void *a, void *b, PetscMPIInt *len, MPI_Datatype *type) 2323 { 2324 PetscMPIInt i; 2325 2326 for (i = 0; i < *len; i++) { 2327 PetscSFNode *A = (PetscSFNode *)a; 2328 PetscSFNode *B = (PetscSFNode *)b; 2329 2330 if (B->rank < 0) *B = *A; 2331 } 2332 } 2333 2334 static PetscErrorCode DMPforestGetTransferSF_Point(DM coarse, DM fine, PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[]) 2335 { 2336 MPI_Comm comm; 2337 PetscMPIInt rank, size; 2338 DM_Forest_pforest *pforestC, *pforestF; 2339 p4est_t *p4estC, *p4estF; 2340 PetscInt numClosureIndices; 2341 PetscInt numClosurePointsC, numClosurePointsF; 2342 PetscSFNode *closurePointsC, *closurePointsF; 2343 p4est_quadrant_t *coverQuads = NULL; 2344 p4est_quadrant_t **treeQuads; 2345 PetscInt *treeQuadCounts; 2346 MPI_Datatype nodeType; 2347 MPI_Datatype nodeClosureType; 2348 MPI_Op sfNodeReduce; 2349 p4est_topidx_t fltF, lltF, t; 2350 DM plexC, plexF; 2351 PetscInt pStartF, pEndF, pStartC, pEndC; 2352 PetscBool saveInCoarse = PETSC_FALSE; 2353 PetscBool saveInFine = PETSC_FALSE; 2354 PetscBool formCids = (childIds != NULL) ? PETSC_TRUE : PETSC_FALSE; 2355 PetscInt *cids = NULL; 2356 2357 PetscFunctionBegin; 2358 pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data; 2359 pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data; 2360 p4estC = pforestC->forest; 2361 p4estF = pforestF->forest; 2362 PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM"); 2363 comm = PetscObjectComm((PetscObject)coarse); 2364 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2365 PetscCallMPI(MPI_Comm_size(comm, &size)); 2366 PetscCall(DMPforestGetPlex(fine, &plexF)); 2367 PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF)); 2368 PetscCall(DMPforestGetPlex(coarse, &plexC)); 2369 PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC)); 2370 { /* check if the results have been cached */ 2371 DM adaptCoarse, adaptFine; 2372 2373 PetscCall(DMForestGetAdaptivityForest(coarse, &adaptCoarse)); 2374 PetscCall(DMForestGetAdaptivityForest(fine, &adaptFine)); 2375 if (adaptCoarse && adaptCoarse->data == fine->data) { /* coarse is adapted from fine */ 2376 if (pforestC->pointSelfToAdaptSF) { 2377 PetscCall(PetscObjectReference((PetscObject)pforestC->pointSelfToAdaptSF)); 2378 *sf = pforestC->pointSelfToAdaptSF; 2379 if (childIds) { 2380 PetscCall(PetscMalloc1(pEndF - pStartF, &cids)); 2381 PetscCall(PetscArraycpy(cids, pforestC->pointSelfToAdaptCids, pEndF - pStartF)); 2382 *childIds = cids; 2383 } 2384 PetscFunctionReturn(PETSC_SUCCESS); 2385 } else { 2386 saveInCoarse = PETSC_TRUE; 2387 formCids = PETSC_TRUE; 2388 } 2389 } else if (adaptFine && adaptFine->data == coarse->data) { /* fine is adapted from coarse */ 2390 if (pforestF->pointAdaptToSelfSF) { 2391 PetscCall(PetscObjectReference((PetscObject)pforestF->pointAdaptToSelfSF)); 2392 *sf = pforestF->pointAdaptToSelfSF; 2393 if (childIds) { 2394 PetscCall(PetscMalloc1(pEndF - pStartF, &cids)); 2395 PetscCall(PetscArraycpy(cids, pforestF->pointAdaptToSelfCids, pEndF - pStartF)); 2396 *childIds = cids; 2397 } 2398 PetscFunctionReturn(PETSC_SUCCESS); 2399 } else { 2400 saveInFine = PETSC_TRUE; 2401 formCids = PETSC_TRUE; 2402 } 2403 } 2404 } 2405 2406 /* count the number of closure points that have dofs and create a list */ 2407 numClosureIndices = P4EST_INSUL; 2408 /* create the datatype */ 2409 PetscCallMPI(MPI_Type_contiguous(2, MPIU_INT, &nodeType)); 2410 PetscCallMPI(MPI_Type_commit(&nodeType)); 2411 PetscCallMPI(MPI_Op_create(DMPforestMaxSFNode, PETSC_FALSE, &sfNodeReduce)); 2412 PetscCallMPI(MPI_Type_contiguous(numClosureIndices * 2, MPIU_INT, &nodeClosureType)); 2413 PetscCallMPI(MPI_Type_commit(&nodeClosureType)); 2414 /* everything has to go through cells: for each cell, create a list of the sfnodes in its closure */ 2415 /* get lists of closure point SF nodes for every cell */ 2416 PetscCall(DMPforestGetCellSFNodes(coarse, numClosureIndices, &numClosurePointsC, &closurePointsC, PETSC_TRUE)); 2417 PetscCall(DMPforestGetCellSFNodes(fine, numClosureIndices, &numClosurePointsF, &closurePointsF, PETSC_FALSE)); 2418 /* create pointers for tree lists */ 2419 fltF = p4estF->first_local_tree; 2420 lltF = p4estF->last_local_tree; 2421 PetscCall(PetscCalloc2(lltF + 1 - fltF, &treeQuads, lltF + 1 - fltF, &treeQuadCounts)); 2422 /* if the partitions don't match, ship the coarse to cover the fine */ 2423 if (size > 1) { 2424 PetscInt p; 2425 2426 for (p = 0; p < size; p++) { 2427 int equal; 2428 2429 PetscCallP4estReturn(equal, p4est_quadrant_is_equal_piggy, (&p4estC->global_first_position[p], &p4estF->global_first_position[p])); 2430 if (!equal) break; 2431 } 2432 if (p < size) { /* non-matching distribution: send the coarse to cover the fine */ 2433 PetscInt cStartC, cEndC; 2434 PetscSF coveringSF; 2435 PetscInt nleaves; 2436 PetscInt count; 2437 PetscSFNode *newClosurePointsC; 2438 p4est_quadrant_t *coverQuadsSend; 2439 p4est_topidx_t fltC = p4estC->first_local_tree; 2440 p4est_topidx_t lltC = p4estC->last_local_tree; 2441 p4est_topidx_t t; 2442 PetscMPIInt blockSizes[4] = {P4EST_DIM, 2, 1, 1}; 2443 MPI_Aint blockOffsets[4] = {offsetof(p4est_quadrant_t, x), offsetof(p4est_quadrant_t, level), offsetof(p4est_quadrant_t, pad16), offsetof(p4est_quadrant_t, p)}; 2444 MPI_Datatype blockTypes[4] = {MPI_INT32_T, MPI_INT8_T, MPI_INT16_T, MPI_INT32_T /* p.which_tree */}; 2445 MPI_Datatype quadStruct, quadType; 2446 2447 PetscCall(DMPlexGetSimplexOrBoxCells(plexC, 0, &cStartC, &cEndC)); 2448 PetscCall(DMPforestGetCellCoveringSF(comm, p4estC, p4estF, pforestC->cLocalStart, pforestC->cLocalEnd, &coveringSF)); 2449 PetscCall(PetscSFGetGraph(coveringSF, NULL, &nleaves, NULL, NULL)); 2450 PetscCall(PetscMalloc1(numClosureIndices * nleaves, &newClosurePointsC)); 2451 PetscCall(PetscMalloc1(nleaves, &coverQuads)); 2452 PetscCall(PetscMalloc1(cEndC - cStartC, &coverQuadsSend)); 2453 count = 0; 2454 for (t = fltC; t <= lltC; t++) { /* unfortunately, we need to pack a send array, since quads are not stored packed in p4est */ 2455 p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]); 2456 PetscInt q; 2457 2458 PetscCall(PetscMemcpy(&coverQuadsSend[count], tree->quadrants.array, tree->quadrants.elem_count * sizeof(p4est_quadrant_t))); 2459 for (q = 0; (size_t)q < tree->quadrants.elem_count; q++) coverQuadsSend[count + q].p.which_tree = t; 2460 count += tree->quadrants.elem_count; 2461 } 2462 /* p is of a union type p4est_quadrant_data, but only the p.which_tree field is active at this time. So, we 2463 have a simple blockTypes[] to use. Note that quadStruct does not count potential padding in array of 2464 p4est_quadrant_t. We have to call MPI_Type_create_resized() to change upper-bound of quadStruct. 2465 */ 2466 PetscCallMPI(MPI_Type_create_struct(4, blockSizes, blockOffsets, blockTypes, &quadStruct)); 2467 PetscCallMPI(MPI_Type_create_resized(quadStruct, 0, sizeof(p4est_quadrant_t), &quadType)); 2468 PetscCallMPI(MPI_Type_commit(&quadType)); 2469 PetscCall(PetscSFBcastBegin(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE)); 2470 PetscCall(PetscSFBcastBegin(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE)); 2471 PetscCall(PetscSFBcastEnd(coveringSF, nodeClosureType, closurePointsC, newClosurePointsC, MPI_REPLACE)); 2472 PetscCall(PetscSFBcastEnd(coveringSF, quadType, coverQuadsSend, coverQuads, MPI_REPLACE)); 2473 PetscCallMPI(MPI_Type_free(&quadStruct)); 2474 PetscCallMPI(MPI_Type_free(&quadType)); 2475 PetscCall(PetscFree(coverQuadsSend)); 2476 PetscCall(PetscFree(closurePointsC)); 2477 PetscCall(PetscSFDestroy(&coveringSF)); 2478 closurePointsC = newClosurePointsC; 2479 2480 /* assign tree quads based on locations in coverQuads */ 2481 { 2482 PetscInt q; 2483 for (q = 0; q < nleaves; q++) { 2484 p4est_locidx_t t = coverQuads[q].p.which_tree; 2485 if (!treeQuadCounts[t - fltF]++) treeQuads[t - fltF] = &coverQuads[q]; 2486 } 2487 } 2488 } 2489 } 2490 if (!coverQuads) { /* matching partitions: assign tree quads based on locations in p4est native arrays */ 2491 for (t = fltF; t <= lltF; t++) { 2492 p4est_tree_t *tree = &(((p4est_tree_t *)p4estC->trees->array)[t]); 2493 2494 treeQuadCounts[t - fltF] = tree->quadrants.elem_count; 2495 treeQuads[t - fltF] = (p4est_quadrant_t *)tree->quadrants.array; 2496 } 2497 } 2498 2499 { 2500 PetscInt p; 2501 PetscInt cLocalStartF; 2502 PetscSF pointSF; 2503 PetscSFNode *roots; 2504 PetscInt *rootType; 2505 DM refTree = NULL; 2506 DMLabel canonical; 2507 PetscInt *childClosures[P4EST_CHILDREN] = {NULL}; 2508 PetscInt *rootClosure = NULL; 2509 PetscInt coarseOffset; 2510 PetscInt numCoarseQuads; 2511 2512 PetscCall(PetscMalloc1(pEndF - pStartF, &roots)); 2513 PetscCall(PetscMalloc1(pEndF - pStartF, &rootType)); 2514 PetscCall(DMGetPointSF(fine, &pointSF)); 2515 for (p = pStartF; p < pEndF; p++) { 2516 roots[p - pStartF].rank = -1; 2517 roots[p - pStartF].index = -1; 2518 rootType[p - pStartF] = -1; 2519 } 2520 if (formCids) { 2521 PetscInt child; 2522 2523 PetscCall(PetscMalloc1(pEndF - pStartF, &cids)); 2524 for (p = pStartF; p < pEndF; p++) cids[p - pStartF] = -2; 2525 PetscCall(DMPlexGetReferenceTree(plexF, &refTree)); 2526 PetscCall(DMPlexGetTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure)); 2527 for (child = 0; child < P4EST_CHILDREN; child++) { /* get the closures of the child cells in the reference tree */ 2528 PetscCall(DMPlexGetTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child])); 2529 } 2530 PetscCall(DMGetLabel(refTree, "canonical", &canonical)); 2531 } 2532 cLocalStartF = pforestF->cLocalStart; 2533 for (t = fltF, coarseOffset = 0, numCoarseQuads = 0; t <= lltF; t++, coarseOffset += numCoarseQuads) { 2534 p4est_tree_t *tree = &(((p4est_tree_t *)p4estF->trees->array)[t]); 2535 PetscInt numFineQuads = tree->quadrants.elem_count; 2536 p4est_quadrant_t *coarseQuads = treeQuads[t - fltF]; 2537 p4est_quadrant_t *fineQuads = (p4est_quadrant_t *)tree->quadrants.array; 2538 PetscInt i, coarseCount = 0; 2539 PetscInt offset = tree->quadrants_offset; 2540 sc_array_t coarseQuadsArray; 2541 2542 numCoarseQuads = treeQuadCounts[t - fltF]; 2543 PetscCallP4est(sc_array_init_data, (&coarseQuadsArray, coarseQuads, sizeof(p4est_quadrant_t), (size_t)numCoarseQuads)); 2544 for (i = 0; i < numFineQuads; i++) { 2545 PetscInt c = i + offset; 2546 p4est_quadrant_t *quad = &fineQuads[i]; 2547 p4est_quadrant_t *quadCoarse = NULL; 2548 ssize_t disjoint = -1; 2549 2550 while (disjoint < 0 && coarseCount < numCoarseQuads) { 2551 quadCoarse = &coarseQuads[coarseCount]; 2552 PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad)); 2553 if (disjoint < 0) coarseCount++; 2554 } 2555 PetscCheck(disjoint == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "did not find overlapping coarse quad"); 2556 if (quadCoarse->level > quad->level || (quadCoarse->level == quad->level && !transferIdent)) { /* the "coarse" mesh is finer than the fine mesh at the point: continue */ 2557 if (transferIdent) { /* find corners */ 2558 PetscInt j = 0; 2559 2560 do { 2561 if (j < P4EST_CHILDREN) { 2562 p4est_quadrant_t cornerQuad; 2563 int equal; 2564 2565 PetscCallP4est(p4est_quadrant_corner_descendant, (quad, &cornerQuad, j, quadCoarse->level)); 2566 PetscCallP4estReturn(equal, p4est_quadrant_is_equal, (&cornerQuad, quadCoarse)); 2567 if (equal) { 2568 PetscInt petscJ = P4estVertToPetscVert[j]; 2569 PetscInt p = closurePointsF[numClosureIndices * c + (P4EST_INSUL - P4EST_CHILDREN) + petscJ].index; 2570 PetscSFNode q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + (P4EST_INSUL - P4EST_CHILDREN) + petscJ]; 2571 2572 roots[p - pStartF] = q; 2573 rootType[p - pStartF] = PETSC_MAX_INT; 2574 cids[p - pStartF] = -1; 2575 j++; 2576 } 2577 } 2578 coarseCount++; 2579 disjoint = 1; 2580 if (coarseCount < numCoarseQuads) { 2581 quadCoarse = &coarseQuads[coarseCount]; 2582 PetscCallP4estReturn(disjoint, p4est_quadrant_disjoint, (quadCoarse, quad)); 2583 } 2584 } while (!disjoint); 2585 } 2586 continue; 2587 } 2588 if (quadCoarse->level == quad->level) { /* same quad present in coarse and fine mesh */ 2589 PetscInt j; 2590 for (j = 0; j < numClosureIndices; j++) { 2591 PetscInt p = closurePointsF[numClosureIndices * c + j].index; 2592 2593 roots[p - pStartF] = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + j]; 2594 rootType[p - pStartF] = PETSC_MAX_INT; /* unconditionally accept */ 2595 cids[p - pStartF] = -1; 2596 } 2597 } else { 2598 PetscInt levelDiff = quad->level - quadCoarse->level; 2599 PetscInt proposedCids[P4EST_INSUL] = {0}; 2600 2601 if (formCids) { 2602 PetscInt cl; 2603 PetscInt *pointClosure = NULL; 2604 int cid; 2605 2606 PetscCheck(levelDiff <= 1, PETSC_COMM_SELF, PETSC_ERR_USER, "Recursive child ids not implemented"); 2607 PetscCallP4estReturn(cid, p4est_quadrant_child_id, (quad)); 2608 PetscCall(DMPlexGetTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure)); 2609 for (cl = 0; cl < P4EST_INSUL; cl++) { 2610 PetscInt p = pointClosure[2 * cl]; 2611 PetscInt point = childClosures[cid][2 * cl]; 2612 PetscInt ornt = childClosures[cid][2 * cl + 1]; 2613 PetscInt newcid = -1; 2614 DMPolytopeType ct; 2615 2616 if (rootType[p - pStartF] == PETSC_MAX_INT) continue; 2617 PetscCall(DMPlexGetCellType(refTree, point, &ct)); 2618 ornt = DMPolytopeConvertNewOrientation_Internal(ct, ornt); 2619 if (!cl) { 2620 newcid = cid + 1; 2621 } else { 2622 PetscInt rcl, parent, parentOrnt = 0; 2623 2624 PetscCall(DMPlexGetTreeParent(refTree, point, &parent, NULL)); 2625 if (parent == point) { 2626 newcid = -1; 2627 } else if (!parent) { /* in the root */ 2628 newcid = point; 2629 } else { 2630 DMPolytopeType rct = DM_POLYTOPE_UNKNOWN; 2631 2632 for (rcl = 1; rcl < P4EST_INSUL; rcl++) { 2633 if (rootClosure[2 * rcl] == parent) { 2634 PetscCall(DMPlexGetCellType(refTree, parent, &rct)); 2635 parentOrnt = DMPolytopeConvertNewOrientation_Internal(rct, rootClosure[2 * rcl + 1]); 2636 break; 2637 } 2638 } 2639 PetscCheck(rcl < P4EST_INSUL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Couldn't find parent in root closure"); 2640 PetscCall(DMPlexReferenceTreeGetChildSymmetry(refTree, parent, parentOrnt, ornt, point, DMPolytopeConvertNewOrientation_Internal(rct, pointClosure[2 * rcl + 1]), NULL, &newcid)); 2641 } 2642 } 2643 if (newcid >= 0) { 2644 if (canonical) PetscCall(DMLabelGetValue(canonical, newcid, &newcid)); 2645 proposedCids[cl] = newcid; 2646 } 2647 } 2648 PetscCall(DMPlexRestoreTransitiveClosure(plexF, c + cLocalStartF, PETSC_TRUE, NULL, &pointClosure)); 2649 } 2650 p4est_qcoord_t coarseBound[2][P4EST_DIM] = { 2651 {quadCoarse->x, quadCoarse->y, 2652 #if defined(P4_TO_P8) 2653 quadCoarse->z 2654 #endif 2655 }, 2656 {0} 2657 }; 2658 p4est_qcoord_t fineBound[2][P4EST_DIM] = { 2659 {quad->x, quad->y, 2660 #if defined(P4_TO_P8) 2661 quad->z 2662 #endif 2663 }, 2664 {0} 2665 }; 2666 PetscInt j; 2667 for (j = 0; j < P4EST_DIM; j++) { /* get the coordinates of cell boundaries in each direction */ 2668 coarseBound[1][j] = coarseBound[0][j] + P4EST_QUADRANT_LEN(quadCoarse->level); 2669 fineBound[1][j] = fineBound[0][j] + P4EST_QUADRANT_LEN(quad->level); 2670 } 2671 for (j = 0; j < numClosureIndices; j++) { 2672 PetscInt l, p; 2673 PetscSFNode q; 2674 2675 p = closurePointsF[numClosureIndices * c + j].index; 2676 if (rootType[p - pStartF] == PETSC_MAX_INT) continue; 2677 if (j == 0) { /* volume: ancestor is volume */ 2678 l = 0; 2679 } else if (j < 1 + P4EST_FACES) { /* facet */ 2680 PetscInt face = PetscFaceToP4estFace[j - 1]; 2681 PetscInt direction = face / 2; 2682 PetscInt coarseFace = -1; 2683 2684 if (coarseBound[face % 2][direction] == fineBound[face % 2][direction]) { 2685 coarseFace = face; 2686 l = 1 + P4estFaceToPetscFace[coarseFace]; 2687 } else { 2688 l = 0; 2689 } 2690 #if defined(P4_TO_P8) 2691 } else if (j < 1 + P4EST_FACES + P8EST_EDGES) { 2692 PetscInt edge = PetscEdgeToP4estEdge[j - (1 + P4EST_FACES)]; 2693 PetscInt direction = edge / 4; 2694 PetscInt mod = edge % 4; 2695 PetscInt coarseEdge = -1, coarseFace = -1; 2696 PetscInt minDir = PetscMin((direction + 1) % 3, (direction + 2) % 3); 2697 PetscInt maxDir = PetscMax((direction + 1) % 3, (direction + 2) % 3); 2698 PetscBool dirTest[2]; 2699 2700 dirTest[0] = (PetscBool)(coarseBound[mod % 2][minDir] == fineBound[mod % 2][minDir]); 2701 dirTest[1] = (PetscBool)(coarseBound[mod / 2][maxDir] == fineBound[mod / 2][maxDir]); 2702 2703 if (dirTest[0] && dirTest[1]) { /* fine edge falls on coarse edge */ 2704 coarseEdge = edge; 2705 l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge]; 2706 } else if (dirTest[0]) { /* fine edge falls on a coarse face in the minDir direction */ 2707 coarseFace = 2 * minDir + (mod % 2); 2708 l = 1 + P4estFaceToPetscFace[coarseFace]; 2709 } else if (dirTest[1]) { /* fine edge falls on a coarse face in the maxDir direction */ 2710 coarseFace = 2 * maxDir + (mod / 2); 2711 l = 1 + P4estFaceToPetscFace[coarseFace]; 2712 } else { 2713 l = 0; 2714 } 2715 #endif 2716 } else { 2717 PetscInt vertex = PetscVertToP4estVert[P4EST_CHILDREN - (P4EST_INSUL - j)]; 2718 PetscBool dirTest[P4EST_DIM]; 2719 PetscInt m; 2720 PetscInt numMatch = 0; 2721 PetscInt coarseVertex = -1, coarseFace = -1; 2722 #if defined(P4_TO_P8) 2723 PetscInt coarseEdge = -1; 2724 #endif 2725 2726 for (m = 0; m < P4EST_DIM; m++) { 2727 dirTest[m] = (PetscBool)(coarseBound[(vertex >> m) & 1][m] == fineBound[(vertex >> m) & 1][m]); 2728 if (dirTest[m]) numMatch++; 2729 } 2730 if (numMatch == P4EST_DIM) { /* vertex on vertex */ 2731 coarseVertex = vertex; 2732 l = P4EST_INSUL - (P4EST_CHILDREN - P4estVertToPetscVert[coarseVertex]); 2733 } else if (numMatch == 1) { /* vertex on face */ 2734 for (m = 0; m < P4EST_DIM; m++) { 2735 if (dirTest[m]) { 2736 coarseFace = 2 * m + ((vertex >> m) & 1); 2737 break; 2738 } 2739 } 2740 l = 1 + P4estFaceToPetscFace[coarseFace]; 2741 #if defined(P4_TO_P8) 2742 } else if (numMatch == 2) { /* vertex on edge */ 2743 for (m = 0; m < P4EST_DIM; m++) { 2744 if (!dirTest[m]) { 2745 PetscInt otherDir1 = (m + 1) % 3; 2746 PetscInt otherDir2 = (m + 2) % 3; 2747 PetscInt minDir = PetscMin(otherDir1, otherDir2); 2748 PetscInt maxDir = PetscMax(otherDir1, otherDir2); 2749 2750 coarseEdge = m * 4 + 2 * ((vertex >> maxDir) & 1) + ((vertex >> minDir) & 1); 2751 break; 2752 } 2753 } 2754 l = 1 + P4EST_FACES + P4estEdgeToPetscEdge[coarseEdge]; 2755 #endif 2756 } else { /* volume */ 2757 l = 0; 2758 } 2759 } 2760 q = closurePointsC[numClosureIndices * (coarseCount + coarseOffset) + l]; 2761 if (l > rootType[p - pStartF]) { 2762 if (l >= P4EST_INSUL - P4EST_CHILDREN) { /* vertex on vertex: unconditional acceptance */ 2763 if (transferIdent) { 2764 roots[p - pStartF] = q; 2765 rootType[p - pStartF] = PETSC_MAX_INT; 2766 if (formCids) cids[p - pStartF] = -1; 2767 } 2768 } else { 2769 PetscInt k, thisp = p, limit; 2770 2771 roots[p - pStartF] = q; 2772 rootType[p - pStartF] = l; 2773 if (formCids) cids[p - pStartF] = proposedCids[j]; 2774 limit = transferIdent ? levelDiff : (levelDiff - 1); 2775 for (k = 0; k < limit; k++) { 2776 PetscInt parent; 2777 2778 PetscCall(DMPlexGetTreeParent(plexF, thisp, &parent, NULL)); 2779 if (parent == thisp) break; 2780 2781 roots[parent - pStartF] = q; 2782 rootType[parent - pStartF] = PETSC_MAX_INT; 2783 if (formCids) cids[parent - pStartF] = -1; 2784 thisp = parent; 2785 } 2786 } 2787 } 2788 } 2789 } 2790 } 2791 } 2792 2793 /* now every cell has labeled the points in its closure, so we first make sure everyone agrees by reducing to roots, and the broadcast the agreements */ 2794 if (size > 1) { 2795 PetscInt *rootTypeCopy, p; 2796 2797 PetscCall(PetscMalloc1(pEndF - pStartF, &rootTypeCopy)); 2798 PetscCall(PetscArraycpy(rootTypeCopy, rootType, pEndF - pStartF)); 2799 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX)); 2800 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_MAX)); 2801 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE)); 2802 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, rootTypeCopy, rootTypeCopy, MPI_REPLACE)); 2803 for (p = pStartF; p < pEndF; p++) { 2804 if (rootTypeCopy[p - pStartF] > rootType[p - pStartF]) { /* another process found a root of higher type (e.g. vertex instead of edge), which we want to accept, so nullify this */ 2805 roots[p - pStartF].rank = -1; 2806 roots[p - pStartF].index = -1; 2807 } 2808 if (formCids && rootTypeCopy[p - pStartF] == PETSC_MAX_INT) { cids[p - pStartF] = -1; /* we have found an antecedent that is the same: no child id */ } 2809 } 2810 PetscCall(PetscFree(rootTypeCopy)); 2811 PetscCall(PetscSFReduceBegin(pointSF, nodeType, roots, roots, sfNodeReduce)); 2812 PetscCall(PetscSFReduceEnd(pointSF, nodeType, roots, roots, sfNodeReduce)); 2813 PetscCall(PetscSFBcastBegin(pointSF, nodeType, roots, roots, MPI_REPLACE)); 2814 PetscCall(PetscSFBcastEnd(pointSF, nodeType, roots, roots, MPI_REPLACE)); 2815 } 2816 PetscCall(PetscFree(rootType)); 2817 2818 { 2819 PetscInt numRoots; 2820 PetscInt numLeaves; 2821 PetscInt *leaves; 2822 PetscSFNode *iremote; 2823 /* count leaves */ 2824 2825 numRoots = pEndC - pStartC; 2826 2827 numLeaves = 0; 2828 for (p = pStartF; p < pEndF; p++) { 2829 if (roots[p - pStartF].index >= 0) numLeaves++; 2830 } 2831 PetscCall(PetscMalloc1(numLeaves, &leaves)); 2832 PetscCall(PetscMalloc1(numLeaves, &iremote)); 2833 numLeaves = 0; 2834 for (p = pStartF; p < pEndF; p++) { 2835 if (roots[p - pStartF].index >= 0) { 2836 leaves[numLeaves] = p - pStartF; 2837 iremote[numLeaves] = roots[p - pStartF]; 2838 numLeaves++; 2839 } 2840 } 2841 PetscCall(PetscFree(roots)); 2842 PetscCall(PetscSFCreate(comm, sf)); 2843 if (numLeaves == (pEndF - pStartF)) { 2844 PetscCall(PetscFree(leaves)); 2845 PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2846 } else { 2847 PetscCall(PetscSFSetGraph(*sf, numRoots, numLeaves, leaves, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2848 } 2849 } 2850 if (formCids) { 2851 PetscSF pointSF; 2852 PetscInt child; 2853 2854 PetscCall(DMPlexGetReferenceTree(plexF, &refTree)); 2855 PetscCall(DMGetPointSF(plexF, &pointSF)); 2856 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, cids, cids, MPI_MAX)); 2857 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, cids, cids, MPI_MAX)); 2858 if (childIds) *childIds = cids; 2859 for (child = 0; child < P4EST_CHILDREN; child++) PetscCall(DMPlexRestoreTransitiveClosure(refTree, child + 1, PETSC_TRUE, NULL, &childClosures[child])); 2860 PetscCall(DMPlexRestoreTransitiveClosure(refTree, 0, PETSC_TRUE, NULL, &rootClosure)); 2861 } 2862 } 2863 if (saveInCoarse) { /* cache results */ 2864 PetscCall(PetscObjectReference((PetscObject)*sf)); 2865 pforestC->pointSelfToAdaptSF = *sf; 2866 if (!childIds) { 2867 pforestC->pointSelfToAdaptCids = cids; 2868 } else { 2869 PetscCall(PetscMalloc1(pEndF - pStartF, &pforestC->pointSelfToAdaptCids)); 2870 PetscCall(PetscArraycpy(pforestC->pointSelfToAdaptCids, cids, pEndF - pStartF)); 2871 } 2872 } else if (saveInFine) { 2873 PetscCall(PetscObjectReference((PetscObject)*sf)); 2874 pforestF->pointAdaptToSelfSF = *sf; 2875 if (!childIds) { 2876 pforestF->pointAdaptToSelfCids = cids; 2877 } else { 2878 PetscCall(PetscMalloc1(pEndF - pStartF, &pforestF->pointAdaptToSelfCids)); 2879 PetscCall(PetscArraycpy(pforestF->pointAdaptToSelfCids, cids, pEndF - pStartF)); 2880 } 2881 } 2882 PetscCall(PetscFree2(treeQuads, treeQuadCounts)); 2883 PetscCall(PetscFree(coverQuads)); 2884 PetscCall(PetscFree(closurePointsC)); 2885 PetscCall(PetscFree(closurePointsF)); 2886 PetscCallMPI(MPI_Type_free(&nodeClosureType)); 2887 PetscCallMPI(MPI_Op_free(&sfNodeReduce)); 2888 PetscCallMPI(MPI_Type_free(&nodeType)); 2889 PetscFunctionReturn(PETSC_SUCCESS); 2890 } 2891 2892 /* children are sf leaves of parents */ 2893 static PetscErrorCode DMPforestGetTransferSF_Internal(DM coarse, DM fine, const PetscInt dofPerDim[], PetscSF *sf, PetscBool transferIdent, PetscInt *childIds[]) 2894 { 2895 MPI_Comm comm; 2896 PetscMPIInt rank; 2897 DM_Forest_pforest *pforestC, *pforestF; 2898 DM plexC, plexF; 2899 PetscInt pStartC, pEndC, pStartF, pEndF; 2900 PetscSF pointTransferSF; 2901 PetscBool allOnes = PETSC_TRUE; 2902 2903 PetscFunctionBegin; 2904 pforestC = (DM_Forest_pforest *)((DM_Forest *)coarse->data)->data; 2905 pforestF = (DM_Forest_pforest *)((DM_Forest *)fine->data)->data; 2906 PetscCheck(pforestC->topo == pforestF->topo, PetscObjectComm((PetscObject)coarse), PETSC_ERR_ARG_INCOMP, "DM's must have the same base DM"); 2907 comm = PetscObjectComm((PetscObject)coarse); 2908 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2909 2910 { 2911 PetscInt i; 2912 for (i = 0; i <= P4EST_DIM; i++) { 2913 if (dofPerDim[i] != 1) { 2914 allOnes = PETSC_FALSE; 2915 break; 2916 } 2917 } 2918 } 2919 PetscCall(DMPforestGetTransferSF_Point(coarse, fine, &pointTransferSF, transferIdent, childIds)); 2920 if (allOnes) { 2921 *sf = pointTransferSF; 2922 PetscFunctionReturn(PETSC_SUCCESS); 2923 } 2924 2925 PetscCall(DMPforestGetPlex(fine, &plexF)); 2926 PetscCall(DMPlexGetChart(plexF, &pStartF, &pEndF)); 2927 PetscCall(DMPforestGetPlex(coarse, &plexC)); 2928 PetscCall(DMPlexGetChart(plexC, &pStartC, &pEndC)); 2929 { 2930 PetscInt numRoots; 2931 PetscInt numLeaves; 2932 const PetscInt *leaves; 2933 const PetscSFNode *iremote; 2934 PetscInt d; 2935 PetscSection leafSection, rootSection; 2936 2937 /* count leaves */ 2938 PetscCall(PetscSFGetGraph(pointTransferSF, &numRoots, &numLeaves, &leaves, &iremote)); 2939 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &rootSection)); 2940 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &leafSection)); 2941 PetscCall(PetscSectionSetChart(rootSection, pStartC, pEndC)); 2942 PetscCall(PetscSectionSetChart(leafSection, pStartF, pEndF)); 2943 2944 for (d = 0; d <= P4EST_DIM; d++) { 2945 PetscInt startC, endC, e; 2946 2947 PetscCall(DMPlexGetSimplexOrBoxCells(plexC, P4EST_DIM - d, &startC, &endC)); 2948 for (e = startC; e < endC; e++) PetscCall(PetscSectionSetDof(rootSection, e, dofPerDim[d])); 2949 } 2950 2951 for (d = 0; d <= P4EST_DIM; d++) { 2952 PetscInt startF, endF, e; 2953 2954 PetscCall(DMPlexGetSimplexOrBoxCells(plexF, P4EST_DIM - d, &startF, &endF)); 2955 for (e = startF; e < endF; e++) PetscCall(PetscSectionSetDof(leafSection, e, dofPerDim[d])); 2956 } 2957 2958 PetscCall(PetscSectionSetUp(rootSection)); 2959 PetscCall(PetscSectionSetUp(leafSection)); 2960 { 2961 PetscInt nroots, nleaves; 2962 PetscInt *mine, i, p; 2963 PetscInt *offsets, *offsetsRoot; 2964 PetscSFNode *remote; 2965 2966 PetscCall(PetscMalloc1(pEndF - pStartF, &offsets)); 2967 PetscCall(PetscMalloc1(pEndC - pStartC, &offsetsRoot)); 2968 for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionGetOffset(rootSection, p, &offsetsRoot[p - pStartC])); 2969 PetscCall(PetscSFBcastBegin(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE)); 2970 PetscCall(PetscSFBcastEnd(pointTransferSF, MPIU_INT, offsetsRoot, offsets, MPI_REPLACE)); 2971 PetscCall(PetscSectionGetStorageSize(rootSection, &nroots)); 2972 nleaves = 0; 2973 for (i = 0; i < numLeaves; i++) { 2974 PetscInt leaf = leaves ? leaves[i] : i; 2975 PetscInt dof; 2976 2977 PetscCall(PetscSectionGetDof(leafSection, leaf, &dof)); 2978 nleaves += dof; 2979 } 2980 PetscCall(PetscMalloc1(nleaves, &mine)); 2981 PetscCall(PetscMalloc1(nleaves, &remote)); 2982 nleaves = 0; 2983 for (i = 0; i < numLeaves; i++) { 2984 PetscInt leaf = leaves ? leaves[i] : i; 2985 PetscInt dof; 2986 PetscInt off, j; 2987 2988 PetscCall(PetscSectionGetDof(leafSection, leaf, &dof)); 2989 PetscCall(PetscSectionGetOffset(leafSection, leaf, &off)); 2990 for (j = 0; j < dof; j++) { 2991 remote[nleaves].rank = iremote[i].rank; 2992 remote[nleaves].index = offsets[leaf] + j; 2993 mine[nleaves++] = off + j; 2994 } 2995 } 2996 PetscCall(PetscFree(offsetsRoot)); 2997 PetscCall(PetscFree(offsets)); 2998 PetscCall(PetscSFCreate(comm, sf)); 2999 PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, mine, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER)); 3000 } 3001 PetscCall(PetscSectionDestroy(&leafSection)); 3002 PetscCall(PetscSectionDestroy(&rootSection)); 3003 PetscCall(PetscSFDestroy(&pointTransferSF)); 3004 } 3005 PetscFunctionReturn(PETSC_SUCCESS); 3006 } 3007 3008 static PetscErrorCode DMPforestGetTransferSF(DM dmA, DM dmB, const PetscInt dofPerDim[], PetscSF *sfAtoB, PetscSF *sfBtoA) 3009 { 3010 DM adaptA, adaptB; 3011 DMAdaptFlag purpose; 3012 3013 PetscFunctionBegin; 3014 PetscCall(DMForestGetAdaptivityForest(dmA, &adaptA)); 3015 PetscCall(DMForestGetAdaptivityForest(dmB, &adaptB)); 3016 /* it is more efficient when the coarser mesh is the first argument: reorder if we know one is coarser than the other */ 3017 if (adaptA && adaptA->data == dmB->data) { /* dmA was adapted from dmB */ 3018 PetscCall(DMForestGetAdaptivityPurpose(dmA, &purpose)); 3019 if (purpose == DM_ADAPT_REFINE) { 3020 PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB)); 3021 PetscFunctionReturn(PETSC_SUCCESS); 3022 } 3023 } else if (adaptB && adaptB->data == dmA->data) { /* dmB was adapted from dmA */ 3024 PetscCall(DMForestGetAdaptivityPurpose(dmB, &purpose)); 3025 if (purpose == DM_ADAPT_COARSEN) { 3026 PetscCall(DMPforestGetTransferSF(dmB, dmA, dofPerDim, sfBtoA, sfAtoB)); 3027 PetscFunctionReturn(PETSC_SUCCESS); 3028 } 3029 } 3030 if (sfAtoB) PetscCall(DMPforestGetTransferSF_Internal(dmA, dmB, dofPerDim, sfAtoB, PETSC_TRUE, NULL)); 3031 if (sfBtoA) PetscCall(DMPforestGetTransferSF_Internal(dmB, dmA, dofPerDim, sfBtoA, (PetscBool)(sfAtoB == NULL), NULL)); 3032 PetscFunctionReturn(PETSC_SUCCESS); 3033 } 3034 3035 static PetscErrorCode DMPforestLabelsInitialize(DM dm, DM plex) 3036 { 3037 DM_Forest *forest = (DM_Forest *)dm->data; 3038 DM_Forest_pforest *pforest = (DM_Forest_pforest *)forest->data; 3039 PetscInt cLocalStart, cLocalEnd, cStart, cEnd, fStart, fEnd, eStart, eEnd, vStart, vEnd; 3040 PetscInt cStartBase, cEndBase, fStartBase, fEndBase, vStartBase, vEndBase, eStartBase, eEndBase; 3041 PetscInt pStart, pEnd, pStartBase, pEndBase, p; 3042 DM base; 3043 PetscInt *star = NULL, starSize; 3044 DMLabelLink next = dm->labels; 3045 PetscInt guess = 0; 3046 p4est_topidx_t num_trees = pforest->topo->conn->num_trees; 3047 3048 PetscFunctionBegin; 3049 pforest->labelsFinalized = PETSC_TRUE; 3050 cLocalStart = pforest->cLocalStart; 3051 cLocalEnd = pforest->cLocalEnd; 3052 PetscCall(DMForestGetBaseDM(dm, &base)); 3053 if (!base) { 3054 if (pforest->ghostName) { /* insert a label to make the boundaries, with stratum values denoting which face of the element touches the boundary */ 3055 p4est_connectivity_t *conn = pforest->topo->conn; 3056 p4est_t *p4est = pforest->forest; 3057 p4est_tree_t *trees = (p4est_tree_t *)p4est->trees->array; 3058 p4est_topidx_t t, flt = p4est->first_local_tree; 3059 p4est_topidx_t llt = pforest->forest->last_local_tree; 3060 DMLabel ghostLabel; 3061 PetscInt c; 3062 3063 PetscCall(DMCreateLabel(plex, pforest->ghostName)); 3064 PetscCall(DMGetLabel(plex, pforest->ghostName, &ghostLabel)); 3065 for (c = cLocalStart, t = flt; t <= llt; t++) { 3066 p4est_tree_t *tree = &trees[t]; 3067 p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array; 3068 PetscInt numQuads = (PetscInt)tree->quadrants.elem_count; 3069 PetscInt q; 3070 3071 for (q = 0; q < numQuads; q++, c++) { 3072 p4est_quadrant_t *quad = &quads[q]; 3073 PetscInt f; 3074 3075 for (f = 0; f < P4EST_FACES; f++) { 3076 p4est_quadrant_t neigh; 3077 int isOutside; 3078 3079 PetscCallP4est(p4est_quadrant_face_neighbor, (quad, f, &neigh)); 3080 PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&neigh)); 3081 if (isOutside) { 3082 p4est_topidx_t nt; 3083 PetscInt nf; 3084 3085 nt = conn->tree_to_tree[t * P4EST_FACES + f]; 3086 nf = (PetscInt)conn->tree_to_face[t * P4EST_FACES + f]; 3087 nf = nf % P4EST_FACES; 3088 if (nt == t && nf == f) { 3089 PetscInt plexF = P4estFaceToPetscFace[f]; 3090 const PetscInt *cone; 3091 3092 PetscCall(DMPlexGetCone(plex, c, &cone)); 3093 PetscCall(DMLabelSetValue(ghostLabel, cone[plexF], plexF + 1)); 3094 } 3095 } 3096 } 3097 } 3098 } 3099 } 3100 PetscFunctionReturn(PETSC_SUCCESS); 3101 } 3102 PetscCall(DMPlexGetSimplexOrBoxCells(base, 0, &cStartBase, &cEndBase)); 3103 PetscCall(DMPlexGetSimplexOrBoxCells(base, 1, &fStartBase, &fEndBase)); 3104 PetscCall(DMPlexGetSimplexOrBoxCells(base, P4EST_DIM - 1, &eStartBase, &eEndBase)); 3105 PetscCall(DMPlexGetDepthStratum(base, 0, &vStartBase, &vEndBase)); 3106 3107 PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd)); 3108 PetscCall(DMPlexGetSimplexOrBoxCells(plex, 1, &fStart, &fEnd)); 3109 PetscCall(DMPlexGetSimplexOrBoxCells(plex, P4EST_DIM - 1, &eStart, &eEnd)); 3110 PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd)); 3111 3112 PetscCall(DMPlexGetChart(plex, &pStart, &pEnd)); 3113 PetscCall(DMPlexGetChart(base, &pStartBase, &pEndBase)); 3114 /* go through the mesh: use star to find a quadrant that borders a point. Use the closure to determine the 3115 * orientation of the quadrant relative to that point. Use that to relate the point to the numbering in the base 3116 * mesh, and extract a label value (since the base mesh is redundantly distributed, can be found locally). */ 3117 while (next) { 3118 DMLabel baseLabel; 3119 DMLabel label = next->label; 3120 PetscBool isDepth, isCellType, isGhost, isVTK, isSpmap; 3121 const char *name; 3122 3123 PetscCall(PetscObjectGetName((PetscObject)label, &name)); 3124 PetscCall(PetscStrcmp(name, "depth", &isDepth)); 3125 if (isDepth) { 3126 next = next->next; 3127 continue; 3128 } 3129 PetscCall(PetscStrcmp(name, "celltype", &isCellType)); 3130 if (isCellType) { 3131 next = next->next; 3132 continue; 3133 } 3134 PetscCall(PetscStrcmp(name, "ghost", &isGhost)); 3135 if (isGhost) { 3136 next = next->next; 3137 continue; 3138 } 3139 PetscCall(PetscStrcmp(name, "vtk", &isVTK)); 3140 if (isVTK) { 3141 next = next->next; 3142 continue; 3143 } 3144 PetscCall(PetscStrcmp(name, "_forest_base_subpoint_map", &isSpmap)); 3145 if (!isSpmap) { 3146 PetscCall(DMGetLabel(base, name, &baseLabel)); 3147 if (!baseLabel) { 3148 next = next->next; 3149 continue; 3150 } 3151 PetscCall(DMLabelCreateIndex(baseLabel, pStartBase, pEndBase)); 3152 } else baseLabel = NULL; 3153 3154 for (p = pStart; p < pEnd; p++) { 3155 PetscInt s, c = -1, l; 3156 PetscInt *closure = NULL, closureSize; 3157 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3158 p4est_tree_t *trees = (p4est_tree_t *)pforest->forest->trees->array; 3159 p4est_quadrant_t *q; 3160 PetscInt t, val; 3161 PetscBool zerosupportpoint = PETSC_FALSE; 3162 3163 PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star)); 3164 for (s = 0; s < starSize; s++) { 3165 PetscInt point = star[2 * s]; 3166 3167 if (cStart <= point && point < cEnd) { 3168 PetscCall(DMPlexGetTransitiveClosure(plex, point, PETSC_TRUE, &closureSize, &closure)); 3169 for (l = 0; l < closureSize; l++) { 3170 PetscInt qParent = closure[2 * l], q, pp = p, pParent = p; 3171 do { /* check parents of q */ 3172 q = qParent; 3173 if (q == p) { 3174 c = point; 3175 break; 3176 } 3177 PetscCall(DMPlexGetTreeParent(plex, q, &qParent, NULL)); 3178 } while (qParent != q); 3179 if (c != -1) break; 3180 PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL)); 3181 q = closure[2 * l]; 3182 while (pParent != pp) { /* check parents of p */ 3183 pp = pParent; 3184 if (pp == q) { 3185 c = point; 3186 break; 3187 } 3188 PetscCall(DMPlexGetTreeParent(plex, pp, &pParent, NULL)); 3189 } 3190 if (c != -1) break; 3191 } 3192 PetscCall(DMPlexRestoreTransitiveClosure(plex, point, PETSC_TRUE, NULL, &closure)); 3193 if (l < closureSize) break; 3194 } else { 3195 PetscInt supportSize; 3196 3197 PetscCall(DMPlexGetSupportSize(plex, point, &supportSize)); 3198 zerosupportpoint = (PetscBool)(zerosupportpoint || !supportSize); 3199 } 3200 } 3201 if (c < 0) { 3202 const char *prefix; 3203 PetscBool print = PETSC_FALSE; 3204 3205 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 3206 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, prefix, "-dm_forest_print_label_error", &print, NULL)); 3207 if (print) { 3208 PetscInt i; 3209 3210 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Failed to find cell with point %" PetscInt_FMT " in its closure for label %s (starSize %" PetscInt_FMT ")\n", PetscGlobalRank, p, baseLabel ? ((PetscObject)baseLabel)->name : "_forest_base_subpoint_map", starSize)); 3211 for (i = 0; i < starSize; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, " star[%" PetscInt_FMT "] = %" PetscInt_FMT ",%" PetscInt_FMT "\n", i, star[2 * i], star[2 * i + 1])); 3212 } 3213 PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star)); 3214 if (zerosupportpoint) continue; 3215 else 3216 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Failed to find cell with point %" PetscInt_FMT " in its closure for label %s. Rerun with -dm_forest_print_label_error for more information", p, baseLabel ? ((PetscObject)baseLabel)->name : "_forest_base_subpoint_map"); 3217 } 3218 PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, NULL, &star)); 3219 3220 if (c < cLocalStart) { 3221 /* get from the beginning of the ghost layer */ 3222 q = &ghosts[c]; 3223 t = (PetscInt)q->p.which_tree; 3224 } else if (c < cLocalEnd) { 3225 PetscInt lo = 0, hi = num_trees; 3226 /* get from local quadrants: have to find the right tree */ 3227 3228 c -= cLocalStart; 3229 3230 do { 3231 p4est_tree_t *tree; 3232 3233 PetscCheck(guess >= lo && guess < num_trees && lo < hi, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed binary search"); 3234 tree = &trees[guess]; 3235 if (c < tree->quadrants_offset) { 3236 hi = guess; 3237 } else if (c < tree->quadrants_offset + (PetscInt)tree->quadrants.elem_count) { 3238 q = &((p4est_quadrant_t *)tree->quadrants.array)[c - (PetscInt)tree->quadrants_offset]; 3239 t = guess; 3240 break; 3241 } else { 3242 lo = guess + 1; 3243 } 3244 guess = lo + (hi - lo) / 2; 3245 } while (1); 3246 } else { 3247 /* get from the end of the ghost layer */ 3248 c -= (cLocalEnd - cLocalStart); 3249 3250 q = &ghosts[c]; 3251 t = (PetscInt)q->p.which_tree; 3252 } 3253 3254 if (l == 0) { /* cell */ 3255 if (baseLabel) { 3256 PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val)); 3257 } else { 3258 val = t + cStartBase; 3259 } 3260 PetscCall(DMLabelSetValue(label, p, val)); 3261 } else if (l >= 1 && l < 1 + P4EST_FACES) { /* facet */ 3262 p4est_quadrant_t nq; 3263 int isInside; 3264 3265 l = PetscFaceToP4estFace[l - 1]; 3266 PetscCallP4est(p4est_quadrant_face_neighbor, (q, l, &nq)); 3267 PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq)); 3268 if (isInside) { 3269 /* this facet is in the interior of a tree, so it inherits the label of the tree */ 3270 if (baseLabel) { 3271 PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val)); 3272 } else { 3273 val = t + cStartBase; 3274 } 3275 PetscCall(DMLabelSetValue(label, p, val)); 3276 } else { 3277 PetscInt f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + l]; 3278 3279 if (baseLabel) { 3280 PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val)); 3281 } else { 3282 val = f + fStartBase; 3283 } 3284 PetscCall(DMLabelSetValue(label, p, val)); 3285 } 3286 #if defined(P4_TO_P8) 3287 } else if (l >= 1 + P4EST_FACES && l < 1 + P4EST_FACES + P8EST_EDGES) { /* edge */ 3288 p4est_quadrant_t nq; 3289 int isInside; 3290 3291 l = PetscEdgeToP4estEdge[l - (1 + P4EST_FACES)]; 3292 PetscCallP4est(p8est_quadrant_edge_neighbor, (q, l, &nq)); 3293 PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq)); 3294 if (isInside) { 3295 /* this edge is in the interior of a tree, so it inherits the label of the tree */ 3296 if (baseLabel) { 3297 PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val)); 3298 } else { 3299 val = t + cStartBase; 3300 } 3301 PetscCall(DMLabelSetValue(label, p, val)); 3302 } else { 3303 int isOutsideFace; 3304 3305 PetscCallP4estReturn(isOutsideFace, p4est_quadrant_is_outside_face, (&nq)); 3306 if (isOutsideFace) { 3307 PetscInt f; 3308 3309 if (nq.x < 0) { 3310 f = 0; 3311 } else if (nq.x >= P4EST_ROOT_LEN) { 3312 f = 1; 3313 } else if (nq.y < 0) { 3314 f = 2; 3315 } else if (nq.y >= P4EST_ROOT_LEN) { 3316 f = 3; 3317 } else if (nq.z < 0) { 3318 f = 4; 3319 } else { 3320 f = 5; 3321 } 3322 f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f]; 3323 if (baseLabel) { 3324 PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val)); 3325 } else { 3326 val = f + fStartBase; 3327 } 3328 PetscCall(DMLabelSetValue(label, p, val)); 3329 } else { /* the quadrant edge corresponds to the tree edge */ 3330 PetscInt e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + l]; 3331 3332 if (baseLabel) { 3333 PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val)); 3334 } else { 3335 val = e + eStartBase; 3336 } 3337 PetscCall(DMLabelSetValue(label, p, val)); 3338 } 3339 } 3340 #endif 3341 } else { /* vertex */ 3342 p4est_quadrant_t nq; 3343 int isInside; 3344 3345 #if defined(P4_TO_P8) 3346 l = PetscVertToP4estVert[l - (1 + P4EST_FACES + P8EST_EDGES)]; 3347 #else 3348 l = PetscVertToP4estVert[l - (1 + P4EST_FACES)]; 3349 #endif 3350 PetscCallP4est(p4est_quadrant_corner_neighbor, (q, l, &nq)); 3351 PetscCallP4estReturn(isInside, p4est_quadrant_is_inside_root, (&nq)); 3352 if (isInside) { 3353 if (baseLabel) { 3354 PetscCall(DMLabelGetValue(baseLabel, t + cStartBase, &val)); 3355 } else { 3356 val = t + cStartBase; 3357 } 3358 PetscCall(DMLabelSetValue(label, p, val)); 3359 } else { 3360 int isOutside; 3361 3362 PetscCallP4estReturn(isOutside, p4est_quadrant_is_outside_face, (&nq)); 3363 if (isOutside) { 3364 PetscInt f = -1; 3365 3366 if (nq.x < 0) { 3367 f = 0; 3368 } else if (nq.x >= P4EST_ROOT_LEN) { 3369 f = 1; 3370 } else if (nq.y < 0) { 3371 f = 2; 3372 } else if (nq.y >= P4EST_ROOT_LEN) { 3373 f = 3; 3374 #if defined(P4_TO_P8) 3375 } else if (nq.z < 0) { 3376 f = 4; 3377 } else { 3378 f = 5; 3379 #endif 3380 } 3381 f = pforest->topo->tree_face_to_uniq[P4EST_FACES * t + f]; 3382 if (baseLabel) { 3383 PetscCall(DMLabelGetValue(baseLabel, f + fStartBase, &val)); 3384 } else { 3385 val = f + fStartBase; 3386 } 3387 PetscCall(DMLabelSetValue(label, p, val)); 3388 continue; 3389 } 3390 #if defined(P4_TO_P8) 3391 PetscCallP4estReturn(isOutside, p8est_quadrant_is_outside_edge, (&nq)); 3392 if (isOutside) { 3393 /* outside edge */ 3394 PetscInt e = -1; 3395 3396 if (nq.x >= 0 && nq.x < P4EST_ROOT_LEN) { 3397 if (nq.z < 0) { 3398 if (nq.y < 0) { 3399 e = 0; 3400 } else { 3401 e = 1; 3402 } 3403 } else { 3404 if (nq.y < 0) { 3405 e = 2; 3406 } else { 3407 e = 3; 3408 } 3409 } 3410 } else if (nq.y >= 0 && nq.y < P4EST_ROOT_LEN) { 3411 if (nq.z < 0) { 3412 if (nq.x < 0) { 3413 e = 4; 3414 } else { 3415 e = 5; 3416 } 3417 } else { 3418 if (nq.x < 0) { 3419 e = 6; 3420 } else { 3421 e = 7; 3422 } 3423 } 3424 } else { 3425 if (nq.y < 0) { 3426 if (nq.x < 0) { 3427 e = 8; 3428 } else { 3429 e = 9; 3430 } 3431 } else { 3432 if (nq.x < 0) { 3433 e = 10; 3434 } else { 3435 e = 11; 3436 } 3437 } 3438 } 3439 3440 e = pforest->topo->conn->tree_to_edge[P8EST_EDGES * t + e]; 3441 if (baseLabel) { 3442 PetscCall(DMLabelGetValue(baseLabel, e + eStartBase, &val)); 3443 } else { 3444 val = e + eStartBase; 3445 } 3446 PetscCall(DMLabelSetValue(label, p, val)); 3447 continue; 3448 } 3449 #endif 3450 { 3451 /* outside vertex: same corner as quadrant corner */ 3452 PetscInt v = pforest->topo->conn->tree_to_corner[P4EST_CHILDREN * t + l]; 3453 3454 if (baseLabel) { 3455 PetscCall(DMLabelGetValue(baseLabel, v + vStartBase, &val)); 3456 } else { 3457 val = v + vStartBase; 3458 } 3459 PetscCall(DMLabelSetValue(label, p, val)); 3460 } 3461 } 3462 } 3463 } 3464 next = next->next; 3465 } 3466 PetscFunctionReturn(PETSC_SUCCESS); 3467 } 3468 3469 static PetscErrorCode DMPforestLabelsFinalize(DM dm, DM plex) 3470 { 3471 DM_Forest_pforest *pforest = (DM_Forest_pforest *)((DM_Forest *)dm->data)->data; 3472 DM adapt; 3473 3474 PetscFunctionBegin; 3475 if (pforest->labelsFinalized) PetscFunctionReturn(PETSC_SUCCESS); 3476 pforest->labelsFinalized = PETSC_TRUE; 3477 PetscCall(DMForestGetAdaptivityForest(dm, &adapt)); 3478 if (!adapt) { 3479 /* Initialize labels from the base dm */ 3480 PetscCall(DMPforestLabelsInitialize(dm, plex)); 3481 } else { 3482 PetscInt dofPerDim[4] = {1, 1, 1, 1}; 3483 PetscSF transferForward, transferBackward, pointSF; 3484 PetscInt pStart, pEnd, pStartA, pEndA; 3485 PetscInt *values, *adaptValues; 3486 DMLabelLink next = adapt->labels; 3487 DMLabel adaptLabel; 3488 DM adaptPlex; 3489 3490 PetscCall(DMForestGetAdaptivityLabel(dm, &adaptLabel)); 3491 PetscCall(DMPforestGetPlex(adapt, &adaptPlex)); 3492 PetscCall(DMPforestGetTransferSF(adapt, dm, dofPerDim, &transferForward, &transferBackward)); 3493 PetscCall(DMPlexGetChart(plex, &pStart, &pEnd)); 3494 PetscCall(DMPlexGetChart(adaptPlex, &pStartA, &pEndA)); 3495 PetscCall(PetscMalloc2(pEnd - pStart, &values, pEndA - pStartA, &adaptValues)); 3496 PetscCall(DMGetPointSF(plex, &pointSF)); 3497 if (PetscDefined(USE_DEBUG)) { 3498 PetscInt p; 3499 for (p = pStartA; p < pEndA; p++) adaptValues[p - pStartA] = -1; 3500 for (p = pStart; p < pEnd; p++) values[p - pStart] = -2; 3501 if (transferForward) { 3502 PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE)); 3503 PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE)); 3504 } 3505 if (transferBackward) { 3506 PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX)); 3507 PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX)); 3508 } 3509 for (p = pStart; p < pEnd; p++) { 3510 PetscInt q = p, parent; 3511 3512 PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL)); 3513 while (parent != q) { 3514 if (values[parent] == -2) values[parent] = values[q]; 3515 q = parent; 3516 PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL)); 3517 } 3518 } 3519 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX)); 3520 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX)); 3521 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE)); 3522 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE)); 3523 for (p = pStart; p < pEnd; p++) PetscCheck(values[p - pStart] != -2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "uncovered point %" PetscInt_FMT, p); 3524 } 3525 while (next) { 3526 DMLabel nextLabel = next->label; 3527 const char *name; 3528 PetscBool isDepth, isCellType, isGhost, isVTK; 3529 DMLabel label; 3530 PetscInt p; 3531 3532 PetscCall(PetscObjectGetName((PetscObject)nextLabel, &name)); 3533 PetscCall(PetscStrcmp(name, "depth", &isDepth)); 3534 if (isDepth) { 3535 next = next->next; 3536 continue; 3537 } 3538 PetscCall(PetscStrcmp(name, "celltype", &isCellType)); 3539 if (isCellType) { 3540 next = next->next; 3541 continue; 3542 } 3543 PetscCall(PetscStrcmp(name, "ghost", &isGhost)); 3544 if (isGhost) { 3545 next = next->next; 3546 continue; 3547 } 3548 PetscCall(PetscStrcmp(name, "vtk", &isVTK)); 3549 if (isVTK) { 3550 next = next->next; 3551 continue; 3552 } 3553 if (nextLabel == adaptLabel) { 3554 next = next->next; 3555 continue; 3556 } 3557 /* label was created earlier */ 3558 PetscCall(DMGetLabel(dm, name, &label)); 3559 for (p = pStartA; p < pEndA; p++) PetscCall(DMLabelGetValue(nextLabel, p, &adaptValues[p])); 3560 for (p = pStart; p < pEnd; p++) values[p] = PETSC_MIN_INT; 3561 3562 if (transferForward) PetscCall(PetscSFBcastBegin(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE)); 3563 if (transferBackward) PetscCall(PetscSFReduceBegin(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX)); 3564 if (transferForward) PetscCall(PetscSFBcastEnd(transferForward, MPIU_INT, adaptValues, values, MPI_REPLACE)); 3565 if (transferBackward) PetscCall(PetscSFReduceEnd(transferBackward, MPIU_INT, adaptValues, values, MPI_MAX)); 3566 for (p = pStart; p < pEnd; p++) { 3567 PetscInt q = p, parent; 3568 3569 PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL)); 3570 while (parent != q) { 3571 if (values[parent] == PETSC_MIN_INT) values[parent] = values[q]; 3572 q = parent; 3573 PetscCall(DMPlexGetTreeParent(plex, q, &parent, NULL)); 3574 } 3575 } 3576 PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, values, values, MPI_MAX)); 3577 PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, values, values, MPI_MAX)); 3578 PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, values, values, MPI_REPLACE)); 3579 PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, values, values, MPI_REPLACE)); 3580 3581 for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(label, p, values[p])); 3582 next = next->next; 3583 } 3584 PetscCall(PetscFree2(values, adaptValues)); 3585 PetscCall(PetscSFDestroy(&transferForward)); 3586 PetscCall(PetscSFDestroy(&transferBackward)); 3587 pforest->labelsFinalized = PETSC_TRUE; 3588 } 3589 PetscFunctionReturn(PETSC_SUCCESS); 3590 } 3591 3592 static PetscErrorCode DMPforestMapCoordinates_Cell(DM plex, p4est_geometry_t *geom, PetscInt cell, p4est_quadrant_t *q, p4est_topidx_t t, p4est_connectivity_t *conn, PetscScalar *coords) 3593 { 3594 PetscInt closureSize, c, coordStart, coordEnd, coordDim; 3595 PetscInt *closure = NULL; 3596 PetscSection coordSec; 3597 3598 PetscFunctionBegin; 3599 PetscCall(DMGetCoordinateSection(plex, &coordSec)); 3600 PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd)); 3601 PetscCall(DMGetCoordinateDim(plex, &coordDim)); 3602 PetscCall(DMPlexGetTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure)); 3603 for (c = 0; c < closureSize; c++) { 3604 PetscInt point = closure[2 * c]; 3605 3606 if (point >= coordStart && point < coordEnd) { 3607 PetscInt dof, off; 3608 PetscInt nCoords, i; 3609 PetscCall(PetscSectionGetDof(coordSec, point, &dof)); 3610 PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout"); 3611 nCoords = dof / coordDim; 3612 PetscCall(PetscSectionGetOffset(coordSec, point, &off)); 3613 for (i = 0; i < nCoords; i++) { 3614 PetscScalar *coord = &coords[off + i * coordDim]; 3615 double coordP4est[3] = {0.}; 3616 double coordP4estMapped[3] = {0.}; 3617 PetscInt j; 3618 PetscReal treeCoords[P4EST_CHILDREN][3] = {{0.}}; 3619 PetscReal eta[3] = {0.}; 3620 PetscInt numRounds = 10; 3621 PetscReal coordGuess[3] = {0.}; 3622 3623 eta[0] = (PetscReal)q->x / (PetscReal)P4EST_ROOT_LEN; 3624 eta[1] = (PetscReal)q->y / (PetscReal)P4EST_ROOT_LEN; 3625 #if defined(P4_TO_P8) 3626 eta[2] = (PetscReal)q->z / (PetscReal)P4EST_ROOT_LEN; 3627 #endif 3628 3629 for (j = 0; j < P4EST_CHILDREN; j++) { 3630 PetscInt k; 3631 3632 for (k = 0; k < 3; k++) treeCoords[j][k] = conn->vertices[3 * conn->tree_to_vertex[P4EST_CHILDREN * t + j] + k]; 3633 } 3634 3635 for (j = 0; j < P4EST_CHILDREN; j++) { 3636 PetscInt k; 3637 PetscReal prod = 1.; 3638 3639 for (k = 0; k < P4EST_DIM; k++) prod *= (j & (1 << k)) ? eta[k] : (1. - eta[k]); 3640 for (k = 0; k < 3; k++) coordGuess[k] += prod * treeCoords[j][k]; 3641 } 3642 3643 for (j = 0; j < numRounds; j++) { 3644 PetscInt dir; 3645 3646 for (dir = 0; dir < P4EST_DIM; dir++) { 3647 PetscInt k; 3648 PetscReal diff[3]; 3649 PetscReal dXdeta[3] = {0.}; 3650 PetscReal rhs, scale, update; 3651 3652 for (k = 0; k < 3; k++) diff[k] = coordP4est[k] - coordGuess[k]; 3653 for (k = 0; k < P4EST_CHILDREN; k++) { 3654 PetscInt l; 3655 PetscReal prod = 1.; 3656 3657 for (l = 0; l < P4EST_DIM; l++) { 3658 if (l == dir) { 3659 prod *= (k & (1 << l)) ? 1. : -1.; 3660 } else { 3661 prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]); 3662 } 3663 } 3664 for (l = 0; l < 3; l++) dXdeta[l] += prod * treeCoords[k][l]; 3665 } 3666 rhs = 0.; 3667 scale = 0; 3668 for (k = 0; k < 3; k++) { 3669 rhs += diff[k] * dXdeta[k]; 3670 scale += dXdeta[k] * dXdeta[k]; 3671 } 3672 update = rhs / scale; 3673 eta[dir] += update; 3674 eta[dir] = PetscMin(eta[dir], 1.); 3675 eta[dir] = PetscMax(eta[dir], 0.); 3676 3677 coordGuess[0] = coordGuess[1] = coordGuess[2] = 0.; 3678 for (k = 0; k < P4EST_CHILDREN; k++) { 3679 PetscInt l; 3680 PetscReal prod = 1.; 3681 3682 for (l = 0; l < P4EST_DIM; l++) prod *= (k & (1 << l)) ? eta[l] : (1. - eta[l]); 3683 for (l = 0; l < 3; l++) coordGuess[l] += prod * treeCoords[k][l]; 3684 } 3685 } 3686 } 3687 for (j = 0; j < 3; j++) coordP4est[j] = (double)eta[j]; 3688 3689 if (geom) { 3690 (geom->X)(geom, t, coordP4est, coordP4estMapped); 3691 for (j = 0; j < coordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j]; 3692 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded"); 3693 } 3694 } 3695 } 3696 PetscCall(DMPlexRestoreTransitiveClosure(plex, cell, PETSC_TRUE, &closureSize, &closure)); 3697 PetscFunctionReturn(PETSC_SUCCESS); 3698 } 3699 3700 static PetscErrorCode DMPforestMapCoordinates(DM dm, DM plex) 3701 { 3702 DM_Forest *forest; 3703 DM_Forest_pforest *pforest; 3704 p4est_geometry_t *geom; 3705 PetscInt cLocalStart, cLocalEnd; 3706 Vec coordLocalVec; 3707 PetscScalar *coords; 3708 p4est_topidx_t flt, llt, t; 3709 p4est_tree_t *trees; 3710 PetscErrorCode (*map)(DM, PetscInt, PetscInt, const PetscReal[], PetscReal[], void *); 3711 void *mapCtx; 3712 3713 PetscFunctionBegin; 3714 forest = (DM_Forest *)dm->data; 3715 pforest = (DM_Forest_pforest *)forest->data; 3716 geom = pforest->topo->geom; 3717 PetscCall(DMForestGetBaseCoordinateMapping(dm, &map, &mapCtx)); 3718 if (!geom && !map) PetscFunctionReturn(PETSC_SUCCESS); 3719 PetscCall(DMGetCoordinatesLocal(plex, &coordLocalVec)); 3720 PetscCall(VecGetArray(coordLocalVec, &coords)); 3721 cLocalStart = pforest->cLocalStart; 3722 cLocalEnd = pforest->cLocalEnd; 3723 flt = pforest->forest->first_local_tree; 3724 llt = pforest->forest->last_local_tree; 3725 trees = (p4est_tree_t *)pforest->forest->trees->array; 3726 if (map) { /* apply the map directly to the existing coordinates */ 3727 PetscSection coordSec; 3728 PetscInt coordStart, coordEnd, p, coordDim, p4estCoordDim, cStart, cEnd, cEndInterior; 3729 DM base; 3730 3731 PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd)); 3732 PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3733 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 3734 PetscCall(DMForestGetBaseDM(dm, &base)); 3735 PetscCall(DMGetCoordinateSection(plex, &coordSec)); 3736 PetscCall(PetscSectionGetChart(coordSec, &coordStart, &coordEnd)); 3737 PetscCall(DMGetCoordinateDim(plex, &coordDim)); 3738 p4estCoordDim = PetscMin(coordDim, 3); 3739 for (p = coordStart; p < coordEnd; p++) { 3740 PetscInt *star = NULL, starSize; 3741 PetscInt dof, off, cell = -1, coarsePoint = -1; 3742 PetscInt nCoords, i; 3743 PetscCall(PetscSectionGetDof(coordSec, p, &dof)); 3744 PetscCheck(dof % coordDim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not understand coordinate layout"); 3745 nCoords = dof / coordDim; 3746 PetscCall(PetscSectionGetOffset(coordSec, p, &off)); 3747 PetscCall(DMPlexGetTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star)); 3748 for (i = 0; i < starSize; i++) { 3749 PetscInt point = star[2 * i]; 3750 3751 if (cStart <= point && point < cEnd) { 3752 cell = point; 3753 break; 3754 } 3755 } 3756 PetscCall(DMPlexRestoreTransitiveClosure(plex, p, PETSC_FALSE, &starSize, &star)); 3757 if (cell >= 0) { 3758 if (cell < cLocalStart) { 3759 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3760 3761 coarsePoint = ghosts[cell].p.which_tree; 3762 } else if (cell < cLocalEnd) { 3763 cell -= cLocalStart; 3764 for (t = flt; t <= llt; t++) { 3765 p4est_tree_t *tree = &trees[t]; 3766 3767 if (cell >= tree->quadrants_offset && (size_t)cell < tree->quadrants_offset + tree->quadrants.elem_count) { 3768 coarsePoint = t; 3769 break; 3770 } 3771 } 3772 } else { 3773 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3774 3775 coarsePoint = ghosts[cell - cLocalEnd].p.which_tree; 3776 } 3777 } 3778 for (i = 0; i < nCoords; i++) { 3779 PetscScalar *coord = &coords[off + i * coordDim]; 3780 PetscReal coordP4est[3] = {0.}; 3781 PetscReal coordP4estMapped[3] = {0.}; 3782 PetscInt j; 3783 3784 for (j = 0; j < p4estCoordDim; j++) coordP4est[j] = PetscRealPart(coord[j]); 3785 PetscCall((map)(base, coarsePoint, p4estCoordDim, coordP4est, coordP4estMapped, mapCtx)); 3786 for (j = 0; j < p4estCoordDim; j++) coord[j] = (PetscScalar)coordP4estMapped[j]; 3787 } 3788 } 3789 } else { /* we have to transform coordinates back to the unit cube (where geom is defined), and then apply geom */ 3790 PetscInt cStart, cEnd, cEndInterior; 3791 3792 PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd)); 3793 PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3794 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 3795 if (cLocalStart > 0) { 3796 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3797 PetscInt count; 3798 3799 for (count = 0; count < cLocalStart; count++) { 3800 p4est_quadrant_t *quad = &ghosts[count]; 3801 p4est_topidx_t t = quad->p.which_tree; 3802 3803 PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, quad, t, pforest->topo->conn, coords)); 3804 } 3805 } 3806 for (t = flt; t <= llt; t++) { 3807 p4est_tree_t *tree = &trees[t]; 3808 PetscInt offset = cLocalStart + tree->quadrants_offset, i; 3809 PetscInt numQuads = (PetscInt)tree->quadrants.elem_count; 3810 p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array; 3811 3812 for (i = 0; i < numQuads; i++) { 3813 PetscInt count = i + offset; 3814 3815 PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count, &quads[i], t, pforest->topo->conn, coords)); 3816 } 3817 } 3818 if (cLocalEnd - cLocalStart < cEnd - cStart) { 3819 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3820 PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count; 3821 PetscInt count; 3822 3823 for (count = 0; count < numGhosts - cLocalStart; count++) { 3824 p4est_quadrant_t *quad = &ghosts[count + cLocalStart]; 3825 p4est_topidx_t t = quad->p.which_tree; 3826 3827 PetscCall(DMPforestMapCoordinates_Cell(plex, geom, count + cLocalEnd, quad, t, pforest->topo->conn, coords)); 3828 } 3829 } 3830 } 3831 PetscCall(VecRestoreArray(coordLocalVec, &coords)); 3832 PetscFunctionReturn(PETSC_SUCCESS); 3833 } 3834 3835 static PetscErrorCode PforestQuadrantIsInterior(p4est_quadrant_t *quad, PetscBool *is_interior) 3836 { 3837 PetscFunctionBegin; 3838 p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level); 3839 if ((quad->x > 0) && (quad->x + h < P4EST_ROOT_LEN) 3840 #ifdef P4_TO_P8 3841 && (quad->z > 0) && (quad->z + h < P4EST_ROOT_LEN) 3842 #endif 3843 && (quad->y > 0) && (quad->y + h < P4EST_ROOT_LEN)) { 3844 *is_interior = PETSC_TRUE; 3845 } else { 3846 *is_interior = PETSC_FALSE; 3847 } 3848 PetscFunctionReturn(PETSC_SUCCESS); 3849 } 3850 3851 /* We always use DG coordinates with p4est: if they do not match the vertex 3852 coordinates, add space for them in the section */ 3853 static PetscErrorCode PforestCheckLocalizeCell(DM plex, PetscInt cDim, Vec cVecOld, DM_Forest_pforest *pforest, PetscSection oldSection, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad) 3854 { 3855 PetscBool is_interior; 3856 3857 PetscFunctionBegin; 3858 PetscCall(PforestQuadrantIsInterior(quad, &is_interior)); 3859 if (is_interior) { // quads in the interior of a coarse cell can't touch periodic interfaces 3860 PetscCall(PetscSectionSetDof(newSection, cell, 0)); 3861 PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0)); 3862 } else { 3863 PetscInt cSize; 3864 PetscScalar *values = NULL; 3865 PetscBool same_coords = PETSC_TRUE; 3866 3867 PetscCall(DMPlexVecGetClosure(plex, oldSection, cVecOld, cell, &cSize, &values)); 3868 PetscAssert(cSize == cDim * P4EST_CHILDREN, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected closure size"); 3869 for (int c = 0; c < P4EST_CHILDREN; c++) { 3870 p4est_qcoord_t quad_coords[3]; 3871 p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level); 3872 double corner_coords[3]; 3873 double vert_coords[3]; 3874 PetscInt corner = PetscVertToP4estVert[c]; 3875 3876 for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) vert_coords[d] = PetscRealPart(values[c * cDim + d]); 3877 3878 quad_coords[0] = quad->x; 3879 quad_coords[1] = quad->y; 3880 #ifdef P4_TO_P8 3881 quad_coords[2] = quad->z; 3882 #endif 3883 for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0; 3884 #ifndef P4_TO_P8 3885 PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords)); 3886 #else 3887 PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords)); 3888 #endif 3889 for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) { 3890 if (fabs(vert_coords[d] - corner_coords[d]) > PETSC_SMALL) { 3891 same_coords = PETSC_FALSE; 3892 break; 3893 } 3894 } 3895 } 3896 if (same_coords) { 3897 PetscCall(PetscSectionSetDof(newSection, cell, 0)); 3898 PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, 0)); 3899 } else { 3900 PetscCall(PetscSectionSetDof(newSection, cell, cSize)); 3901 PetscCall(PetscSectionSetFieldDof(newSection, cell, 0, cSize)); 3902 } 3903 PetscCall(DMPlexVecRestoreClosure(plex, oldSection, cVecOld, cell, &cSize, &values)); 3904 } 3905 PetscFunctionReturn(PETSC_SUCCESS); 3906 } 3907 3908 static PetscErrorCode PforestLocalizeCell(DM plex, PetscInt cDim, DM_Forest_pforest *pforest, PetscSection newSection, PetscInt cell, PetscInt coarsePoint, p4est_quadrant_t *quad, PetscScalar coords[]) 3909 { 3910 PetscInt cdof, off; 3911 3912 PetscFunctionBegin; 3913 PetscCall(PetscSectionGetDof(newSection, cell, &cdof)); 3914 if (!cdof) PetscFunctionReturn(PETSC_SUCCESS); 3915 3916 PetscCall(PetscSectionGetOffset(newSection, cell, &off)); 3917 for (PetscInt c = 0, pos = off; c < P4EST_CHILDREN; c++) { 3918 p4est_qcoord_t quad_coords[3]; 3919 p4est_qcoord_t h = P4EST_QUADRANT_LEN(quad->level); 3920 double corner_coords[3]; 3921 PetscInt corner = PetscVertToP4estVert[c]; 3922 3923 quad_coords[0] = quad->x; 3924 quad_coords[1] = quad->y; 3925 #ifdef P4_TO_P8 3926 quad_coords[2] = quad->z; 3927 #endif 3928 for (int d = 0; d < 3; d++) quad_coords[d] += (corner & (1 << d)) ? h : 0; 3929 #ifndef P4_TO_P8 3930 PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], corner_coords)); 3931 #else 3932 PetscCallP4est(p4est_qcoord_to_vertex, (pforest->forest->connectivity, coarsePoint, quad_coords[0], quad_coords[1], quad_coords[2], corner_coords)); 3933 #endif 3934 for (PetscInt d = 0; d < PetscMin(cDim, 3); d++) coords[pos++] = corner_coords[d]; 3935 for (PetscInt d = PetscMin(cDim, 3); d < cDim; d++) coords[pos++] = 0.; 3936 } 3937 PetscFunctionReturn(PETSC_SUCCESS); 3938 } 3939 3940 static PetscErrorCode DMPforestLocalizeCoordinates(DM dm, DM plex) 3941 { 3942 DM_Forest *forest; 3943 DM_Forest_pforest *pforest; 3944 DM base, cdm, cdmCell; 3945 Vec cVec, cVecOld; 3946 PetscSection oldSection, newSection; 3947 PetscScalar *coords2; 3948 const PetscReal *L; 3949 PetscInt cLocalStart, cLocalEnd, coarsePoint; 3950 PetscInt cDim, newStart, newEnd; 3951 PetscInt v, vStart, vEnd, cp, cStart, cEnd, cEndInterior; 3952 p4est_topidx_t flt, llt, t; 3953 p4est_tree_t *trees; 3954 PetscBool baseLocalized = PETSC_FALSE; 3955 3956 PetscFunctionBegin; 3957 PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L)); 3958 /* we localize on all cells if we don't have a base DM or the base DM coordinates have not been localized */ 3959 PetscCall(DMGetCoordinateDim(dm, &cDim)); 3960 PetscCall(DMForestGetBaseDM(dm, &base)); 3961 if (base) PetscCall(DMGetCoordinatesLocalized(base, &baseLocalized)); 3962 if (!baseLocalized) base = NULL; 3963 if (!baseLocalized && !L) PetscFunctionReturn(PETSC_SUCCESS); 3964 PetscCall(DMPlexGetChart(plex, &newStart, &newEnd)); 3965 3966 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newSection)); 3967 PetscCall(PetscSectionSetNumFields(newSection, 1)); 3968 PetscCall(PetscSectionSetFieldComponents(newSection, 0, cDim)); 3969 PetscCall(PetscSectionSetChart(newSection, newStart, newEnd)); 3970 3971 PetscCall(DMGetCoordinateSection(plex, &oldSection)); 3972 PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd)); 3973 PetscCall(DMGetCoordinatesLocal(plex, &cVecOld)); 3974 3975 forest = (DM_Forest *)dm->data; 3976 pforest = (DM_Forest_pforest *)forest->data; 3977 cLocalStart = pforest->cLocalStart; 3978 cLocalEnd = pforest->cLocalEnd; 3979 flt = pforest->forest->first_local_tree; 3980 llt = pforest->forest->last_local_tree; 3981 trees = (p4est_tree_t *)pforest->forest->trees->array; 3982 3983 PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd)); 3984 PetscCall(DMPlexGetCellTypeStratum(plex, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL)); 3985 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 3986 cp = 0; 3987 if (cLocalStart > 0) { 3988 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 3989 PetscInt cell; 3990 3991 for (cell = 0; cell < cLocalStart; ++cell, cp++) { 3992 p4est_quadrant_t *quad = &ghosts[cell]; 3993 3994 coarsePoint = quad->p.which_tree; 3995 PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad)); 3996 } 3997 } 3998 for (t = flt; t <= llt; t++) { 3999 p4est_tree_t *tree = &trees[t]; 4000 PetscInt offset = cLocalStart + tree->quadrants_offset; 4001 PetscInt numQuads = (PetscInt)tree->quadrants.elem_count; 4002 p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array; 4003 PetscInt i; 4004 4005 if (!numQuads) continue; 4006 coarsePoint = t; 4007 for (i = 0; i < numQuads; i++, cp++) { 4008 PetscInt cell = i + offset; 4009 p4est_quadrant_t *quad = &quads[i]; 4010 4011 PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad)); 4012 } 4013 } 4014 if (cLocalEnd - cLocalStart < cEnd - cStart) { 4015 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 4016 PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count; 4017 PetscInt count; 4018 4019 for (count = 0; count < numGhosts - cLocalStart; count++, cp++) { 4020 p4est_quadrant_t *quad = &ghosts[count + cLocalStart]; 4021 coarsePoint = quad->p.which_tree; 4022 PetscInt cell = count + cLocalEnd; 4023 4024 PetscCall(PforestCheckLocalizeCell(plex, cDim, cVecOld, pforest, oldSection, newSection, cell, coarsePoint, quad)); 4025 } 4026 } 4027 PetscAssert(cp == cEnd - cStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of fine cells %" PetscInt_FMT " != %" PetscInt_FMT, cp, cEnd - cStart); 4028 4029 PetscCall(PetscSectionSetUp(newSection)); 4030 PetscCall(DMGetCoordinateDM(plex, &cdm)); 4031 PetscCall(DMClone(cdm, &cdmCell)); 4032 PetscCall(DMSetCellCoordinateDM(plex, cdmCell)); 4033 PetscCall(DMDestroy(&cdmCell)); 4034 PetscCall(DMSetCellCoordinateSection(plex, cDim, newSection)); 4035 PetscCall(PetscSectionGetStorageSize(newSection, &v)); 4036 PetscCall(VecCreate(PETSC_COMM_SELF, &cVec)); 4037 PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates")); 4038 PetscCall(VecSetBlockSize(cVec, cDim)); 4039 PetscCall(VecSetSizes(cVec, v, PETSC_DETERMINE)); 4040 PetscCall(VecSetType(cVec, VECSTANDARD)); 4041 PetscCall(VecSet(cVec, PETSC_MIN_REAL)); 4042 4043 /* Localize coordinates on cells if needed */ 4044 PetscCall(VecGetArray(cVec, &coords2)); 4045 cp = 0; 4046 if (cLocalStart > 0) { 4047 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 4048 PetscInt cell; 4049 4050 for (cell = 0; cell < cLocalStart; ++cell, cp++) { 4051 p4est_quadrant_t *quad = &ghosts[cell]; 4052 4053 coarsePoint = quad->p.which_tree; 4054 PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2)); 4055 } 4056 } 4057 for (t = flt; t <= llt; t++) { 4058 p4est_tree_t *tree = &trees[t]; 4059 PetscInt offset = cLocalStart + tree->quadrants_offset; 4060 PetscInt numQuads = (PetscInt)tree->quadrants.elem_count; 4061 p4est_quadrant_t *quads = (p4est_quadrant_t *)tree->quadrants.array; 4062 PetscInt i; 4063 4064 if (!numQuads) continue; 4065 coarsePoint = t; 4066 for (i = 0; i < numQuads; i++, cp++) { 4067 PetscInt cell = i + offset; 4068 p4est_quadrant_t *quad = &quads[i]; 4069 4070 PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2)); 4071 } 4072 } 4073 if (cLocalEnd - cLocalStart < cEnd - cStart) { 4074 p4est_quadrant_t *ghosts = (p4est_quadrant_t *)pforest->ghost->ghosts.array; 4075 PetscInt numGhosts = (PetscInt)pforest->ghost->ghosts.elem_count; 4076 PetscInt count; 4077 4078 for (count = 0; count < numGhosts - cLocalStart; count++, cp++) { 4079 p4est_quadrant_t *quad = &ghosts[count + cLocalStart]; 4080 coarsePoint = quad->p.which_tree; 4081 PetscInt cell = count + cLocalEnd; 4082 4083 PetscCall(PforestLocalizeCell(plex, cDim, pforest, newSection, cell, coarsePoint, quad, coords2)); 4084 } 4085 } 4086 PetscCall(VecRestoreArray(cVec, &coords2)); 4087 PetscCall(DMSetCellCoordinatesLocal(plex, cVec)); 4088 PetscCall(VecDestroy(&cVec)); 4089 PetscCall(PetscSectionDestroy(&newSection)); 4090 PetscFunctionReturn(PETSC_SUCCESS); 4091 } 4092 4093 #define DMForestClearAdaptivityForest_pforest _append_pforest(DMForestClearAdaptivityForest) 4094 static PetscErrorCode DMForestClearAdaptivityForest_pforest(DM dm) 4095 { 4096 DM_Forest *forest; 4097 DM_Forest_pforest *pforest; 4098 4099 PetscFunctionBegin; 4100 forest = (DM_Forest *)dm->data; 4101 pforest = (DM_Forest_pforest *)forest->data; 4102 PetscCall(PetscSFDestroy(&pforest->pointAdaptToSelfSF)); 4103 PetscCall(PetscSFDestroy(&pforest->pointSelfToAdaptSF)); 4104 PetscCall(PetscFree(pforest->pointAdaptToSelfCids)); 4105 PetscCall(PetscFree(pforest->pointSelfToAdaptCids)); 4106 PetscFunctionReturn(PETSC_SUCCESS); 4107 } 4108 4109 static PetscErrorCode DMConvert_pforest_plex(DM dm, DMType newtype, DM *plex) 4110 { 4111 DM_Forest *forest; 4112 DM_Forest_pforest *pforest; 4113 DM refTree, newPlex, base; 4114 PetscInt adjDim, adjCodim, coordDim; 4115 MPI_Comm comm; 4116 PetscBool isPforest; 4117 PetscInt dim; 4118 PetscInt overlap; 4119 p4est_connect_type_t ctype; 4120 p4est_locidx_t first_local_quad = -1; 4121 sc_array_t *points_per_dim, *cone_sizes, *cones, *cone_orientations, *coords, *children, *parents, *childids, *leaves, *remotes; 4122 PetscSection parentSection; 4123 PetscSF pointSF; 4124 size_t zz, count; 4125 PetscInt pStart, pEnd; 4126 DMLabel ghostLabelBase = NULL; 4127 4128 PetscFunctionBegin; 4129 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4130 comm = PetscObjectComm((PetscObject)dm); 4131 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPFOREST, &isPforest)); 4132 PetscCheck(isPforest, comm, PETSC_ERR_ARG_WRONG, "Expected DM type %s, got %s", DMPFOREST, ((PetscObject)dm)->type_name); 4133 PetscCall(DMSetUp(dm)); 4134 PetscCall(DMGetDimension(dm, &dim)); 4135 PetscCheck(dim == P4EST_DIM, comm, PETSC_ERR_ARG_WRONG, "Expected DM dimension %d, got %" PetscInt_FMT, P4EST_DIM, dim); 4136 forest = (DM_Forest *)dm->data; 4137 pforest = (DM_Forest_pforest *)forest->data; 4138 PetscCall(DMForestGetBaseDM(dm, &base)); 4139 if (base) PetscCall(DMGetLabel(base, "ghost", &ghostLabelBase)); 4140 if (!pforest->plex) { 4141 PetscMPIInt size; 4142 4143 PetscCallMPI(MPI_Comm_size(comm, &size)); 4144 PetscCall(DMCreate(comm, &newPlex)); 4145 PetscCall(DMSetType(newPlex, DMPLEX)); 4146 PetscCall(DMSetMatType(newPlex, dm->mattype)); 4147 /* share labels */ 4148 PetscCall(DMCopyLabels(dm, newPlex, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL)); 4149 PetscCall(DMForestGetAdjacencyDimension(dm, &adjDim)); 4150 PetscCall(DMForestGetAdjacencyCodimension(dm, &adjCodim)); 4151 PetscCall(DMGetCoordinateDim(dm, &coordDim)); 4152 if (adjDim == 0) { 4153 ctype = P4EST_CONNECT_FULL; 4154 } else if (adjCodim == 1) { 4155 ctype = P4EST_CONNECT_FACE; 4156 #if defined(P4_TO_P8) 4157 } else if (adjDim == 1) { 4158 ctype = P8EST_CONNECT_EDGE; 4159 #endif 4160 } else { 4161 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid adjacency dimension %" PetscInt_FMT, adjDim); 4162 } 4163 PetscCheck(ctype == P4EST_CONNECT_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Adjacency dimension %" PetscInt_FMT " / codimension %" PetscInt_FMT " not supported yet", adjDim, adjCodim); 4164 PetscCall(DMForestGetPartitionOverlap(dm, &overlap)); 4165 PetscCall(DMPlexSetOverlap_Plex(newPlex, NULL, overlap)); 4166 4167 points_per_dim = sc_array_new(sizeof(p4est_locidx_t)); 4168 cone_sizes = sc_array_new(sizeof(p4est_locidx_t)); 4169 cones = sc_array_new(sizeof(p4est_locidx_t)); 4170 cone_orientations = sc_array_new(sizeof(p4est_locidx_t)); 4171 coords = sc_array_new(3 * sizeof(double)); 4172 children = sc_array_new(sizeof(p4est_locidx_t)); 4173 parents = sc_array_new(sizeof(p4est_locidx_t)); 4174 childids = sc_array_new(sizeof(p4est_locidx_t)); 4175 leaves = sc_array_new(sizeof(p4est_locidx_t)); 4176 remotes = sc_array_new(2 * sizeof(p4est_locidx_t)); 4177 4178 PetscCallP4est(p4est_get_plex_data_ext, (pforest->forest, &pforest->ghost, &pforest->lnodes, ctype, (int)((size > 1) ? overlap : 0), &first_local_quad, points_per_dim, cone_sizes, cones, cone_orientations, coords, children, parents, childids, leaves, remotes, 1)); 4179 4180 pforest->cLocalStart = (PetscInt)first_local_quad; 4181 pforest->cLocalEnd = pforest->cLocalStart + (PetscInt)pforest->forest->local_num_quadrants; 4182 PetscCall(locidx_to_PetscInt(points_per_dim)); 4183 PetscCall(locidx_to_PetscInt(cone_sizes)); 4184 PetscCall(locidx_to_PetscInt(cones)); 4185 PetscCall(locidx_to_PetscInt(cone_orientations)); 4186 PetscCall(coords_double_to_PetscScalar(coords, coordDim)); 4187 PetscCall(locidx_to_PetscInt(children)); 4188 PetscCall(locidx_to_PetscInt(parents)); 4189 PetscCall(locidx_to_PetscInt(childids)); 4190 PetscCall(locidx_to_PetscInt(leaves)); 4191 PetscCall(locidx_pair_to_PetscSFNode(remotes)); 4192 4193 PetscCall(DMSetDimension(newPlex, P4EST_DIM)); 4194 PetscCall(DMSetCoordinateDim(newPlex, coordDim)); 4195 PetscCall(DMPlexSetMaxProjectionHeight(newPlex, P4EST_DIM - 1)); 4196 PetscCall(DMPlexCreateFromDAG(newPlex, P4EST_DIM, (PetscInt *)points_per_dim->array, (PetscInt *)cone_sizes->array, (PetscInt *)cones->array, (PetscInt *)cone_orientations->array, (PetscScalar *)coords->array)); 4197 PetscCall(DMPlexConvertOldOrientations_Internal(newPlex)); 4198 PetscCall(DMCreateReferenceTree_pforest(comm, &refTree)); 4199 PetscCall(DMPlexSetReferenceTree(newPlex, refTree)); 4200 PetscCall(PetscSectionCreate(comm, &parentSection)); 4201 PetscCall(DMPlexGetChart(newPlex, &pStart, &pEnd)); 4202 PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd)); 4203 count = children->elem_count; 4204 for (zz = 0; zz < count; zz++) { 4205 PetscInt child = *((PetscInt *)sc_array_index(children, zz)); 4206 4207 PetscCall(PetscSectionSetDof(parentSection, child, 1)); 4208 } 4209 PetscCall(PetscSectionSetUp(parentSection)); 4210 PetscCall(DMPlexSetTree(newPlex, parentSection, (PetscInt *)parents->array, (PetscInt *)childids->array)); 4211 PetscCall(PetscSectionDestroy(&parentSection)); 4212 PetscCall(PetscSFCreate(comm, &pointSF)); 4213 /* 4214 These arrays defining the sf are from the p4est library, but the code there shows the leaves being populated in increasing order. 4215 https://gitlab.com/petsc/petsc/merge_requests/2248#note_240186391 4216 */ 4217 PetscCall(PetscSFSetGraph(pointSF, pEnd - pStart, (PetscInt)leaves->elem_count, (PetscInt *)leaves->array, PETSC_COPY_VALUES, (PetscSFNode *)remotes->array, PETSC_COPY_VALUES)); 4218 PetscCall(DMSetPointSF(newPlex, pointSF)); 4219 PetscCall(DMSetPointSF(dm, pointSF)); 4220 { 4221 DM coordDM; 4222 4223 PetscCall(DMGetCoordinateDM(newPlex, &coordDM)); 4224 PetscCall(DMSetPointSF(coordDM, pointSF)); 4225 } 4226 PetscCall(PetscSFDestroy(&pointSF)); 4227 sc_array_destroy(points_per_dim); 4228 sc_array_destroy(cone_sizes); 4229 sc_array_destroy(cones); 4230 sc_array_destroy(cone_orientations); 4231 sc_array_destroy(coords); 4232 sc_array_destroy(children); 4233 sc_array_destroy(parents); 4234 sc_array_destroy(childids); 4235 sc_array_destroy(leaves); 4236 sc_array_destroy(remotes); 4237 4238 { 4239 const PetscReal *maxCell, *Lstart, *L; 4240 4241 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 4242 PetscCall(DMSetPeriodicity(newPlex, maxCell, Lstart, L)); 4243 PetscCall(DMPforestLocalizeCoordinates(dm, newPlex)); 4244 } 4245 4246 if (overlap > 0) { /* the p4est routine can't set all of the coordinates in its routine if there is overlap */ 4247 Vec coordsGlobal, coordsLocal; 4248 const PetscScalar *globalArray; 4249 PetscScalar *localArray; 4250 PetscSF coordSF; 4251 DM coordDM; 4252 4253 PetscCall(DMGetCoordinateDM(newPlex, &coordDM)); 4254 PetscCall(DMGetSectionSF(coordDM, &coordSF)); 4255 PetscCall(DMGetCoordinates(newPlex, &coordsGlobal)); 4256 PetscCall(DMGetCoordinatesLocal(newPlex, &coordsLocal)); 4257 PetscCall(VecGetArrayRead(coordsGlobal, &globalArray)); 4258 PetscCall(VecGetArray(coordsLocal, &localArray)); 4259 PetscCall(PetscSFBcastBegin(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE)); 4260 PetscCall(PetscSFBcastEnd(coordSF, MPIU_SCALAR, globalArray, localArray, MPI_REPLACE)); 4261 PetscCall(VecRestoreArray(coordsLocal, &localArray)); 4262 PetscCall(VecRestoreArrayRead(coordsGlobal, &globalArray)); 4263 PetscCall(DMSetCoordinatesLocal(newPlex, coordsLocal)); 4264 } 4265 PetscCall(DMPforestMapCoordinates(dm, newPlex)); 4266 4267 pforest->plex = newPlex; 4268 4269 /* copy labels */ 4270 PetscCall(DMPforestLabelsFinalize(dm, newPlex)); 4271 4272 if (ghostLabelBase || pforest->ghostName) { /* we have to do this after copying labels because the labels drive the construction of ghost cells */ 4273 PetscInt numAdded; 4274 DM newPlexGhosted; 4275 void *ctx; 4276 4277 PetscCall(DMPlexConstructGhostCells(newPlex, pforest->ghostName, &numAdded, &newPlexGhosted)); 4278 PetscCall(DMGetApplicationContext(newPlex, &ctx)); 4279 PetscCall(DMSetApplicationContext(newPlexGhosted, ctx)); 4280 /* we want the sf for the ghost dm to be the one for the p4est dm as well */ 4281 PetscCall(DMGetPointSF(newPlexGhosted, &pointSF)); 4282 PetscCall(DMSetPointSF(dm, pointSF)); 4283 PetscCall(DMDestroy(&newPlex)); 4284 PetscCall(DMPlexSetReferenceTree(newPlexGhosted, refTree)); 4285 PetscCall(DMForestClearAdaptivityForest_pforest(dm)); 4286 newPlex = newPlexGhosted; 4287 4288 /* share the labels back */ 4289 PetscCall(DMDestroyLabelLinkList_Internal(dm)); 4290 PetscCall(DMCopyLabels(newPlex, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL)); 4291 pforest->plex = newPlex; 4292 } 4293 PetscCall(DMDestroy(&refTree)); 4294 if (dm->setfromoptionscalled) { 4295 PetscObjectOptionsBegin((PetscObject)newPlex); 4296 PetscCall(DMSetFromOptions_NonRefinement_Plex(newPlex, PetscOptionsObject)); 4297 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)newPlex, PetscOptionsObject)); 4298 PetscOptionsEnd(); 4299 } 4300 PetscCall(DMViewFromOptions(newPlex, NULL, "-dm_p4est_plex_view")); 4301 { 4302 DM cdm; 4303 PetscSection coordsSec; 4304 Vec coords; 4305 PetscInt cDim; 4306 4307 PetscCall(DMGetCoordinateDim(newPlex, &cDim)); 4308 PetscCall(DMGetCoordinateSection(newPlex, &coordsSec)); 4309 PetscCall(DMSetCoordinateSection(dm, cDim, coordsSec)); 4310 PetscCall(DMGetCoordinatesLocal(newPlex, &coords)); 4311 PetscCall(DMSetCoordinatesLocal(dm, coords)); 4312 PetscCall(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 PetscCall(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->createglobalvector = DMCreateGlobalVector_pforest; 5134 dm->ops->createlocalvector = DMCreateLocalVector_pforest; 5135 dm->ops->creatematrix = DMCreateMatrix_pforest; 5136 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_pforest; 5137 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_pforest; 5138 dm->ops->projectfieldlocal = DMProjectFieldLocal_pforest; 5139 dm->ops->createlocalsection = DMCreatelocalsection_pforest; 5140 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_pforest; 5141 dm->ops->computel2diff = DMComputeL2Diff_pforest; 5142 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_pforest; 5143 dm->ops->getdimpoints = DMGetDimPoints_pforest; 5144 dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_pforest; 5145 dm->ops->createddscatters = DMCreateDomainDecompositionScatters_pforest; 5146 5147 PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_plex_pforest) "_C", DMConvert_plex_pforest)); 5148 PetscCall(PetscObjectComposeFunction((PetscObject)dm, PetscStringize(DMConvert_pforest_plex) "_C", DMConvert_pforest_plex)); 5149 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", DMCreateNeumannOverlap_pforest)); 5150 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", DMForestGetPartitionOverlap)); 5151 PetscFunctionReturn(PETSC_SUCCESS); 5152 } 5153 5154 #define DMCreate_pforest _append_pforest(DMCreate) 5155 PETSC_EXTERN PetscErrorCode DMCreate_pforest(DM dm) 5156 { 5157 DM_Forest *forest; 5158 DM_Forest_pforest *pforest; 5159 5160 PetscFunctionBegin; 5161 PetscCall(PetscP4estInitialize()); 5162 PetscCall(DMCreate_Forest(dm)); 5163 PetscCall(DMInitialize_pforest(dm)); 5164 PetscCall(DMSetDimension(dm, P4EST_DIM)); 5165 5166 /* set forest defaults */ 5167 PetscCall(DMForestSetTopology(dm, "unit")); 5168 PetscCall(DMForestSetMinimumRefinement(dm, 0)); 5169 PetscCall(DMForestSetInitialRefinement(dm, 0)); 5170 PetscCall(DMForestSetMaximumRefinement(dm, P4EST_QMAXLEVEL)); 5171 PetscCall(DMForestSetGradeFactor(dm, 2)); 5172 PetscCall(DMForestSetAdjacencyDimension(dm, 0)); 5173 PetscCall(DMForestSetPartitionOverlap(dm, 0)); 5174 5175 /* create p4est data */ 5176 PetscCall(PetscNew(&pforest)); 5177 5178 forest = (DM_Forest *)dm->data; 5179 forest->data = pforest; 5180 forest->destroy = DMForestDestroy_pforest; 5181 forest->ftemplate = DMForestTemplate_pforest; 5182 forest->transfervec = DMForestTransferVec_pforest; 5183 forest->transfervecfrombase = DMForestTransferVecFromBase_pforest; 5184 forest->createcellchart = DMForestCreateCellChart_pforest; 5185 forest->createcellsf = DMForestCreateCellSF_pforest; 5186 forest->clearadaptivityforest = DMForestClearAdaptivityForest_pforest; 5187 forest->getadaptivitysuccess = DMForestGetAdaptivitySuccess_pforest; 5188 pforest->topo = NULL; 5189 pforest->forest = NULL; 5190 pforest->ghost = NULL; 5191 pforest->lnodes = NULL; 5192 pforest->partition_for_coarsening = PETSC_TRUE; 5193 pforest->coarsen_hierarchy = PETSC_FALSE; 5194 pforest->cLocalStart = -1; 5195 pforest->cLocalEnd = -1; 5196 pforest->labelsFinalized = PETSC_FALSE; 5197 pforest->ghostName = NULL; 5198 PetscFunctionReturn(PETSC_SUCCESS); 5199 } 5200 5201 #endif /* defined(PETSC_HAVE_P4EST) */ 5202