1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/ 3 4 /*@C 5 DMPlexSetAdjacencyUser - Define adjacency in the mesh using a user-provided callback 6 7 Input Parameters: 8 + dm - The DM object 9 . user - The user callback, may be NULL (to clear the callback) 10 - ctx - context for callback evaluation, may be NULL 11 12 Level: advanced 13 14 Notes: 15 The caller of DMPlexGetAdjacency may need to arrange that a large enough array is available for the adjacency. 16 17 Any setting here overrides other configuration of DMPlex adjacency determination. 18 19 .seealso: DMSetAdjacency(), DMPlexDistribute(), DMPlexPreallocateOperator(), DMPlexGetAdjacency(), DMPlexGetAdjacencyUser() 20 @*/ 21 PetscErrorCode DMPlexSetAdjacencyUser(DM dm,PetscErrorCode (*user)(DM,PetscInt,PetscInt*,PetscInt[],void*),void *ctx) 22 { 23 DM_Plex *mesh = (DM_Plex *)dm->data; 24 25 PetscFunctionBegin; 26 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 27 mesh->useradjacency = user; 28 mesh->useradjacencyctx = ctx; 29 PetscFunctionReturn(0); 30 } 31 32 /*@C 33 DMPlexGetAdjacencyUser - get the user-defined adjacency callback 34 35 Input Parameter: 36 . dm - The DM object 37 38 Output Parameters: 39 + user - The user callback 40 - ctx - context for callback evaluation 41 42 Level: advanced 43 44 .seealso: DMSetAdjacency(), DMPlexDistribute(), DMPlexPreallocateOperator(), DMPlexGetAdjacency(), DMPlexSetAdjacencyUser() 45 @*/ 46 PetscErrorCode DMPlexGetAdjacencyUser(DM dm, PetscErrorCode (**user)(DM,PetscInt,PetscInt*,PetscInt[],void*), void **ctx) 47 { 48 DM_Plex *mesh = (DM_Plex *)dm->data; 49 50 PetscFunctionBegin; 51 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 52 if (user) *user = mesh->useradjacency; 53 if (ctx) *ctx = mesh->useradjacencyctx; 54 PetscFunctionReturn(0); 55 } 56 57 /*@ 58 DMPlexSetAdjacencyUseAnchors - Define adjacency in the mesh using the point-to-point constraints. 59 60 Input Parameters: 61 + dm - The DM object 62 - useAnchors - Flag to use the constraints. If PETSC_TRUE, then constrained points are omitted from DMPlexGetAdjacency(), and their anchor points appear in their place. 63 64 Level: intermediate 65 66 .seealso: DMGetAdjacency(), DMSetAdjacency(), DMPlexDistribute(), DMPlexPreallocateOperator(), DMPlexSetAnchors() 67 @*/ 68 PetscErrorCode DMPlexSetAdjacencyUseAnchors(DM dm, PetscBool useAnchors) 69 { 70 DM_Plex *mesh = (DM_Plex *) dm->data; 71 72 PetscFunctionBegin; 73 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 74 mesh->useAnchors = useAnchors; 75 PetscFunctionReturn(0); 76 } 77 78 /*@ 79 DMPlexGetAdjacencyUseAnchors - Query whether adjacency in the mesh uses the point-to-point constraints. 80 81 Input Parameter: 82 . dm - The DM object 83 84 Output Parameter: 85 . useAnchors - Flag to use the closure. If PETSC_TRUE, then constrained points are omitted from DMPlexGetAdjacency(), and their anchor points appear in their place. 86 87 Level: intermediate 88 89 .seealso: DMPlexSetAdjacencyUseAnchors(), DMSetAdjacency(), DMGetAdjacency(), DMPlexDistribute(), DMPlexPreallocateOperator(), DMPlexSetAnchors() 90 @*/ 91 PetscErrorCode DMPlexGetAdjacencyUseAnchors(DM dm, PetscBool *useAnchors) 92 { 93 DM_Plex *mesh = (DM_Plex *) dm->data; 94 95 PetscFunctionBegin; 96 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 97 PetscValidBoolPointer(useAnchors, 2); 98 *useAnchors = mesh->useAnchors; 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode DMPlexGetAdjacency_Cone_Internal(DM dm, PetscInt p, PetscInt *adjSize, PetscInt adj[]) 103 { 104 const PetscInt *cone = NULL; 105 PetscInt numAdj = 0, maxAdjSize = *adjSize, coneSize, c; 106 107 PetscFunctionBeginHot; 108 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 109 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 110 for (c = 0; c <= coneSize; ++c) { 111 const PetscInt point = !c ? p : cone[c-1]; 112 const PetscInt *support = NULL; 113 PetscInt supportSize, s, q; 114 115 CHKERRQ(DMPlexGetSupportSize(dm, point, &supportSize)); 116 CHKERRQ(DMPlexGetSupport(dm, point, &support)); 117 for (s = 0; s < supportSize; ++s) { 118 for (q = 0; q < numAdj || ((void)(adj[numAdj++] = support[s]),0); ++q) { 119 if (support[s] == adj[q]) break; 120 } 121 PetscCheckFalse(numAdj > maxAdjSize,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize); 122 } 123 } 124 *adjSize = numAdj; 125 PetscFunctionReturn(0); 126 } 127 128 static PetscErrorCode DMPlexGetAdjacency_Support_Internal(DM dm, PetscInt p, PetscInt *adjSize, PetscInt adj[]) 129 { 130 const PetscInt *support = NULL; 131 PetscInt numAdj = 0, maxAdjSize = *adjSize, supportSize, s; 132 133 PetscFunctionBeginHot; 134 CHKERRQ(DMPlexGetSupportSize(dm, p, &supportSize)); 135 CHKERRQ(DMPlexGetSupport(dm, p, &support)); 136 for (s = 0; s <= supportSize; ++s) { 137 const PetscInt point = !s ? p : support[s-1]; 138 const PetscInt *cone = NULL; 139 PetscInt coneSize, c, q; 140 141 CHKERRQ(DMPlexGetConeSize(dm, point, &coneSize)); 142 CHKERRQ(DMPlexGetCone(dm, point, &cone)); 143 for (c = 0; c < coneSize; ++c) { 144 for (q = 0; q < numAdj || ((void)(adj[numAdj++] = cone[c]),0); ++q) { 145 if (cone[c] == adj[q]) break; 146 } 147 PetscCheckFalse(numAdj > maxAdjSize,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize); 148 } 149 } 150 *adjSize = numAdj; 151 PetscFunctionReturn(0); 152 } 153 154 static PetscErrorCode DMPlexGetAdjacency_Transitive_Internal(DM dm, PetscInt p, PetscBool useClosure, PetscInt *adjSize, PetscInt adj[]) 155 { 156 PetscInt *star = NULL; 157 PetscInt numAdj = 0, maxAdjSize = *adjSize, starSize, s; 158 159 PetscFunctionBeginHot; 160 CHKERRQ(DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, &star)); 161 for (s = 0; s < starSize*2; s += 2) { 162 const PetscInt *closure = NULL; 163 PetscInt closureSize, c, q; 164 165 CHKERRQ(DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure)); 166 for (c = 0; c < closureSize*2; c += 2) { 167 for (q = 0; q < numAdj || ((void)(adj[numAdj++] = closure[c]),0); ++q) { 168 if (closure[c] == adj[q]) break; 169 } 170 PetscCheckFalse(numAdj > maxAdjSize,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize); 171 } 172 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure)); 173 } 174 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, p, useClosure, &starSize, &star)); 175 *adjSize = numAdj; 176 PetscFunctionReturn(0); 177 } 178 179 PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useCone, PetscBool useTransitiveClosure, PetscBool useAnchors, PetscInt *adjSize, PetscInt *adj[]) 180 { 181 static PetscInt asiz = 0; 182 PetscInt maxAnchors = 1; 183 PetscInt aStart = -1, aEnd = -1; 184 PetscInt maxAdjSize; 185 PetscSection aSec = NULL; 186 IS aIS = NULL; 187 const PetscInt *anchors; 188 DM_Plex *mesh = (DM_Plex *)dm->data; 189 190 PetscFunctionBeginHot; 191 if (useAnchors) { 192 CHKERRQ(DMPlexGetAnchors(dm,&aSec,&aIS)); 193 if (aSec) { 194 CHKERRQ(PetscSectionGetMaxDof(aSec,&maxAnchors)); 195 maxAnchors = PetscMax(1,maxAnchors); 196 CHKERRQ(PetscSectionGetChart(aSec,&aStart,&aEnd)); 197 CHKERRQ(ISGetIndices(aIS,&anchors)); 198 } 199 } 200 if (!*adj) { 201 PetscInt depth, coneSeries, supportSeries, maxC, maxS, pStart, pEnd; 202 203 CHKERRQ(DMPlexGetChart(dm, &pStart,&pEnd)); 204 CHKERRQ(DMPlexGetDepth(dm, &depth)); 205 depth = PetscMax(depth, -depth); 206 CHKERRQ(DMPlexGetMaxSizes(dm, &maxC, &maxS)); 207 coneSeries = (maxC > 1) ? ((PetscPowInt(maxC,depth+1)-1)/(maxC-1)) : depth+1; 208 supportSeries = (maxS > 1) ? ((PetscPowInt(maxS,depth+1)-1)/(maxS-1)) : depth+1; 209 asiz = PetscMax(PetscPowInt(maxS,depth)*coneSeries,PetscPowInt(maxC,depth)*supportSeries); 210 asiz *= maxAnchors; 211 asiz = PetscMin(asiz,pEnd-pStart); 212 CHKERRQ(PetscMalloc1(asiz,adj)); 213 } 214 if (*adjSize < 0) *adjSize = asiz; 215 maxAdjSize = *adjSize; 216 if (mesh->useradjacency) { 217 CHKERRQ((*mesh->useradjacency)(dm, p, adjSize, *adj, mesh->useradjacencyctx)); 218 } else if (useTransitiveClosure) { 219 CHKERRQ(DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, *adj)); 220 } else if (useCone) { 221 CHKERRQ(DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, *adj)); 222 } else { 223 CHKERRQ(DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, *adj)); 224 } 225 if (useAnchors && aSec) { 226 PetscInt origSize = *adjSize; 227 PetscInt numAdj = origSize; 228 PetscInt i = 0, j; 229 PetscInt *orig = *adj; 230 231 while (i < origSize) { 232 PetscInt p = orig[i]; 233 PetscInt aDof = 0; 234 235 if (p >= aStart && p < aEnd) { 236 CHKERRQ(PetscSectionGetDof(aSec,p,&aDof)); 237 } 238 if (aDof) { 239 PetscInt aOff; 240 PetscInt s, q; 241 242 for (j = i + 1; j < numAdj; j++) { 243 orig[j - 1] = orig[j]; 244 } 245 origSize--; 246 numAdj--; 247 CHKERRQ(PetscSectionGetOffset(aSec,p,&aOff)); 248 for (s = 0; s < aDof; ++s) { 249 for (q = 0; q < numAdj || ((void)(orig[numAdj++] = anchors[aOff+s]),0); ++q) { 250 if (anchors[aOff+s] == orig[q]) break; 251 } 252 PetscCheckFalse(numAdj > maxAdjSize,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize); 253 } 254 } 255 else { 256 i++; 257 } 258 } 259 *adjSize = numAdj; 260 CHKERRQ(ISRestoreIndices(aIS,&anchors)); 261 } 262 PetscFunctionReturn(0); 263 } 264 265 /*@ 266 DMPlexGetAdjacency - Return all points adjacent to the given point 267 268 Input Parameters: 269 + dm - The DM object 270 - p - The point 271 272 Input/Output Parameters: 273 + adjSize - The maximum size of adj if it is non-NULL, or PETSC_DETERMINE; 274 on output the number of adjacent points 275 - adj - Either NULL so that the array is allocated, or an existing array with size adjSize; 276 on output contains the adjacent points 277 278 Level: advanced 279 280 Notes: 281 The user must PetscFree the adj array if it was not passed in. 282 283 .seealso: DMSetAdjacency(), DMPlexDistribute(), DMCreateMatrix(), DMPlexPreallocateOperator() 284 @*/ 285 PetscErrorCode DMPlexGetAdjacency(DM dm, PetscInt p, PetscInt *adjSize, PetscInt *adj[]) 286 { 287 PetscBool useCone, useClosure, useAnchors; 288 289 PetscFunctionBeginHot; 290 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 291 PetscValidPointer(adjSize,3); 292 PetscValidPointer(adj,4); 293 CHKERRQ(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 294 CHKERRQ(DMPlexGetAdjacencyUseAnchors(dm, &useAnchors)); 295 CHKERRQ(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, adjSize, adj)); 296 PetscFunctionReturn(0); 297 } 298 299 /*@ 300 DMPlexCreateTwoSidedProcessSF - Create an SF which just has process connectivity 301 302 Collective on dm 303 304 Input Parameters: 305 + dm - The DM 306 . sfPoint - The PetscSF which encodes point connectivity 307 . rootRankSection - 308 . rootRanks - 309 . leftRankSection - 310 - leafRanks - 311 312 Output Parameters: 313 + processRanks - A list of process neighbors, or NULL 314 - sfProcess - An SF encoding the two-sided process connectivity, or NULL 315 316 Level: developer 317 318 .seealso: PetscSFCreate(), DMPlexCreateProcessSF() 319 @*/ 320 PetscErrorCode DMPlexCreateTwoSidedProcessSF(DM dm, PetscSF sfPoint, PetscSection rootRankSection, IS rootRanks, PetscSection leafRankSection, IS leafRanks, IS *processRanks, PetscSF *sfProcess) 321 { 322 const PetscSFNode *remotePoints; 323 PetscInt *localPointsNew; 324 PetscSFNode *remotePointsNew; 325 const PetscInt *nranks; 326 PetscInt *ranksNew; 327 PetscBT neighbors; 328 PetscInt pStart, pEnd, p, numLeaves, l, numNeighbors, n; 329 PetscMPIInt size, proc, rank; 330 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 333 PetscValidHeaderSpecific(sfPoint, PETSCSF_CLASSID, 2); 334 if (processRanks) {PetscValidPointer(processRanks, 7);} 335 if (sfProcess) {PetscValidPointer(sfProcess, 8);} 336 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 337 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 338 CHKERRQ(PetscSFGetGraph(sfPoint, NULL, &numLeaves, NULL, &remotePoints)); 339 CHKERRQ(PetscBTCreate(size, &neighbors)); 340 CHKERRQ(PetscBTMemzero(size, neighbors)); 341 /* Compute root-to-leaf process connectivity */ 342 CHKERRQ(PetscSectionGetChart(rootRankSection, &pStart, &pEnd)); 343 CHKERRQ(ISGetIndices(rootRanks, &nranks)); 344 for (p = pStart; p < pEnd; ++p) { 345 PetscInt ndof, noff, n; 346 347 CHKERRQ(PetscSectionGetDof(rootRankSection, p, &ndof)); 348 CHKERRQ(PetscSectionGetOffset(rootRankSection, p, &noff)); 349 for (n = 0; n < ndof; ++n) CHKERRQ(PetscBTSet(neighbors, nranks[noff+n])); 350 } 351 CHKERRQ(ISRestoreIndices(rootRanks, &nranks)); 352 /* Compute leaf-to-neighbor process connectivity */ 353 CHKERRQ(PetscSectionGetChart(leafRankSection, &pStart, &pEnd)); 354 CHKERRQ(ISGetIndices(leafRanks, &nranks)); 355 for (p = pStart; p < pEnd; ++p) { 356 PetscInt ndof, noff, n; 357 358 CHKERRQ(PetscSectionGetDof(leafRankSection, p, &ndof)); 359 CHKERRQ(PetscSectionGetOffset(leafRankSection, p, &noff)); 360 for (n = 0; n < ndof; ++n) CHKERRQ(PetscBTSet(neighbors, nranks[noff+n])); 361 } 362 CHKERRQ(ISRestoreIndices(leafRanks, &nranks)); 363 /* Compute leaf-to-root process connectivity */ 364 for (l = 0; l < numLeaves; ++l) {PetscBTSet(neighbors, remotePoints[l].rank);} 365 /* Calculate edges */ 366 PetscBTClear(neighbors, rank); 367 for (proc = 0, numNeighbors = 0; proc < size; ++proc) {if (PetscBTLookup(neighbors, proc)) ++numNeighbors;} 368 CHKERRQ(PetscMalloc1(numNeighbors, &ranksNew)); 369 CHKERRQ(PetscMalloc1(numNeighbors, &localPointsNew)); 370 CHKERRQ(PetscMalloc1(numNeighbors, &remotePointsNew)); 371 for (proc = 0, n = 0; proc < size; ++proc) { 372 if (PetscBTLookup(neighbors, proc)) { 373 ranksNew[n] = proc; 374 localPointsNew[n] = proc; 375 remotePointsNew[n].index = rank; 376 remotePointsNew[n].rank = proc; 377 ++n; 378 } 379 } 380 CHKERRQ(PetscBTDestroy(&neighbors)); 381 if (processRanks) CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numNeighbors, ranksNew, PETSC_OWN_POINTER, processRanks)); 382 else CHKERRQ(PetscFree(ranksNew)); 383 if (sfProcess) { 384 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess)); 385 CHKERRQ(PetscObjectSetName((PetscObject) *sfProcess, "Two-Sided Process SF")); 386 CHKERRQ(PetscSFSetFromOptions(*sfProcess)); 387 CHKERRQ(PetscSFSetGraph(*sfProcess, size, numNeighbors, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 388 } 389 PetscFunctionReturn(0); 390 } 391 392 /*@ 393 DMPlexDistributeOwnership - Compute owner information for shared points. This basically gets two-sided for an SF. 394 395 Collective on dm 396 397 Input Parameter: 398 . dm - The DM 399 400 Output Parameters: 401 + rootSection - The number of leaves for a given root point 402 . rootrank - The rank of each edge into the root point 403 . leafSection - The number of processes sharing a given leaf point 404 - leafrank - The rank of each process sharing a leaf point 405 406 Level: developer 407 408 .seealso: DMPlexCreateOverlapLabel(), DMPlexDistribute(), DMPlexDistributeOverlap() 409 @*/ 410 PetscErrorCode DMPlexDistributeOwnership(DM dm, PetscSection rootSection, IS *rootrank, PetscSection leafSection, IS *leafrank) 411 { 412 MPI_Comm comm; 413 PetscSF sfPoint; 414 const PetscInt *rootdegree; 415 PetscInt *myrank, *remoterank; 416 PetscInt pStart, pEnd, p, nedges; 417 PetscMPIInt rank; 418 419 PetscFunctionBegin; 420 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 421 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 422 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 423 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 424 /* Compute number of leaves for each root */ 425 CHKERRQ(PetscObjectSetName((PetscObject) rootSection, "Root Section")); 426 CHKERRQ(PetscSectionSetChart(rootSection, pStart, pEnd)); 427 CHKERRQ(PetscSFComputeDegreeBegin(sfPoint, &rootdegree)); 428 CHKERRQ(PetscSFComputeDegreeEnd(sfPoint, &rootdegree)); 429 for (p = pStart; p < pEnd; ++p) CHKERRQ(PetscSectionSetDof(rootSection, p, rootdegree[p-pStart])); 430 CHKERRQ(PetscSectionSetUp(rootSection)); 431 /* Gather rank of each leaf to root */ 432 CHKERRQ(PetscSectionGetStorageSize(rootSection, &nedges)); 433 CHKERRQ(PetscMalloc1(pEnd-pStart, &myrank)); 434 CHKERRQ(PetscMalloc1(nedges, &remoterank)); 435 for (p = 0; p < pEnd-pStart; ++p) myrank[p] = rank; 436 CHKERRQ(PetscSFGatherBegin(sfPoint, MPIU_INT, myrank, remoterank)); 437 CHKERRQ(PetscSFGatherEnd(sfPoint, MPIU_INT, myrank, remoterank)); 438 CHKERRQ(PetscFree(myrank)); 439 CHKERRQ(ISCreateGeneral(comm, nedges, remoterank, PETSC_OWN_POINTER, rootrank)); 440 /* Distribute remote ranks to leaves */ 441 CHKERRQ(PetscObjectSetName((PetscObject) leafSection, "Leaf Section")); 442 CHKERRQ(DMPlexDistributeFieldIS(dm, sfPoint, rootSection, *rootrank, leafSection, leafrank)); 443 PetscFunctionReturn(0); 444 } 445 446 /*@C 447 DMPlexCreateOverlapLabel - Compute owner information for shared points. This basically gets two-sided for an SF. 448 449 Collective on dm 450 451 Input Parameters: 452 + dm - The DM 453 . levels - Number of overlap levels 454 . rootSection - The number of leaves for a given root point 455 . rootrank - The rank of each edge into the root point 456 . leafSection - The number of processes sharing a given leaf point 457 - leafrank - The rank of each process sharing a leaf point 458 459 Output Parameter: 460 . ovLabel - DMLabel containing remote overlap contributions as point/rank pairings 461 462 Level: developer 463 464 .seealso: DMPlexDistributeOwnership(), DMPlexDistribute() 465 @*/ 466 PetscErrorCode DMPlexCreateOverlapLabel(DM dm, PetscInt levels, PetscSection rootSection, IS rootrank, PetscSection leafSection, IS leafrank, DMLabel *ovLabel) 467 { 468 MPI_Comm comm; 469 DMLabel ovAdjByRank; /* A DMLabel containing all points adjacent to shared points, separated by rank (value in label) */ 470 PetscSF sfPoint; 471 const PetscSFNode *remote; 472 const PetscInt *local; 473 const PetscInt *nrank, *rrank; 474 PetscInt *adj = NULL; 475 PetscInt pStart, pEnd, p, sStart, sEnd, nleaves, l; 476 PetscMPIInt rank, size; 477 PetscBool flg; 478 479 PetscFunctionBegin; 480 *ovLabel = NULL; 481 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 482 CHKERRMPI(MPI_Comm_size(comm, &size)); 483 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 484 if (size == 1) PetscFunctionReturn(0); 485 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 486 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 487 if (!levels) { 488 /* Add owned points */ 489 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel)); 490 for (p = pStart; p < pEnd; ++p) CHKERRQ(DMLabelSetValue(*ovLabel, p, rank)); 491 PetscFunctionReturn(0); 492 } 493 CHKERRQ(PetscSectionGetChart(leafSection, &sStart, &sEnd)); 494 CHKERRQ(PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote)); 495 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Overlap adjacency", &ovAdjByRank)); 496 /* Handle leaves: shared with the root point */ 497 for (l = 0; l < nleaves; ++l) { 498 PetscInt adjSize = PETSC_DETERMINE, a; 499 500 CHKERRQ(DMPlexGetAdjacency(dm, local ? local[l] : l, &adjSize, &adj)); 501 for (a = 0; a < adjSize; ++a) CHKERRQ(DMLabelSetValue(ovAdjByRank, adj[a], remote[l].rank)); 502 } 503 CHKERRQ(ISGetIndices(rootrank, &rrank)); 504 CHKERRQ(ISGetIndices(leafrank, &nrank)); 505 /* Handle roots */ 506 for (p = pStart; p < pEnd; ++p) { 507 PetscInt adjSize = PETSC_DETERMINE, neighbors = 0, noff, n, a; 508 509 if ((p >= sStart) && (p < sEnd)) { 510 /* Some leaves share a root with other leaves on different processes */ 511 CHKERRQ(PetscSectionGetDof(leafSection, p, &neighbors)); 512 if (neighbors) { 513 CHKERRQ(PetscSectionGetOffset(leafSection, p, &noff)); 514 CHKERRQ(DMPlexGetAdjacency(dm, p, &adjSize, &adj)); 515 for (n = 0; n < neighbors; ++n) { 516 const PetscInt remoteRank = nrank[noff+n]; 517 518 if (remoteRank == rank) continue; 519 for (a = 0; a < adjSize; ++a) CHKERRQ(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank)); 520 } 521 } 522 } 523 /* Roots are shared with leaves */ 524 CHKERRQ(PetscSectionGetDof(rootSection, p, &neighbors)); 525 if (!neighbors) continue; 526 CHKERRQ(PetscSectionGetOffset(rootSection, p, &noff)); 527 CHKERRQ(DMPlexGetAdjacency(dm, p, &adjSize, &adj)); 528 for (n = 0; n < neighbors; ++n) { 529 const PetscInt remoteRank = rrank[noff+n]; 530 531 if (remoteRank == rank) continue; 532 for (a = 0; a < adjSize; ++a) CHKERRQ(DMLabelSetValue(ovAdjByRank, adj[a], remoteRank)); 533 } 534 } 535 CHKERRQ(PetscFree(adj)); 536 CHKERRQ(ISRestoreIndices(rootrank, &rrank)); 537 CHKERRQ(ISRestoreIndices(leafrank, &nrank)); 538 /* Add additional overlap levels */ 539 for (l = 1; l < levels; l++) { 540 /* Propagate point donations over SF to capture remote connections */ 541 CHKERRQ(DMPlexPartitionLabelPropagate(dm, ovAdjByRank)); 542 /* Add next level of point donations to the label */ 543 CHKERRQ(DMPlexPartitionLabelAdjacency(dm, ovAdjByRank)); 544 } 545 /* We require the closure in the overlap */ 546 CHKERRQ(DMPlexPartitionLabelClosure(dm, ovAdjByRank)); 547 CHKERRQ(PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-overlap_view", &flg)); 548 if (flg) { 549 PetscViewer viewer; 550 CHKERRQ(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &viewer)); 551 CHKERRQ(DMLabelView(ovAdjByRank, viewer)); 552 } 553 /* Invert sender to receiver label */ 554 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Overlap label", ovLabel)); 555 CHKERRQ(DMPlexPartitionLabelInvert(dm, ovAdjByRank, NULL, *ovLabel)); 556 /* Add owned points, except for shared local points */ 557 for (p = pStart; p < pEnd; ++p) CHKERRQ(DMLabelSetValue(*ovLabel, p, rank)); 558 for (l = 0; l < nleaves; ++l) { 559 CHKERRQ(DMLabelClearValue(*ovLabel, local[l], rank)); 560 CHKERRQ(DMLabelSetValue(*ovLabel, remote[l].index, remote[l].rank)); 561 } 562 /* Clean up */ 563 CHKERRQ(DMLabelDestroy(&ovAdjByRank)); 564 PetscFunctionReturn(0); 565 } 566 567 /*@C 568 DMPlexCreateOverlapMigrationSF - Create an SF describing the new mesh distribution to make the overlap described by the input SF 569 570 Collective on dm 571 572 Input Parameters: 573 + dm - The DM 574 - overlapSF - The SF mapping ghost points in overlap to owner points on other processes 575 576 Output Parameter: 577 . migrationSF - An SF that maps original points in old locations to points in new locations 578 579 Level: developer 580 581 .seealso: DMPlexCreateOverlapLabel(), DMPlexDistributeOverlap(), DMPlexDistribute() 582 @*/ 583 PetscErrorCode DMPlexCreateOverlapMigrationSF(DM dm, PetscSF overlapSF, PetscSF *migrationSF) 584 { 585 MPI_Comm comm; 586 PetscMPIInt rank, size; 587 PetscInt d, dim, p, pStart, pEnd, nroots, nleaves, newLeaves, point, numSharedPoints; 588 PetscInt *pointDepths, *remoteDepths, *ilocal; 589 PetscInt *depthRecv, *depthShift, *depthIdx; 590 PetscSFNode *iremote; 591 PetscSF pointSF; 592 const PetscInt *sharedLocal; 593 const PetscSFNode *overlapRemote, *sharedRemote; 594 595 PetscFunctionBegin; 596 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 597 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 598 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 599 CHKERRMPI(MPI_Comm_size(comm, &size)); 600 CHKERRQ(DMGetDimension(dm, &dim)); 601 602 /* Before building the migration SF we need to know the new stratum offsets */ 603 CHKERRQ(PetscSFGetGraph(overlapSF, &nroots, &nleaves, NULL, &overlapRemote)); 604 CHKERRQ(PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths)); 605 for (d=0; d<dim+1; d++) { 606 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 607 for (p=pStart; p<pEnd; p++) pointDepths[p] = d; 608 } 609 for (p=0; p<nleaves; p++) remoteDepths[p] = -1; 610 CHKERRQ(PetscSFBcastBegin(overlapSF, MPIU_INT, pointDepths, remoteDepths,MPI_REPLACE)); 611 CHKERRQ(PetscSFBcastEnd(overlapSF, MPIU_INT, pointDepths, remoteDepths,MPI_REPLACE)); 612 613 /* Count received points in each stratum and compute the internal strata shift */ 614 CHKERRQ(PetscMalloc3(dim+1, &depthRecv, dim+1, &depthShift, dim+1, &depthIdx)); 615 for (d=0; d<dim+1; d++) depthRecv[d]=0; 616 for (p=0; p<nleaves; p++) depthRecv[remoteDepths[p]]++; 617 depthShift[dim] = 0; 618 for (d=0; d<dim; d++) depthShift[d] = depthRecv[dim]; 619 for (d=1; d<dim; d++) depthShift[d] += depthRecv[0]; 620 for (d=dim-2; d>0; d--) depthShift[d] += depthRecv[d+1]; 621 for (d=0; d<dim+1; d++) { 622 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 623 depthIdx[d] = pStart + depthShift[d]; 624 } 625 626 /* Form the overlap SF build an SF that describes the full overlap migration SF */ 627 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 628 newLeaves = pEnd - pStart + nleaves; 629 CHKERRQ(PetscMalloc1(newLeaves, &ilocal)); 630 CHKERRQ(PetscMalloc1(newLeaves, &iremote)); 631 /* First map local points to themselves */ 632 for (d=0; d<dim+1; d++) { 633 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 634 for (p=pStart; p<pEnd; p++) { 635 point = p + depthShift[d]; 636 ilocal[point] = point; 637 iremote[point].index = p; 638 iremote[point].rank = rank; 639 depthIdx[d]++; 640 } 641 } 642 643 /* Add in the remote roots for currently shared points */ 644 CHKERRQ(DMGetPointSF(dm, &pointSF)); 645 CHKERRQ(PetscSFGetGraph(pointSF, NULL, &numSharedPoints, &sharedLocal, &sharedRemote)); 646 for (d=0; d<dim+1; d++) { 647 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 648 for (p=0; p<numSharedPoints; p++) { 649 if (pStart <= sharedLocal[p] && sharedLocal[p] < pEnd) { 650 point = sharedLocal[p] + depthShift[d]; 651 iremote[point].index = sharedRemote[p].index; 652 iremote[point].rank = sharedRemote[p].rank; 653 } 654 } 655 } 656 657 /* Now add the incoming overlap points */ 658 for (p=0; p<nleaves; p++) { 659 point = depthIdx[remoteDepths[p]]; 660 ilocal[point] = point; 661 iremote[point].index = overlapRemote[p].index; 662 iremote[point].rank = overlapRemote[p].rank; 663 depthIdx[remoteDepths[p]]++; 664 } 665 CHKERRQ(PetscFree2(pointDepths,remoteDepths)); 666 667 CHKERRQ(PetscSFCreate(comm, migrationSF)); 668 CHKERRQ(PetscObjectSetName((PetscObject) *migrationSF, "Overlap Migration SF")); 669 CHKERRQ(PetscSFSetFromOptions(*migrationSF)); 670 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 671 CHKERRQ(PetscSFSetGraph(*migrationSF, pEnd-pStart, newLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 672 673 CHKERRQ(PetscFree3(depthRecv, depthShift, depthIdx)); 674 PetscFunctionReturn(0); 675 } 676 677 /*@ 678 DMPlexStratifyMigrationSF - Rearrange the leaves of a migration sf for stratification. 679 680 Input Parameters: 681 + dm - The DM 682 - sf - A star forest with non-ordered leaves, usually defining a DM point migration 683 684 Output Parameter: 685 . migrationSF - A star forest with added leaf indirection that ensures the resulting DM is stratified 686 687 Note: 688 This lexicographically sorts by (depth, cellType) 689 690 Level: developer 691 692 .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexDistributeOverlap() 693 @*/ 694 PetscErrorCode DMPlexStratifyMigrationSF(DM dm, PetscSF sf, PetscSF *migrationSF) 695 { 696 MPI_Comm comm; 697 PetscMPIInt rank, size; 698 PetscInt d, ldepth, depth, dim, p, pStart, pEnd, nroots, nleaves; 699 PetscSFNode *pointDepths, *remoteDepths; 700 PetscInt *ilocal; 701 PetscInt *depthRecv, *depthShift, *depthIdx; 702 PetscInt *ctRecv, *ctShift, *ctIdx; 703 const PetscSFNode *iremote; 704 705 PetscFunctionBegin; 706 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 707 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 708 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 709 CHKERRMPI(MPI_Comm_size(comm, &size)); 710 CHKERRQ(DMPlexGetDepth(dm, &ldepth)); 711 CHKERRQ(DMGetDimension(dm, &dim)); 712 CHKERRMPI(MPIU_Allreduce(&ldepth, &depth, 1, MPIU_INT, MPI_MAX, comm)); 713 PetscCheckFalse((ldepth >= 0) && (depth != ldepth),PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", ldepth, depth); 714 CHKERRQ(PetscLogEventBegin(DMPLEX_PartStratSF,dm,0,0,0)); 715 716 /* Before building the migration SF we need to know the new stratum offsets */ 717 CHKERRQ(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, &iremote)); 718 CHKERRQ(PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths)); 719 for (d = 0; d < depth+1; ++d) { 720 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 721 for (p = pStart; p < pEnd; ++p) { 722 DMPolytopeType ct; 723 724 CHKERRQ(DMPlexGetCellType(dm, p, &ct)); 725 pointDepths[p].index = d; 726 pointDepths[p].rank = ct; 727 } 728 } 729 for (p = 0; p < nleaves; ++p) {remoteDepths[p].index = -1; remoteDepths[p].rank = -1;} 730 CHKERRQ(PetscSFBcastBegin(sf, MPIU_2INT, pointDepths, remoteDepths,MPI_REPLACE)); 731 CHKERRQ(PetscSFBcastEnd(sf, MPIU_2INT, pointDepths, remoteDepths,MPI_REPLACE)); 732 /* Count received points in each stratum and compute the internal strata shift */ 733 CHKERRQ(PetscCalloc6(depth+1, &depthRecv, depth+1, &depthShift, depth+1, &depthIdx, DM_NUM_POLYTOPES, &ctRecv, DM_NUM_POLYTOPES, &ctShift, DM_NUM_POLYTOPES, &ctIdx)); 734 for (p = 0; p < nleaves; ++p) { 735 if (remoteDepths[p].rank < 0) { 736 ++depthRecv[remoteDepths[p].index]; 737 } else { 738 ++ctRecv[remoteDepths[p].rank]; 739 } 740 } 741 { 742 PetscInt depths[4], dims[4], shift = 0, i, c; 743 744 /* Cells (depth), Vertices (0), Faces (depth-1), Edges (1) 745 Consider DM_POLYTOPE_FV_GHOST and DM_POLYTOPE_INTERIOR_GHOST as cells 746 */ 747 depths[0] = depth; depths[1] = 0; depths[2] = depth-1; depths[3] = 1; 748 dims[0] = dim; dims[1] = 0; dims[2] = dim-1; dims[3] = 1; 749 for (i = 0; i <= depth; ++i) { 750 const PetscInt dep = depths[i]; 751 const PetscInt dim = dims[i]; 752 753 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 754 if (DMPolytopeTypeGetDim((DMPolytopeType) c) != dim && !(i == 0 && (c == DM_POLYTOPE_FV_GHOST || c == DM_POLYTOPE_INTERIOR_GHOST))) continue; 755 ctShift[c] = shift; 756 shift += ctRecv[c]; 757 } 758 depthShift[dep] = shift; 759 shift += depthRecv[dep]; 760 } 761 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 762 const PetscInt ctDim = DMPolytopeTypeGetDim((DMPolytopeType) c); 763 764 if ((ctDim < 0 || ctDim > dim) && (c != DM_POLYTOPE_FV_GHOST && c != DM_POLYTOPE_INTERIOR_GHOST)) { 765 ctShift[c] = shift; 766 shift += ctRecv[c]; 767 } 768 } 769 } 770 /* Derive a new local permutation based on stratified indices */ 771 CHKERRQ(PetscMalloc1(nleaves, &ilocal)); 772 for (p = 0; p < nleaves; ++p) { 773 const PetscInt dep = remoteDepths[p].index; 774 const DMPolytopeType ct = (DMPolytopeType) remoteDepths[p].rank; 775 776 if ((PetscInt) ct < 0) { 777 ilocal[p] = depthShift[dep] + depthIdx[dep]; 778 ++depthIdx[dep]; 779 } else { 780 ilocal[p] = ctShift[ct] + ctIdx[ct]; 781 ++ctIdx[ct]; 782 } 783 } 784 CHKERRQ(PetscSFCreate(comm, migrationSF)); 785 CHKERRQ(PetscObjectSetName((PetscObject) *migrationSF, "Migration SF")); 786 CHKERRQ(PetscSFSetGraph(*migrationSF, nroots, nleaves, ilocal, PETSC_OWN_POINTER, (PetscSFNode*)iremote, PETSC_COPY_VALUES)); 787 CHKERRQ(PetscFree2(pointDepths,remoteDepths)); 788 CHKERRQ(PetscFree6(depthRecv, depthShift, depthIdx, ctRecv, ctShift, ctIdx)); 789 CHKERRQ(PetscLogEventEnd(DMPLEX_PartStratSF,dm,0,0,0)); 790 PetscFunctionReturn(0); 791 } 792 793 /*@ 794 DMPlexDistributeField - Distribute field data to match a given PetscSF, usually the SF from mesh distribution 795 796 Collective on dm 797 798 Input Parameters: 799 + dm - The DMPlex object 800 . pointSF - The PetscSF describing the communication pattern 801 . originalSection - The PetscSection for existing data layout 802 - originalVec - The existing data in a local vector 803 804 Output Parameters: 805 + newSection - The PetscSF describing the new data layout 806 - newVec - The new data in a local vector 807 808 Level: developer 809 810 .seealso: DMPlexDistribute(), DMPlexDistributeFieldIS(), DMPlexDistributeData() 811 @*/ 812 PetscErrorCode DMPlexDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec) 813 { 814 PetscSF fieldSF; 815 PetscInt *remoteOffsets, fieldSize; 816 PetscScalar *originalValues, *newValues; 817 818 PetscFunctionBegin; 819 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0)); 820 CHKERRQ(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection)); 821 822 CHKERRQ(PetscSectionGetStorageSize(newSection, &fieldSize)); 823 CHKERRQ(VecSetSizes(newVec, fieldSize, PETSC_DETERMINE)); 824 CHKERRQ(VecSetType(newVec,dm->vectype)); 825 826 CHKERRQ(VecGetArray(originalVec, &originalValues)); 827 CHKERRQ(VecGetArray(newVec, &newValues)); 828 CHKERRQ(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF)); 829 CHKERRQ(PetscFree(remoteOffsets)); 830 CHKERRQ(PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues,MPI_REPLACE)); 831 CHKERRQ(PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues,MPI_REPLACE)); 832 CHKERRQ(PetscSFDestroy(&fieldSF)); 833 CHKERRQ(VecRestoreArray(newVec, &newValues)); 834 CHKERRQ(VecRestoreArray(originalVec, &originalValues)); 835 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0)); 836 PetscFunctionReturn(0); 837 } 838 839 /*@ 840 DMPlexDistributeFieldIS - Distribute field data to match a given PetscSF, usually the SF from mesh distribution 841 842 Collective on dm 843 844 Input Parameters: 845 + dm - The DMPlex object 846 . pointSF - The PetscSF describing the communication pattern 847 . originalSection - The PetscSection for existing data layout 848 - originalIS - The existing data 849 850 Output Parameters: 851 + newSection - The PetscSF describing the new data layout 852 - newIS - The new data 853 854 Level: developer 855 856 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexDistributeData() 857 @*/ 858 PetscErrorCode DMPlexDistributeFieldIS(DM dm, PetscSF pointSF, PetscSection originalSection, IS originalIS, PetscSection newSection, IS *newIS) 859 { 860 PetscSF fieldSF; 861 PetscInt *newValues, *remoteOffsets, fieldSize; 862 const PetscInt *originalValues; 863 864 PetscFunctionBegin; 865 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0)); 866 CHKERRQ(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection)); 867 868 CHKERRQ(PetscSectionGetStorageSize(newSection, &fieldSize)); 869 CHKERRQ(PetscMalloc1(fieldSize, &newValues)); 870 871 CHKERRQ(ISGetIndices(originalIS, &originalValues)); 872 CHKERRQ(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF)); 873 CHKERRQ(PetscFree(remoteOffsets)); 874 CHKERRQ(PetscSFBcastBegin(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues,MPI_REPLACE)); 875 CHKERRQ(PetscSFBcastEnd(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues,MPI_REPLACE)); 876 CHKERRQ(PetscSFDestroy(&fieldSF)); 877 CHKERRQ(ISRestoreIndices(originalIS, &originalValues)); 878 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject) pointSF), fieldSize, newValues, PETSC_OWN_POINTER, newIS)); 879 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0)); 880 PetscFunctionReturn(0); 881 } 882 883 /*@ 884 DMPlexDistributeData - Distribute field data to match a given PetscSF, usually the SF from mesh distribution 885 886 Collective on dm 887 888 Input Parameters: 889 + dm - The DMPlex object 890 . pointSF - The PetscSF describing the communication pattern 891 . originalSection - The PetscSection for existing data layout 892 . datatype - The type of data 893 - originalData - The existing data 894 895 Output Parameters: 896 + newSection - The PetscSection describing the new data layout 897 - newData - The new data 898 899 Level: developer 900 901 .seealso: DMPlexDistribute(), DMPlexDistributeField() 902 @*/ 903 PetscErrorCode DMPlexDistributeData(DM dm, PetscSF pointSF, PetscSection originalSection, MPI_Datatype datatype, void *originalData, PetscSection newSection, void **newData) 904 { 905 PetscSF fieldSF; 906 PetscInt *remoteOffsets, fieldSize; 907 PetscMPIInt dataSize; 908 909 PetscFunctionBegin; 910 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeData,dm,0,0,0)); 911 CHKERRQ(PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection)); 912 913 CHKERRQ(PetscSectionGetStorageSize(newSection, &fieldSize)); 914 CHKERRMPI(MPI_Type_size(datatype, &dataSize)); 915 CHKERRQ(PetscMalloc(fieldSize * dataSize, newData)); 916 917 CHKERRQ(PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF)); 918 CHKERRQ(PetscFree(remoteOffsets)); 919 CHKERRQ(PetscSFBcastBegin(fieldSF, datatype, originalData, *newData,MPI_REPLACE)); 920 CHKERRQ(PetscSFBcastEnd(fieldSF, datatype, originalData, *newData,MPI_REPLACE)); 921 CHKERRQ(PetscSFDestroy(&fieldSF)); 922 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeData,dm,0,0,0)); 923 PetscFunctionReturn(0); 924 } 925 926 static PetscErrorCode DMPlexDistributeCones(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel) 927 { 928 DM_Plex *pmesh = (DM_Plex*) (dmParallel)->data; 929 MPI_Comm comm; 930 PetscSF coneSF; 931 PetscSection originalConeSection, newConeSection; 932 PetscInt *remoteOffsets, *cones, *globCones, *newCones, newConesSize; 933 PetscBool flg; 934 935 PetscFunctionBegin; 936 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 937 PetscValidHeaderSpecific(dmParallel, DM_CLASSID, 5); 938 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeCones,dm,0,0,0)); 939 /* Distribute cone section */ 940 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 941 CHKERRQ(DMPlexGetConeSection(dm, &originalConeSection)); 942 CHKERRQ(DMPlexGetConeSection(dmParallel, &newConeSection)); 943 CHKERRQ(PetscSFDistributeSection(migrationSF, originalConeSection, &remoteOffsets, newConeSection)); 944 CHKERRQ(DMSetUp(dmParallel)); 945 { 946 PetscInt pStart, pEnd, p; 947 948 CHKERRQ(PetscSectionGetChart(newConeSection, &pStart, &pEnd)); 949 for (p = pStart; p < pEnd; ++p) { 950 PetscInt coneSize; 951 CHKERRQ(PetscSectionGetDof(newConeSection, p, &coneSize)); 952 pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize); 953 } 954 } 955 /* Communicate and renumber cones */ 956 CHKERRQ(PetscSFCreateSectionSF(migrationSF, originalConeSection, remoteOffsets, newConeSection, &coneSF)); 957 CHKERRQ(PetscFree(remoteOffsets)); 958 CHKERRQ(DMPlexGetCones(dm, &cones)); 959 if (original) { 960 PetscInt numCones; 961 962 CHKERRQ(PetscSectionGetStorageSize(originalConeSection,&numCones)); 963 CHKERRQ(PetscMalloc1(numCones,&globCones)); 964 CHKERRQ(ISLocalToGlobalMappingApplyBlock(original, numCones, cones, globCones)); 965 } else { 966 globCones = cones; 967 } 968 CHKERRQ(DMPlexGetCones(dmParallel, &newCones)); 969 CHKERRQ(PetscSFBcastBegin(coneSF, MPIU_INT, globCones, newCones,MPI_REPLACE)); 970 CHKERRQ(PetscSFBcastEnd(coneSF, MPIU_INT, globCones, newCones,MPI_REPLACE)); 971 if (original) { 972 CHKERRQ(PetscFree(globCones)); 973 } 974 CHKERRQ(PetscSectionGetStorageSize(newConeSection, &newConesSize)); 975 CHKERRQ(ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, newConesSize, newCones, NULL, newCones)); 976 if (PetscDefined(USE_DEBUG)) { 977 PetscInt p; 978 PetscBool valid = PETSC_TRUE; 979 for (p = 0; p < newConesSize; ++p) { 980 if (newCones[p] < 0) {valid = PETSC_FALSE; CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "[%d] Point %D not in overlap SF\n", PetscGlobalRank,p));} 981 } 982 PetscCheckFalse(!valid,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map"); 983 } 984 CHKERRQ(PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-cones_view", &flg)); 985 if (flg) { 986 CHKERRQ(PetscPrintf(comm, "Serial Cone Section:\n")); 987 CHKERRQ(PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_(comm))); 988 CHKERRQ(PetscPrintf(comm, "Parallel Cone Section:\n")); 989 CHKERRQ(PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_(comm))); 990 CHKERRQ(PetscSFView(coneSF, NULL)); 991 } 992 CHKERRQ(DMPlexGetConeOrientations(dm, &cones)); 993 CHKERRQ(DMPlexGetConeOrientations(dmParallel, &newCones)); 994 CHKERRQ(PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones,MPI_REPLACE)); 995 CHKERRQ(PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones,MPI_REPLACE)); 996 CHKERRQ(PetscSFDestroy(&coneSF)); 997 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeCones,dm,0,0,0)); 998 /* Create supports and stratify DMPlex */ 999 { 1000 PetscInt pStart, pEnd; 1001 1002 CHKERRQ(PetscSectionGetChart(pmesh->coneSection, &pStart, &pEnd)); 1003 CHKERRQ(PetscSectionSetChart(pmesh->supportSection, pStart, pEnd)); 1004 } 1005 CHKERRQ(DMPlexSymmetrize(dmParallel)); 1006 CHKERRQ(DMPlexStratify(dmParallel)); 1007 { 1008 PetscBool useCone, useClosure, useAnchors; 1009 1010 CHKERRQ(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 1011 CHKERRQ(DMSetBasicAdjacency(dmParallel, useCone, useClosure)); 1012 CHKERRQ(DMPlexGetAdjacencyUseAnchors(dm, &useAnchors)); 1013 CHKERRQ(DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors)); 1014 } 1015 PetscFunctionReturn(0); 1016 } 1017 1018 static PetscErrorCode DMPlexDistributeCoordinates(DM dm, PetscSF migrationSF, DM dmParallel) 1019 { 1020 MPI_Comm comm; 1021 DM cdm, cdmParallel; 1022 PetscSection originalCoordSection, newCoordSection; 1023 Vec originalCoordinates, newCoordinates; 1024 PetscInt bs; 1025 PetscBool isper; 1026 const char *name; 1027 const PetscReal *maxCell, *L; 1028 const DMBoundaryType *bd; 1029 1030 PetscFunctionBegin; 1031 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1032 PetscValidHeaderSpecific(dmParallel, DM_CLASSID, 3); 1033 1034 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1035 CHKERRQ(DMGetCoordinateSection(dm, &originalCoordSection)); 1036 CHKERRQ(DMGetCoordinateSection(dmParallel, &newCoordSection)); 1037 CHKERRQ(DMGetCoordinatesLocal(dm, &originalCoordinates)); 1038 if (originalCoordinates) { 1039 CHKERRQ(VecCreate(PETSC_COMM_SELF, &newCoordinates)); 1040 CHKERRQ(PetscObjectGetName((PetscObject) originalCoordinates, &name)); 1041 CHKERRQ(PetscObjectSetName((PetscObject) newCoordinates, name)); 1042 1043 CHKERRQ(DMPlexDistributeField(dm, migrationSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates)); 1044 CHKERRQ(DMSetCoordinatesLocal(dmParallel, newCoordinates)); 1045 CHKERRQ(VecGetBlockSize(originalCoordinates, &bs)); 1046 CHKERRQ(VecSetBlockSize(newCoordinates, bs)); 1047 CHKERRQ(VecDestroy(&newCoordinates)); 1048 } 1049 CHKERRQ(DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd)); 1050 CHKERRQ(DMSetPeriodicity(dmParallel, isper, maxCell, L, bd)); 1051 CHKERRQ(DMGetCoordinateDM(dm, &cdm)); 1052 CHKERRQ(DMGetCoordinateDM(dmParallel, &cdmParallel)); 1053 CHKERRQ(DMCopyDisc(cdm, cdmParallel)); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 static PetscErrorCode DMPlexDistributeLabels(DM dm, PetscSF migrationSF, DM dmParallel) 1058 { 1059 DM_Plex *mesh = (DM_Plex*) dm->data; 1060 MPI_Comm comm; 1061 DMLabel depthLabel; 1062 PetscMPIInt rank; 1063 PetscInt depth, d, numLabels, numLocalLabels, l; 1064 PetscBool hasLabels = PETSC_FALSE, lsendDepth, sendDepth; 1065 PetscObjectState depthState = -1; 1066 1067 PetscFunctionBegin; 1068 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1069 PetscValidHeaderSpecific(dmParallel, DM_CLASSID, 3); 1070 1071 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeLabels,dm,0,0,0)); 1072 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1073 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1074 1075 /* If the user has changed the depth label, communicate it instead */ 1076 CHKERRQ(DMPlexGetDepth(dm, &depth)); 1077 CHKERRQ(DMPlexGetDepthLabel(dm, &depthLabel)); 1078 if (depthLabel) CHKERRQ(PetscObjectStateGet((PetscObject) depthLabel, &depthState)); 1079 lsendDepth = mesh->depthState != depthState ? PETSC_TRUE : PETSC_FALSE; 1080 CHKERRMPI(MPIU_Allreduce(&lsendDepth, &sendDepth, 1, MPIU_BOOL, MPI_LOR, comm)); 1081 if (sendDepth) { 1082 CHKERRQ(DMPlexGetDepthLabel(dmParallel, &dmParallel->depthLabel)); 1083 CHKERRQ(DMRemoveLabelBySelf(dmParallel, &dmParallel->depthLabel, PETSC_FALSE)); 1084 } 1085 /* Everyone must have either the same number of labels, or none */ 1086 CHKERRQ(DMGetNumLabels(dm, &numLocalLabels)); 1087 numLabels = numLocalLabels; 1088 CHKERRMPI(MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm)); 1089 if (numLabels == numLocalLabels) hasLabels = PETSC_TRUE; 1090 for (l = 0; l < numLabels; ++l) { 1091 DMLabel label = NULL, labelNew = NULL; 1092 PetscBool isDepth, lisOutput = PETSC_TRUE, isOutput; 1093 const char *name = NULL; 1094 1095 if (hasLabels) { 1096 CHKERRQ(DMGetLabelByNum(dm, l, &label)); 1097 /* Skip "depth" because it is recreated */ 1098 CHKERRQ(PetscObjectGetName((PetscObject) label, &name)); 1099 CHKERRQ(PetscStrcmp(name, "depth", &isDepth)); 1100 } 1101 CHKERRMPI(MPI_Bcast(&isDepth, 1, MPIU_BOOL, 0, comm)); 1102 if (isDepth && !sendDepth) continue; 1103 CHKERRQ(DMLabelDistribute(label, migrationSF, &labelNew)); 1104 if (isDepth) { 1105 /* Put in any missing strata which can occur if users are managing the depth label themselves */ 1106 PetscInt gdepth; 1107 1108 CHKERRMPI(MPIU_Allreduce(&depth, &gdepth, 1, MPIU_INT, MPI_MAX, comm)); 1109 PetscCheckFalse((depth >= 0) && (gdepth != depth),PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", depth, gdepth); 1110 for (d = 0; d <= gdepth; ++d) { 1111 PetscBool has; 1112 1113 CHKERRQ(DMLabelHasStratum(labelNew, d, &has)); 1114 if (!has) CHKERRQ(DMLabelAddStratum(labelNew, d)); 1115 } 1116 } 1117 CHKERRQ(DMAddLabel(dmParallel, labelNew)); 1118 /* Put the output flag in the new label */ 1119 if (hasLabels) CHKERRQ(DMGetLabelOutput(dm, name, &lisOutput)); 1120 CHKERRMPI(MPIU_Allreduce(&lisOutput, &isOutput, 1, MPIU_BOOL, MPI_LAND, comm)); 1121 CHKERRQ(PetscObjectGetName((PetscObject) labelNew, &name)); 1122 CHKERRQ(DMSetLabelOutput(dmParallel, name, isOutput)); 1123 CHKERRQ(DMLabelDestroy(&labelNew)); 1124 } 1125 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeLabels,dm,0,0,0)); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 static PetscErrorCode DMPlexDistributeSetupTree(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel) 1130 { 1131 DM_Plex *mesh = (DM_Plex*) dm->data; 1132 DM_Plex *pmesh = (DM_Plex*) (dmParallel)->data; 1133 MPI_Comm comm; 1134 DM refTree; 1135 PetscSection origParentSection, newParentSection; 1136 PetscInt *origParents, *origChildIDs; 1137 PetscBool flg; 1138 1139 PetscFunctionBegin; 1140 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1141 PetscValidHeaderSpecific(dmParallel, DM_CLASSID, 5); 1142 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1143 1144 /* Set up tree */ 1145 CHKERRQ(DMPlexGetReferenceTree(dm,&refTree)); 1146 CHKERRQ(DMPlexSetReferenceTree(dmParallel,refTree)); 1147 CHKERRQ(DMPlexGetTree(dm,&origParentSection,&origParents,&origChildIDs,NULL,NULL)); 1148 if (origParentSection) { 1149 PetscInt pStart, pEnd; 1150 PetscInt *newParents, *newChildIDs, *globParents; 1151 PetscInt *remoteOffsetsParents, newParentSize; 1152 PetscSF parentSF; 1153 1154 CHKERRQ(DMPlexGetChart(dmParallel, &pStart, &pEnd)); 1155 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject)dmParallel),&newParentSection)); 1156 CHKERRQ(PetscSectionSetChart(newParentSection,pStart,pEnd)); 1157 CHKERRQ(PetscSFDistributeSection(migrationSF, origParentSection, &remoteOffsetsParents, newParentSection)); 1158 CHKERRQ(PetscSFCreateSectionSF(migrationSF, origParentSection, remoteOffsetsParents, newParentSection, &parentSF)); 1159 CHKERRQ(PetscFree(remoteOffsetsParents)); 1160 CHKERRQ(PetscSectionGetStorageSize(newParentSection,&newParentSize)); 1161 CHKERRQ(PetscMalloc2(newParentSize,&newParents,newParentSize,&newChildIDs)); 1162 if (original) { 1163 PetscInt numParents; 1164 1165 CHKERRQ(PetscSectionGetStorageSize(origParentSection,&numParents)); 1166 CHKERRQ(PetscMalloc1(numParents,&globParents)); 1167 CHKERRQ(ISLocalToGlobalMappingApplyBlock(original, numParents, origParents, globParents)); 1168 } 1169 else { 1170 globParents = origParents; 1171 } 1172 CHKERRQ(PetscSFBcastBegin(parentSF, MPIU_INT, globParents, newParents,MPI_REPLACE)); 1173 CHKERRQ(PetscSFBcastEnd(parentSF, MPIU_INT, globParents, newParents,MPI_REPLACE)); 1174 if (original) { 1175 CHKERRQ(PetscFree(globParents)); 1176 } 1177 CHKERRQ(PetscSFBcastBegin(parentSF, MPIU_INT, origChildIDs, newChildIDs,MPI_REPLACE)); 1178 CHKERRQ(PetscSFBcastEnd(parentSF, MPIU_INT, origChildIDs, newChildIDs,MPI_REPLACE)); 1179 CHKERRQ(ISGlobalToLocalMappingApplyBlock(renumbering,IS_GTOLM_MASK, newParentSize, newParents, NULL, newParents)); 1180 if (PetscDefined(USE_DEBUG)) { 1181 PetscInt p; 1182 PetscBool valid = PETSC_TRUE; 1183 for (p = 0; p < newParentSize; ++p) { 1184 if (newParents[p] < 0) {valid = PETSC_FALSE; CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "Point %d not in overlap SF\n", p));} 1185 } 1186 PetscCheckFalse(!valid,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map"); 1187 } 1188 CHKERRQ(PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-parents_view", &flg)); 1189 if (flg) { 1190 CHKERRQ(PetscPrintf(comm, "Serial Parent Section: \n")); 1191 CHKERRQ(PetscSectionView(origParentSection, PETSC_VIEWER_STDOUT_(comm))); 1192 CHKERRQ(PetscPrintf(comm, "Parallel Parent Section: \n")); 1193 CHKERRQ(PetscSectionView(newParentSection, PETSC_VIEWER_STDOUT_(comm))); 1194 CHKERRQ(PetscSFView(parentSF, NULL)); 1195 } 1196 CHKERRQ(DMPlexSetTree(dmParallel,newParentSection,newParents,newChildIDs)); 1197 CHKERRQ(PetscSectionDestroy(&newParentSection)); 1198 CHKERRQ(PetscFree2(newParents,newChildIDs)); 1199 CHKERRQ(PetscSFDestroy(&parentSF)); 1200 } 1201 pmesh->useAnchors = mesh->useAnchors; 1202 PetscFunctionReturn(0); 1203 } 1204 1205 PETSC_UNUSED static PetscErrorCode DMPlexDistributeSF(DM dm, PetscSF migrationSF, DM dmParallel) 1206 { 1207 PetscMPIInt rank, size; 1208 MPI_Comm comm; 1209 1210 PetscFunctionBegin; 1211 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1212 PetscValidHeaderSpecific(dmParallel, DM_CLASSID, 3); 1213 1214 /* Create point SF for parallel mesh */ 1215 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeSF,dm,0,0,0)); 1216 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1217 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1218 CHKERRMPI(MPI_Comm_size(comm, &size)); 1219 { 1220 const PetscInt *leaves; 1221 PetscSFNode *remotePoints, *rowners, *lowners; 1222 PetscInt numRoots, numLeaves, numGhostPoints = 0, p, gp, *ghostPoints; 1223 PetscInt pStart, pEnd; 1224 1225 CHKERRQ(DMPlexGetChart(dmParallel, &pStart, &pEnd)); 1226 CHKERRQ(PetscSFGetGraph(migrationSF, &numRoots, &numLeaves, &leaves, NULL)); 1227 CHKERRQ(PetscMalloc2(numRoots,&rowners,numLeaves,&lowners)); 1228 for (p=0; p<numRoots; p++) { 1229 rowners[p].rank = -1; 1230 rowners[p].index = -1; 1231 } 1232 CHKERRQ(PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE)); 1233 CHKERRQ(PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE)); 1234 for (p = 0; p < numLeaves; ++p) { 1235 if (lowners[p].rank < 0 || lowners[p].rank == rank) { /* Either put in a bid or we know we own it */ 1236 lowners[p].rank = rank; 1237 lowners[p].index = leaves ? leaves[p] : p; 1238 } else if (lowners[p].rank >= 0) { /* Point already claimed so flag so that MAXLOC does not listen to us */ 1239 lowners[p].rank = -2; 1240 lowners[p].index = -2; 1241 } 1242 } 1243 for (p=0; p<numRoots; p++) { /* Root must not participate in the rediction, flag so that MAXLOC does not use */ 1244 rowners[p].rank = -3; 1245 rowners[p].index = -3; 1246 } 1247 CHKERRQ(PetscSFReduceBegin(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC)); 1248 CHKERRQ(PetscSFReduceEnd(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC)); 1249 CHKERRQ(PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE)); 1250 CHKERRQ(PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners,MPI_REPLACE)); 1251 for (p = 0; p < numLeaves; ++p) { 1252 PetscCheckFalse(lowners[p].rank < 0 || lowners[p].index < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cell partition corrupt: point not claimed"); 1253 if (lowners[p].rank != rank) ++numGhostPoints; 1254 } 1255 CHKERRQ(PetscMalloc1(numGhostPoints, &ghostPoints)); 1256 CHKERRQ(PetscMalloc1(numGhostPoints, &remotePoints)); 1257 for (p = 0, gp = 0; p < numLeaves; ++p) { 1258 if (lowners[p].rank != rank) { 1259 ghostPoints[gp] = leaves ? leaves[p] : p; 1260 remotePoints[gp].rank = lowners[p].rank; 1261 remotePoints[gp].index = lowners[p].index; 1262 ++gp; 1263 } 1264 } 1265 CHKERRQ(PetscFree2(rowners,lowners)); 1266 CHKERRQ(PetscSFSetGraph((dmParallel)->sf, pEnd - pStart, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER)); 1267 CHKERRQ(PetscSFSetFromOptions((dmParallel)->sf)); 1268 } 1269 { 1270 PetscBool useCone, useClosure, useAnchors; 1271 1272 CHKERRQ(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 1273 CHKERRQ(DMSetBasicAdjacency(dmParallel, useCone, useClosure)); 1274 CHKERRQ(DMPlexGetAdjacencyUseAnchors(dm, &useAnchors)); 1275 CHKERRQ(DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors)); 1276 } 1277 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeSF,dm,0,0,0)); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 /*@ 1282 DMPlexSetPartitionBalance - Should distribution of the DM attempt to balance the shared point partition? 1283 1284 Input Parameters: 1285 + dm - The DMPlex object 1286 - flg - Balance the partition? 1287 1288 Level: intermediate 1289 1290 .seealso: DMPlexDistribute(), DMPlexGetPartitionBalance() 1291 @*/ 1292 PetscErrorCode DMPlexSetPartitionBalance(DM dm, PetscBool flg) 1293 { 1294 DM_Plex *mesh = (DM_Plex *)dm->data; 1295 1296 PetscFunctionBegin; 1297 mesh->partitionBalance = flg; 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /*@ 1302 DMPlexGetPartitionBalance - Does distribution of the DM attempt to balance the shared point partition? 1303 1304 Input Parameter: 1305 . dm - The DMPlex object 1306 1307 Output Parameter: 1308 . flg - Balance the partition? 1309 1310 Level: intermediate 1311 1312 .seealso: DMPlexDistribute(), DMPlexSetPartitionBalance() 1313 @*/ 1314 PetscErrorCode DMPlexGetPartitionBalance(DM dm, PetscBool *flg) 1315 { 1316 DM_Plex *mesh = (DM_Plex *)dm->data; 1317 1318 PetscFunctionBegin; 1319 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1320 PetscValidBoolPointer(flg, 2); 1321 *flg = mesh->partitionBalance; 1322 PetscFunctionReturn(0); 1323 } 1324 1325 typedef struct { 1326 PetscInt vote, rank, index; 1327 } Petsc3Int; 1328 1329 /* MaxLoc, but carry a third piece of information around */ 1330 static void MPIAPI MaxLocCarry(void *in_, void *inout_, PetscMPIInt *len_, MPI_Datatype *dtype) 1331 { 1332 Petsc3Int *a = (Petsc3Int *)inout_; 1333 Petsc3Int *b = (Petsc3Int *)in_; 1334 PetscInt i, len = *len_; 1335 for (i = 0; i < len; i++) { 1336 if (a[i].vote < b[i].vote) { 1337 a[i].vote = b[i].vote; 1338 a[i].rank = b[i].rank; 1339 a[i].index = b[i].index; 1340 } else if (a[i].vote <= b[i].vote) { 1341 if (a[i].rank >= b[i].rank) { 1342 a[i].rank = b[i].rank; 1343 a[i].index = b[i].index; 1344 } 1345 } 1346 } 1347 } 1348 1349 /*@C 1350 DMPlexCreatePointSF - Build a point SF from an SF describing a point migration 1351 1352 Input Parameters: 1353 + dm - The source DMPlex object 1354 . migrationSF - The star forest that describes the parallel point remapping 1355 - ownership - Flag causing a vote to determine point ownership 1356 1357 Output Parameter: 1358 . pointSF - The star forest describing the point overlap in the remapped DM 1359 1360 Notes: 1361 Output pointSF is guaranteed to have the array of local indices (leaves) sorted. 1362 1363 Level: developer 1364 1365 .seealso: DMPlexDistribute(), DMPlexDistributeOverlap() 1366 @*/ 1367 PetscErrorCode DMPlexCreatePointSF(DM dm, PetscSF migrationSF, PetscBool ownership, PetscSF *pointSF) 1368 { 1369 PetscMPIInt rank, size; 1370 PetscInt p, nroots, nleaves, idx, npointLeaves; 1371 PetscInt *pointLocal; 1372 const PetscInt *leaves; 1373 const PetscSFNode *roots; 1374 PetscSFNode *rootNodes, *leafNodes, *pointRemote; 1375 Vec shifts; 1376 const PetscInt numShifts = 13759; 1377 const PetscScalar *shift = NULL; 1378 const PetscBool shiftDebug = PETSC_FALSE; 1379 PetscBool balance; 1380 1381 PetscFunctionBegin; 1382 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1383 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 1384 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 1385 CHKERRQ(PetscLogEventBegin(DMPLEX_CreatePointSF,dm,0,0,0)); 1386 1387 CHKERRQ(DMPlexGetPartitionBalance(dm, &balance)); 1388 CHKERRQ(PetscSFGetGraph(migrationSF, &nroots, &nleaves, &leaves, &roots)); 1389 CHKERRQ(PetscMalloc2(nroots, &rootNodes, nleaves, &leafNodes)); 1390 if (ownership) { 1391 MPI_Op op; 1392 MPI_Datatype datatype; 1393 Petsc3Int *rootVote = NULL, *leafVote = NULL; 1394 /* If balancing, we compute a random cyclic shift of the rank for each remote point. That way, the max will evenly distribute among ranks. */ 1395 if (balance) { 1396 PetscRandom r; 1397 1398 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF, &r)); 1399 CHKERRQ(PetscRandomSetInterval(r, 0, 2467*size)); 1400 CHKERRQ(VecCreate(PETSC_COMM_SELF, &shifts)); 1401 CHKERRQ(VecSetSizes(shifts, numShifts, numShifts)); 1402 CHKERRQ(VecSetType(shifts, VECSTANDARD)); 1403 CHKERRQ(VecSetRandom(shifts, r)); 1404 CHKERRQ(PetscRandomDestroy(&r)); 1405 CHKERRQ(VecGetArrayRead(shifts, &shift)); 1406 } 1407 1408 CHKERRQ(PetscMalloc1(nroots, &rootVote)); 1409 CHKERRQ(PetscMalloc1(nleaves, &leafVote)); 1410 /* Point ownership vote: Process with highest rank owns shared points */ 1411 for (p = 0; p < nleaves; ++p) { 1412 if (shiftDebug) { 1413 CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Point %D RemotePoint %D Shift %D MyRank %D\n", rank, leaves ? leaves[p] : p, roots[p].index, (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]), (rank + (shift ? (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]) : 0))%size)); 1414 } 1415 /* Either put in a bid or we know we own it */ 1416 leafVote[p].vote = (rank + (shift ? (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]) : 0))%size; 1417 leafVote[p].rank = rank; 1418 leafVote[p].index = p; 1419 } 1420 for (p = 0; p < nroots; p++) { 1421 /* Root must not participate in the reduction, flag so that MAXLOC does not use */ 1422 rootVote[p].vote = -3; 1423 rootVote[p].rank = -3; 1424 rootVote[p].index = -3; 1425 } 1426 CHKERRMPI(MPI_Type_contiguous(3, MPIU_INT, &datatype)); 1427 CHKERRMPI(MPI_Type_commit(&datatype)); 1428 CHKERRMPI(MPI_Op_create(&MaxLocCarry, 1, &op)); 1429 CHKERRQ(PetscSFReduceBegin(migrationSF, datatype, leafVote, rootVote, op)); 1430 CHKERRQ(PetscSFReduceEnd(migrationSF, datatype, leafVote, rootVote, op)); 1431 CHKERRMPI(MPI_Op_free(&op)); 1432 CHKERRMPI(MPI_Type_free(&datatype)); 1433 for (p = 0; p < nroots; p++) { 1434 rootNodes[p].rank = rootVote[p].rank; 1435 rootNodes[p].index = rootVote[p].index; 1436 } 1437 CHKERRQ(PetscFree(leafVote)); 1438 CHKERRQ(PetscFree(rootVote)); 1439 } else { 1440 for (p = 0; p < nroots; p++) { 1441 rootNodes[p].index = -1; 1442 rootNodes[p].rank = rank; 1443 } 1444 for (p = 0; p < nleaves; p++) { 1445 /* Write new local id into old location */ 1446 if (roots[p].rank == rank) { 1447 rootNodes[roots[p].index].index = leaves ? leaves[p] : p; 1448 } 1449 } 1450 } 1451 CHKERRQ(PetscSFBcastBegin(migrationSF, MPIU_2INT, rootNodes, leafNodes,MPI_REPLACE)); 1452 CHKERRQ(PetscSFBcastEnd(migrationSF, MPIU_2INT, rootNodes, leafNodes,MPI_REPLACE)); 1453 1454 for (npointLeaves = 0, p = 0; p < nleaves; p++) { 1455 if (leafNodes[p].rank != rank) npointLeaves++; 1456 } 1457 CHKERRQ(PetscMalloc1(npointLeaves, &pointLocal)); 1458 CHKERRQ(PetscMalloc1(npointLeaves, &pointRemote)); 1459 for (idx = 0, p = 0; p < nleaves; p++) { 1460 if (leafNodes[p].rank != rank) { 1461 /* Note that pointLocal is automatically sorted as it is sublist of 0, ..., nleaves-1 */ 1462 pointLocal[idx] = p; 1463 pointRemote[idx] = leafNodes[p]; 1464 idx++; 1465 } 1466 } 1467 if (shift) { 1468 CHKERRQ(VecRestoreArrayRead(shifts, &shift)); 1469 CHKERRQ(VecDestroy(&shifts)); 1470 } 1471 if (shiftDebug) CHKERRQ(PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), PETSC_STDOUT)); 1472 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject) dm), pointSF)); 1473 CHKERRQ(PetscSFSetFromOptions(*pointSF)); 1474 CHKERRQ(PetscSFSetGraph(*pointSF, nleaves, npointLeaves, pointLocal, PETSC_OWN_POINTER, pointRemote, PETSC_OWN_POINTER)); 1475 CHKERRQ(PetscFree2(rootNodes, leafNodes)); 1476 CHKERRQ(PetscLogEventEnd(DMPLEX_CreatePointSF,dm,0,0,0)); 1477 PetscFunctionReturn(0); 1478 } 1479 1480 /*@C 1481 DMPlexMigrate - Migrates internal DM data over the supplied star forest 1482 1483 Collective on dm 1484 1485 Input Parameters: 1486 + dm - The source DMPlex object 1487 - sf - The star forest communication context describing the migration pattern 1488 1489 Output Parameter: 1490 . targetDM - The target DMPlex object 1491 1492 Level: intermediate 1493 1494 .seealso: DMPlexDistribute(), DMPlexDistributeOverlap() 1495 @*/ 1496 PetscErrorCode DMPlexMigrate(DM dm, PetscSF sf, DM targetDM) 1497 { 1498 MPI_Comm comm; 1499 PetscInt dim, cdim, nroots; 1500 PetscSF sfPoint; 1501 ISLocalToGlobalMapping ltogMigration; 1502 ISLocalToGlobalMapping ltogOriginal = NULL; 1503 PetscBool flg; 1504 1505 PetscFunctionBegin; 1506 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1507 CHKERRQ(PetscLogEventBegin(DMPLEX_Migrate, dm, 0, 0, 0)); 1508 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 1509 CHKERRQ(DMGetDimension(dm, &dim)); 1510 CHKERRQ(DMSetDimension(targetDM, dim)); 1511 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 1512 CHKERRQ(DMSetCoordinateDim(targetDM, cdim)); 1513 1514 /* Check for a one-to-all distribution pattern */ 1515 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 1516 CHKERRQ(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1517 if (nroots >= 0) { 1518 IS isOriginal; 1519 PetscInt n, size, nleaves; 1520 PetscInt *numbering_orig, *numbering_new; 1521 1522 /* Get the original point numbering */ 1523 CHKERRQ(DMPlexCreatePointNumbering(dm, &isOriginal)); 1524 CHKERRQ(ISLocalToGlobalMappingCreateIS(isOriginal, <ogOriginal)); 1525 CHKERRQ(ISLocalToGlobalMappingGetSize(ltogOriginal, &size)); 1526 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(ltogOriginal, (const PetscInt**)&numbering_orig)); 1527 /* Convert to positive global numbers */ 1528 for (n=0; n<size; n++) {if (numbering_orig[n] < 0) numbering_orig[n] = -(numbering_orig[n]+1);} 1529 /* Derive the new local-to-global mapping from the old one */ 1530 CHKERRQ(PetscSFGetGraph(sf, NULL, &nleaves, NULL, NULL)); 1531 CHKERRQ(PetscMalloc1(nleaves, &numbering_new)); 1532 CHKERRQ(PetscSFBcastBegin(sf, MPIU_INT, numbering_orig, numbering_new,MPI_REPLACE)); 1533 CHKERRQ(PetscSFBcastEnd(sf, MPIU_INT, numbering_orig, numbering_new,MPI_REPLACE)); 1534 CHKERRQ(ISLocalToGlobalMappingCreate(comm, 1, nleaves, numbering_new, PETSC_OWN_POINTER, <ogMigration)); 1535 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(ltogOriginal, (const PetscInt**)&numbering_orig)); 1536 CHKERRQ(ISDestroy(&isOriginal)); 1537 } else { 1538 /* One-to-all distribution pattern: We can derive LToG from SF */ 1539 CHKERRQ(ISLocalToGlobalMappingCreateSF(sf, 0, <ogMigration)); 1540 } 1541 CHKERRQ(PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg)); 1542 if (flg) { 1543 CHKERRQ(PetscPrintf(comm, "Point renumbering for DM migration:\n")); 1544 CHKERRQ(ISLocalToGlobalMappingView(ltogMigration, NULL)); 1545 } 1546 /* Migrate DM data to target DM */ 1547 CHKERRQ(DMPlexDistributeCones(dm, sf, ltogOriginal, ltogMigration, targetDM)); 1548 CHKERRQ(DMPlexDistributeLabels(dm, sf, targetDM)); 1549 CHKERRQ(DMPlexDistributeCoordinates(dm, sf, targetDM)); 1550 CHKERRQ(DMPlexDistributeSetupTree(dm, sf, ltogOriginal, ltogMigration, targetDM)); 1551 CHKERRQ(ISLocalToGlobalMappingDestroy(<ogOriginal)); 1552 CHKERRQ(ISLocalToGlobalMappingDestroy(<ogMigration)); 1553 CHKERRQ(PetscLogEventEnd(DMPLEX_Migrate, dm, 0, 0, 0)); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 /*@C 1558 DMPlexDistribute - Distributes the mesh and any associated sections. 1559 1560 Collective on dm 1561 1562 Input Parameters: 1563 + dm - The original DMPlex object 1564 - overlap - The overlap of partitions, 0 is the default 1565 1566 Output Parameters: 1567 + sf - The PetscSF used for point distribution, or NULL if not needed 1568 - dmParallel - The distributed DMPlex object 1569 1570 Note: If the mesh was not distributed, the output dmParallel will be NULL. 1571 1572 The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function 1573 representation on the mesh. 1574 1575 Level: intermediate 1576 1577 .seealso: DMPlexCreate(), DMSetAdjacency(), DMPlexGetOverlap() 1578 @*/ 1579 PetscErrorCode DMPlexDistribute(DM dm, PetscInt overlap, PetscSF *sf, DM *dmParallel) 1580 { 1581 MPI_Comm comm; 1582 PetscPartitioner partitioner; 1583 IS cellPart; 1584 PetscSection cellPartSection; 1585 DM dmCoord; 1586 DMLabel lblPartition, lblMigration; 1587 PetscSF sfMigration, sfStratified, sfPoint; 1588 PetscBool flg, balance; 1589 PetscMPIInt rank, size; 1590 1591 PetscFunctionBegin; 1592 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1593 PetscValidLogicalCollectiveInt(dm, overlap, 2); 1594 if (sf) PetscValidPointer(sf,3); 1595 PetscValidPointer(dmParallel,4); 1596 1597 if (sf) *sf = NULL; 1598 *dmParallel = NULL; 1599 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 1600 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1601 CHKERRMPI(MPI_Comm_size(comm, &size)); 1602 if (size == 1) PetscFunctionReturn(0); 1603 1604 CHKERRQ(PetscLogEventBegin(DMPLEX_Distribute,dm,0,0,0)); 1605 /* Create cell partition */ 1606 CHKERRQ(PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0)); 1607 CHKERRQ(PetscSectionCreate(comm, &cellPartSection)); 1608 CHKERRQ(DMPlexGetPartitioner(dm, &partitioner)); 1609 CHKERRQ(PetscPartitionerDMPlexPartition(partitioner, dm, NULL, cellPartSection, &cellPart)); 1610 CHKERRQ(PetscLogEventBegin(DMPLEX_PartSelf,dm,0,0,0)); 1611 { 1612 /* Convert partition to DMLabel */ 1613 IS is; 1614 PetscHSetI ht; 1615 const PetscInt *points; 1616 PetscInt *iranks; 1617 PetscInt pStart, pEnd, proc, npoints, poff = 0, nranks; 1618 1619 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Point Partition", &lblPartition)); 1620 /* Preallocate strata */ 1621 CHKERRQ(PetscHSetICreate(&ht)); 1622 CHKERRQ(PetscSectionGetChart(cellPartSection, &pStart, &pEnd)); 1623 for (proc = pStart; proc < pEnd; proc++) { 1624 CHKERRQ(PetscSectionGetDof(cellPartSection, proc, &npoints)); 1625 if (npoints) CHKERRQ(PetscHSetIAdd(ht, proc)); 1626 } 1627 CHKERRQ(PetscHSetIGetSize(ht, &nranks)); 1628 CHKERRQ(PetscMalloc1(nranks, &iranks)); 1629 CHKERRQ(PetscHSetIGetElems(ht, &poff, iranks)); 1630 CHKERRQ(PetscHSetIDestroy(&ht)); 1631 CHKERRQ(DMLabelAddStrata(lblPartition, nranks, iranks)); 1632 CHKERRQ(PetscFree(iranks)); 1633 /* Inline DMPlexPartitionLabelClosure() */ 1634 CHKERRQ(ISGetIndices(cellPart, &points)); 1635 CHKERRQ(PetscSectionGetChart(cellPartSection, &pStart, &pEnd)); 1636 for (proc = pStart; proc < pEnd; proc++) { 1637 CHKERRQ(PetscSectionGetDof(cellPartSection, proc, &npoints)); 1638 if (!npoints) continue; 1639 CHKERRQ(PetscSectionGetOffset(cellPartSection, proc, &poff)); 1640 CHKERRQ(DMPlexClosurePoints_Private(dm, npoints, points+poff, &is)); 1641 CHKERRQ(DMLabelSetStratumIS(lblPartition, proc, is)); 1642 CHKERRQ(ISDestroy(&is)); 1643 } 1644 CHKERRQ(ISRestoreIndices(cellPart, &points)); 1645 } 1646 CHKERRQ(PetscLogEventEnd(DMPLEX_PartSelf,dm,0,0,0)); 1647 1648 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "Point migration", &lblMigration)); 1649 CHKERRQ(DMPlexPartitionLabelInvert(dm, lblPartition, NULL, lblMigration)); 1650 CHKERRQ(DMPlexPartitionLabelCreateSF(dm, lblMigration, &sfMigration)); 1651 CHKERRQ(DMPlexStratifyMigrationSF(dm, sfMigration, &sfStratified)); 1652 CHKERRQ(PetscSFDestroy(&sfMigration)); 1653 sfMigration = sfStratified; 1654 CHKERRQ(PetscSFSetUp(sfMigration)); 1655 CHKERRQ(PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0)); 1656 CHKERRQ(PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg)); 1657 if (flg) { 1658 CHKERRQ(DMLabelView(lblPartition, PETSC_VIEWER_STDOUT_(comm))); 1659 CHKERRQ(PetscSFView(sfMigration, PETSC_VIEWER_STDOUT_(comm))); 1660 } 1661 1662 /* Create non-overlapping parallel DM and migrate internal data */ 1663 CHKERRQ(DMPlexCreate(comm, dmParallel)); 1664 CHKERRQ(PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh")); 1665 CHKERRQ(DMPlexMigrate(dm, sfMigration, *dmParallel)); 1666 1667 /* Build the point SF without overlap */ 1668 CHKERRQ(DMPlexGetPartitionBalance(dm, &balance)); 1669 CHKERRQ(DMPlexSetPartitionBalance(*dmParallel, balance)); 1670 CHKERRQ(DMPlexCreatePointSF(*dmParallel, sfMigration, PETSC_TRUE, &sfPoint)); 1671 CHKERRQ(DMSetPointSF(*dmParallel, sfPoint)); 1672 CHKERRQ(DMGetCoordinateDM(*dmParallel, &dmCoord)); 1673 if (dmCoord) CHKERRQ(DMSetPointSF(dmCoord, sfPoint)); 1674 if (flg) CHKERRQ(PetscSFView(sfPoint, NULL)); 1675 1676 if (overlap > 0) { 1677 DM dmOverlap; 1678 PetscInt nroots, nleaves, noldleaves, l; 1679 const PetscInt *oldLeaves; 1680 PetscSFNode *newRemote, *permRemote; 1681 const PetscSFNode *oldRemote; 1682 PetscSF sfOverlap, sfOverlapPoint; 1683 1684 /* Add the partition overlap to the distributed DM */ 1685 CHKERRQ(DMPlexDistributeOverlap(*dmParallel, overlap, &sfOverlap, &dmOverlap)); 1686 CHKERRQ(DMDestroy(dmParallel)); 1687 *dmParallel = dmOverlap; 1688 if (flg) { 1689 CHKERRQ(PetscPrintf(comm, "Overlap Migration SF:\n")); 1690 CHKERRQ(PetscSFView(sfOverlap, NULL)); 1691 } 1692 1693 /* Re-map the migration SF to establish the full migration pattern */ 1694 CHKERRQ(PetscSFGetGraph(sfMigration, &nroots, &noldleaves, &oldLeaves, &oldRemote)); 1695 CHKERRQ(PetscSFGetGraph(sfOverlap, NULL, &nleaves, NULL, NULL)); 1696 CHKERRQ(PetscMalloc1(nleaves, &newRemote)); 1697 /* oldRemote: original root point mapping to original leaf point 1698 newRemote: original leaf point mapping to overlapped leaf point */ 1699 if (oldLeaves) { 1700 /* After stratification, the migration remotes may not be in root (canonical) order, so we reorder using the leaf numbering */ 1701 CHKERRQ(PetscMalloc1(noldleaves, &permRemote)); 1702 for (l = 0; l < noldleaves; ++l) permRemote[oldLeaves[l]] = oldRemote[l]; 1703 oldRemote = permRemote; 1704 } 1705 CHKERRQ(PetscSFBcastBegin(sfOverlap, MPIU_2INT, oldRemote, newRemote,MPI_REPLACE)); 1706 CHKERRQ(PetscSFBcastEnd(sfOverlap, MPIU_2INT, oldRemote, newRemote,MPI_REPLACE)); 1707 if (oldLeaves) CHKERRQ(PetscFree(oldRemote)); 1708 CHKERRQ(PetscSFCreate(comm, &sfOverlapPoint)); 1709 CHKERRQ(PetscSFSetGraph(sfOverlapPoint, nroots, nleaves, NULL, PETSC_OWN_POINTER, newRemote, PETSC_OWN_POINTER)); 1710 CHKERRQ(PetscSFDestroy(&sfOverlap)); 1711 CHKERRQ(PetscSFDestroy(&sfMigration)); 1712 sfMigration = sfOverlapPoint; 1713 } 1714 /* Cleanup Partition */ 1715 CHKERRQ(DMLabelDestroy(&lblPartition)); 1716 CHKERRQ(DMLabelDestroy(&lblMigration)); 1717 CHKERRQ(PetscSectionDestroy(&cellPartSection)); 1718 CHKERRQ(ISDestroy(&cellPart)); 1719 /* Copy discretization */ 1720 CHKERRQ(DMCopyDisc(dm, *dmParallel)); 1721 /* Create sfNatural */ 1722 if (dm->useNatural) { 1723 PetscSection section; 1724 1725 CHKERRQ(DMGetLocalSection(dm, §ion)); 1726 CHKERRQ(DMPlexCreateGlobalToNaturalSF(*dmParallel, section, sfMigration, &(*dmParallel)->sfNatural)); 1727 CHKERRQ(DMSetUseNatural(*dmParallel, PETSC_TRUE)); 1728 } 1729 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *dmParallel)); 1730 /* Cleanup */ 1731 if (sf) {*sf = sfMigration;} 1732 else CHKERRQ(PetscSFDestroy(&sfMigration)); 1733 CHKERRQ(PetscSFDestroy(&sfPoint)); 1734 CHKERRQ(PetscLogEventEnd(DMPLEX_Distribute,dm,0,0,0)); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 /*@C 1739 DMPlexDistributeOverlap - Add partition overlap to a distributed non-overlapping DM. 1740 1741 Collective on dm 1742 1743 Input Parameters: 1744 + dm - The non-overlapping distributed DMPlex object 1745 - overlap - The overlap of partitions (the same on all ranks) 1746 1747 Output Parameters: 1748 + sf - The PetscSF used for point distribution 1749 - dmOverlap - The overlapping distributed DMPlex object, or NULL 1750 1751 Notes: 1752 If the mesh was not distributed, the return value is NULL. 1753 1754 The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function 1755 representation on the mesh. 1756 1757 Level: advanced 1758 1759 .seealso: DMPlexCreate(), DMSetAdjacency(), DMPlexDistribute(), DMPlexCreateOverlapLabel(), DMPlexGetOverlap() 1760 @*/ 1761 PetscErrorCode DMPlexDistributeOverlap(DM dm, PetscInt overlap, PetscSF *sf, DM *dmOverlap) 1762 { 1763 MPI_Comm comm; 1764 PetscMPIInt size, rank; 1765 PetscSection rootSection, leafSection; 1766 IS rootrank, leafrank; 1767 DM dmCoord; 1768 DMLabel lblOverlap; 1769 PetscSF sfOverlap, sfStratified, sfPoint; 1770 1771 PetscFunctionBegin; 1772 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1773 PetscValidLogicalCollectiveInt(dm, overlap, 2); 1774 if (sf) PetscValidPointer(sf, 3); 1775 PetscValidPointer(dmOverlap, 4); 1776 1777 if (sf) *sf = NULL; 1778 *dmOverlap = NULL; 1779 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 1780 CHKERRMPI(MPI_Comm_size(comm, &size)); 1781 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1782 if (size == 1) PetscFunctionReturn(0); 1783 1784 CHKERRQ(PetscLogEventBegin(DMPLEX_DistributeOverlap, dm, 0, 0, 0)); 1785 /* Compute point overlap with neighbouring processes on the distributed DM */ 1786 CHKERRQ(PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0)); 1787 CHKERRQ(PetscSectionCreate(comm, &rootSection)); 1788 CHKERRQ(PetscSectionCreate(comm, &leafSection)); 1789 CHKERRQ(DMPlexDistributeOwnership(dm, rootSection, &rootrank, leafSection, &leafrank)); 1790 CHKERRQ(DMPlexCreateOverlapLabel(dm, overlap, rootSection, rootrank, leafSection, leafrank, &lblOverlap)); 1791 /* Convert overlap label to stratified migration SF */ 1792 CHKERRQ(DMPlexPartitionLabelCreateSF(dm, lblOverlap, &sfOverlap)); 1793 CHKERRQ(DMPlexStratifyMigrationSF(dm, sfOverlap, &sfStratified)); 1794 CHKERRQ(PetscSFDestroy(&sfOverlap)); 1795 sfOverlap = sfStratified; 1796 CHKERRQ(PetscObjectSetName((PetscObject) sfOverlap, "Overlap SF")); 1797 CHKERRQ(PetscSFSetFromOptions(sfOverlap)); 1798 1799 CHKERRQ(PetscSectionDestroy(&rootSection)); 1800 CHKERRQ(PetscSectionDestroy(&leafSection)); 1801 CHKERRQ(ISDestroy(&rootrank)); 1802 CHKERRQ(ISDestroy(&leafrank)); 1803 CHKERRQ(PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0)); 1804 1805 /* Build the overlapping DM */ 1806 CHKERRQ(DMPlexCreate(comm, dmOverlap)); 1807 CHKERRQ(PetscObjectSetName((PetscObject) *dmOverlap, "Parallel Mesh")); 1808 CHKERRQ(DMPlexMigrate(dm, sfOverlap, *dmOverlap)); 1809 /* Store the overlap in the new DM */ 1810 ((DM_Plex*)(*dmOverlap)->data)->overlap = overlap + ((DM_Plex*)dm->data)->overlap; 1811 /* Build the new point SF */ 1812 CHKERRQ(DMPlexCreatePointSF(*dmOverlap, sfOverlap, PETSC_FALSE, &sfPoint)); 1813 CHKERRQ(DMSetPointSF(*dmOverlap, sfPoint)); 1814 CHKERRQ(DMGetCoordinateDM(*dmOverlap, &dmCoord)); 1815 if (dmCoord) CHKERRQ(DMSetPointSF(dmCoord, sfPoint)); 1816 CHKERRQ(PetscSFDestroy(&sfPoint)); 1817 /* Cleanup overlap partition */ 1818 CHKERRQ(DMLabelDestroy(&lblOverlap)); 1819 if (sf) *sf = sfOverlap; 1820 else CHKERRQ(PetscSFDestroy(&sfOverlap)); 1821 CHKERRQ(PetscLogEventEnd(DMPLEX_DistributeOverlap, dm, 0, 0, 0)); 1822 PetscFunctionReturn(0); 1823 } 1824 1825 PetscErrorCode DMPlexGetOverlap_Plex(DM dm, PetscInt *overlap) 1826 { 1827 DM_Plex *mesh = (DM_Plex*) dm->data; 1828 1829 PetscFunctionBegin; 1830 *overlap = mesh->overlap; 1831 PetscFunctionReturn(0); 1832 } 1833 1834 /*@ 1835 DMPlexGetOverlap - Get the DMPlex partition overlap. 1836 1837 Not collective 1838 1839 Input Parameter: 1840 . dm - The DM 1841 1842 Output Parameter: 1843 . overlap - The overlap of this DM 1844 1845 Level: intermediate 1846 1847 .seealso: DMPlexDistribute(), DMPlexDistributeOverlap(), DMPlexCreateOverlapLabel() 1848 @*/ 1849 PetscErrorCode DMPlexGetOverlap(DM dm, PetscInt *overlap) 1850 { 1851 PetscFunctionBegin; 1852 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1853 CHKERRQ(PetscUseMethod(dm,"DMPlexGetOverlap_C",(DM,PetscInt*),(dm,overlap))); 1854 PetscFunctionReturn(0); 1855 } 1856 1857 PetscErrorCode DMPlexDistributeSetDefault_Plex(DM dm, PetscBool dist) 1858 { 1859 DM_Plex *mesh = (DM_Plex *) dm->data; 1860 1861 PetscFunctionBegin; 1862 mesh->distDefault = dist; 1863 PetscFunctionReturn(0); 1864 } 1865 1866 /*@ 1867 DMPlexDistributeSetDefault - Set flag indicating whether the DM should be distributed by default 1868 1869 Logically collective 1870 1871 Input Parameters: 1872 + dm - The DM 1873 - dist - Flag for distribution 1874 1875 Level: intermediate 1876 1877 .seealso: DMDistributeGetDefault(), DMPlexDistribute() 1878 @*/ 1879 PetscErrorCode DMPlexDistributeSetDefault(DM dm, PetscBool dist) 1880 { 1881 PetscFunctionBegin; 1882 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1883 PetscValidLogicalCollectiveBool(dm, dist, 2); 1884 CHKERRQ(PetscTryMethod(dm,"DMPlexDistributeSetDefault_C",(DM,PetscBool),(dm,dist))); 1885 PetscFunctionReturn(0); 1886 } 1887 1888 PetscErrorCode DMPlexDistributeGetDefault_Plex(DM dm, PetscBool *dist) 1889 { 1890 DM_Plex *mesh = (DM_Plex *) dm->data; 1891 1892 PetscFunctionBegin; 1893 *dist = mesh->distDefault; 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /*@ 1898 DMPlexDistributeGetDefault - Get flag indicating whether the DM should be distributed by default 1899 1900 Not collective 1901 1902 Input Parameter: 1903 . dm - The DM 1904 1905 Output Parameter: 1906 . dist - Flag for distribution 1907 1908 Level: intermediate 1909 1910 .seealso: DMDistributeSetDefault(), DMPlexDistribute() 1911 @*/ 1912 PetscErrorCode DMPlexDistributeGetDefault(DM dm, PetscBool *dist) 1913 { 1914 PetscFunctionBegin; 1915 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1916 PetscValidBoolPointer(dist, 2); 1917 CHKERRQ(PetscUseMethod(dm,"DMPlexDistributeGetDefault_C",(DM,PetscBool*),(dm,dist))); 1918 PetscFunctionReturn(0); 1919 } 1920 1921 /*@C 1922 DMPlexGetGatherDM - Get a copy of the DMPlex that gathers all points on the 1923 root process of the original's communicator. 1924 1925 Collective on dm 1926 1927 Input Parameter: 1928 . dm - the original DMPlex object 1929 1930 Output Parameters: 1931 + sf - the PetscSF used for point distribution (optional) 1932 - gatherMesh - the gathered DM object, or NULL 1933 1934 Level: intermediate 1935 1936 .seealso: DMPlexDistribute(), DMPlexGetRedundantDM() 1937 @*/ 1938 PetscErrorCode DMPlexGetGatherDM(DM dm, PetscSF *sf, DM *gatherMesh) 1939 { 1940 MPI_Comm comm; 1941 PetscMPIInt size; 1942 PetscPartitioner oldPart, gatherPart; 1943 1944 PetscFunctionBegin; 1945 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1946 PetscValidPointer(gatherMesh,3); 1947 *gatherMesh = NULL; 1948 if (sf) *sf = NULL; 1949 comm = PetscObjectComm((PetscObject)dm); 1950 CHKERRMPI(MPI_Comm_size(comm,&size)); 1951 if (size == 1) PetscFunctionReturn(0); 1952 CHKERRQ(DMPlexGetPartitioner(dm,&oldPart)); 1953 CHKERRQ(PetscObjectReference((PetscObject)oldPart)); 1954 CHKERRQ(PetscPartitionerCreate(comm,&gatherPart)); 1955 CHKERRQ(PetscPartitionerSetType(gatherPart,PETSCPARTITIONERGATHER)); 1956 CHKERRQ(DMPlexSetPartitioner(dm,gatherPart)); 1957 CHKERRQ(DMPlexDistribute(dm,0,sf,gatherMesh)); 1958 1959 CHKERRQ(DMPlexSetPartitioner(dm,oldPart)); 1960 CHKERRQ(PetscPartitionerDestroy(&gatherPart)); 1961 CHKERRQ(PetscPartitionerDestroy(&oldPart)); 1962 PetscFunctionReturn(0); 1963 } 1964 1965 /*@C 1966 DMPlexGetRedundantDM - Get a copy of the DMPlex that is completely copied on each process. 1967 1968 Collective on dm 1969 1970 Input Parameter: 1971 . dm - the original DMPlex object 1972 1973 Output Parameters: 1974 + sf - the PetscSF used for point distribution (optional) 1975 - redundantMesh - the redundant DM object, or NULL 1976 1977 Level: intermediate 1978 1979 .seealso: DMPlexDistribute(), DMPlexGetGatherDM() 1980 @*/ 1981 PetscErrorCode DMPlexGetRedundantDM(DM dm, PetscSF *sf, DM *redundantMesh) 1982 { 1983 MPI_Comm comm; 1984 PetscMPIInt size, rank; 1985 PetscInt pStart, pEnd, p; 1986 PetscInt numPoints = -1; 1987 PetscSF migrationSF, sfPoint, gatherSF; 1988 DM gatherDM, dmCoord; 1989 PetscSFNode *points; 1990 1991 PetscFunctionBegin; 1992 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1993 PetscValidPointer(redundantMesh,3); 1994 *redundantMesh = NULL; 1995 comm = PetscObjectComm((PetscObject)dm); 1996 CHKERRMPI(MPI_Comm_size(comm,&size)); 1997 if (size == 1) { 1998 CHKERRQ(PetscObjectReference((PetscObject) dm)); 1999 *redundantMesh = dm; 2000 if (sf) *sf = NULL; 2001 PetscFunctionReturn(0); 2002 } 2003 CHKERRQ(DMPlexGetGatherDM(dm,&gatherSF,&gatherDM)); 2004 if (!gatherDM) PetscFunctionReturn(0); 2005 CHKERRMPI(MPI_Comm_rank(comm,&rank)); 2006 CHKERRQ(DMPlexGetChart(gatherDM,&pStart,&pEnd)); 2007 numPoints = pEnd - pStart; 2008 CHKERRMPI(MPI_Bcast(&numPoints,1,MPIU_INT,0,comm)); 2009 CHKERRQ(PetscMalloc1(numPoints,&points)); 2010 CHKERRQ(PetscSFCreate(comm,&migrationSF)); 2011 for (p = 0; p < numPoints; p++) { 2012 points[p].index = p; 2013 points[p].rank = 0; 2014 } 2015 CHKERRQ(PetscSFSetGraph(migrationSF,pEnd-pStart,numPoints,NULL,PETSC_OWN_POINTER,points,PETSC_OWN_POINTER)); 2016 CHKERRQ(DMPlexCreate(comm, redundantMesh)); 2017 CHKERRQ(PetscObjectSetName((PetscObject) *redundantMesh, "Redundant Mesh")); 2018 CHKERRQ(DMPlexMigrate(gatherDM, migrationSF, *redundantMesh)); 2019 CHKERRQ(DMPlexCreatePointSF(*redundantMesh, migrationSF, PETSC_FALSE, &sfPoint)); 2020 CHKERRQ(DMSetPointSF(*redundantMesh, sfPoint)); 2021 CHKERRQ(DMGetCoordinateDM(*redundantMesh, &dmCoord)); 2022 if (dmCoord) CHKERRQ(DMSetPointSF(dmCoord, sfPoint)); 2023 CHKERRQ(PetscSFDestroy(&sfPoint)); 2024 if (sf) { 2025 PetscSF tsf; 2026 2027 CHKERRQ(PetscSFCompose(gatherSF,migrationSF,&tsf)); 2028 CHKERRQ(DMPlexStratifyMigrationSF(dm, tsf, sf)); 2029 CHKERRQ(PetscSFDestroy(&tsf)); 2030 } 2031 CHKERRQ(PetscSFDestroy(&migrationSF)); 2032 CHKERRQ(PetscSFDestroy(&gatherSF)); 2033 CHKERRQ(DMDestroy(&gatherDM)); 2034 CHKERRQ(DMCopyDisc(dm, *redundantMesh)); 2035 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *redundantMesh)); 2036 PetscFunctionReturn(0); 2037 } 2038 2039 /*@ 2040 DMPlexIsDistributed - Find out whether this DM is distributed, i.e. more than one rank owns some points. 2041 2042 Collective 2043 2044 Input Parameter: 2045 . dm - The DM object 2046 2047 Output Parameter: 2048 . distributed - Flag whether the DM is distributed 2049 2050 Level: intermediate 2051 2052 Notes: 2053 This currently finds out whether at least two ranks have any DAG points. 2054 This involves MPI_Allreduce() with one integer. 2055 The result is currently not stashed so every call to this routine involves this global communication. 2056 2057 .seealso: DMPlexDistribute(), DMPlexGetOverlap(), DMPlexIsInterpolated() 2058 @*/ 2059 PetscErrorCode DMPlexIsDistributed(DM dm, PetscBool *distributed) 2060 { 2061 PetscInt pStart, pEnd, count; 2062 MPI_Comm comm; 2063 PetscMPIInt size; 2064 2065 PetscFunctionBegin; 2066 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2067 PetscValidPointer(distributed,2); 2068 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 2069 CHKERRMPI(MPI_Comm_size(comm,&size)); 2070 if (size == 1) { *distributed = PETSC_FALSE; PetscFunctionReturn(0); } 2071 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 2072 count = (pEnd - pStart) > 0 ? 1 : 0; 2073 CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE, &count, 1, MPIU_INT, MPI_SUM, comm)); 2074 *distributed = count > 1 ? PETSC_TRUE : PETSC_FALSE; 2075 PetscFunctionReturn(0); 2076 } 2077