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