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