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, PetscBool useInitialGuess) 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]-pStart] = 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 ierr = PetscSortInt(counter, adjncy); CHKERRQ(ierr); 1513 1514 for(i=0; i<pEnd-pStart; i++) { 1515 if(toBalance[i]) { 1516 if(isNonExclusivelyOwned[i]) { 1517 adjncy[counter] = cumSumVertices[rank]; 1518 counter++; 1519 for(j=0; j<degree[i]; j++) { 1520 adjncy[counter] = cumSumVertices[locationsOfLeafs[cumSumDegrees[i]+j]]; 1521 counter++; 1522 } 1523 if(degree[i]>0) { 1524 ierr = PetscSortInt(degree[i]+1, &adjncy[counter-degree[i]-1]); CHKERRQ(ierr); 1525 } 1526 } 1527 } 1528 } 1529 1530 nparts = size; 1531 wgtflag = 2; 1532 1533 numflag = 0; 1534 ncon = 1; 1535 ncommonnodes = 2; 1536 real_t *tpwgts; 1537 ierr = PetscMalloc1(ncon * nparts, &tpwgts);CHKERRQ(ierr); 1538 for(i=0; i<ncon*nparts; i++) { 1539 tpwgts[i] = 1./(ncon*nparts); 1540 } 1541 ierr = PetscMalloc1(1, &ubvec);CHKERRQ(ierr); 1542 ubvec[0] = 1.01; 1543 ierr = PetscMalloc1(1, &options);CHKERRQ(ierr); 1544 options[0] = 0; 1545 1546 ierr = PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);CHKERRQ(ierr); 1547 if(useInitialGuess) { 1548 PetscStackPush("ParMETIS_V3_RefineKway"); 1549 ierr = ParMETIS_V3_RefineKway(cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm); 1550 PetscStackPop; 1551 } else { 1552 PetscStackPush("ParMETIS_V3_PartKway"); 1553 ierr = ParMETIS_V3_PartKway(cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm); 1554 PetscStackPop; 1555 } 1556 1557 ierr = PetscFree(ubvec);CHKERRQ(ierr); 1558 ierr = PetscFree(tpwgts);CHKERRQ(ierr); 1559 ierr = PetscFree(options);CHKERRQ(ierr); 1560 1561 /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */ 1562 1563 ierr = PetscMalloc1(size, &firstVertices);CHKERRQ(ierr); 1564 ierr = PetscMalloc1(size, &renumbering);CHKERRQ(ierr); 1565 firstVertices[rank] = part[0]; 1566 MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);CHKERRQ(ierr); 1567 for(i=0; i<size; i++) { 1568 renumbering[firstVertices[i]] = i; 1569 } 1570 for(i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) { 1571 part[i] = renumbering[part[i]]; 1572 } 1573 /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */ 1574 failed = (PetscInt)(part[0] != rank); 1575 MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm); 1576 1577 if(failedGlobal > 0) { 1578 PetscFunctionReturn(1); 1579 } 1580 1581 ierr = PetscFree(firstVertices);CHKERRQ(ierr); 1582 ierr = PetscFree(renumbering);CHKERRQ(ierr); 1583 ierr = PetscFree(xadj);CHKERRQ(ierr); 1584 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1585 1586 1587 /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */ 1588 PetscInt *newOwners, *newNumbers; 1589 ierr = PetscMalloc1(pEnd-pStart, &newOwners);CHKERRQ(ierr); 1590 ierr = PetscMalloc1(pEnd-pStart, &newNumbers);CHKERRQ(ierr); 1591 for(i=0; i<pEnd-pStart; i++) { 1592 newOwners[i] = -1; 1593 newNumbers[i] = -1; 1594 } 1595 counter = 1; 1596 for(i=0; i<pEnd-pStart; i++) { 1597 if(toBalance[i]) { 1598 if(isNonExclusivelyOwned[i]) { 1599 PetscInt oldNumber, newNumber, oldOwner, newOwner; 1600 oldNumber = i; 1601 oldOwner = rank; 1602 newOwner = part[counter]; 1603 if(oldOwner == newOwner) { 1604 newNumber = oldNumber; 1605 } else { 1606 for(j=0; j<degree[i]; j++) { 1607 if(locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) { 1608 newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j]; 1609 break; 1610 } 1611 } 1612 } 1613 newOwners[oldNumber] = newOwner; 1614 newNumbers[oldNumber] = newNumber; 1615 counter++; 1616 } 1617 } 1618 } 1619 1620 PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners); 1621 PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners); 1622 PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers); 1623 PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers); 1624 1625 /* Now count how many leafs we have on each processor. */ 1626 leafCounter=0; 1627 for(i=0; i<pEnd-pStart; i++) { 1628 if(toBalance[i]){ 1629 if(newOwners[i] >= 0 && newOwners[i] != rank) { 1630 leafCounter++; 1631 } 1632 } else { 1633 if(isLeaf[i]) { 1634 leafCounter++; 1635 } 1636 } 1637 } 1638 1639 ierr = PetscMalloc1(leafCounter, &leafsNew);CHKERRQ(ierr); 1640 ierr = PetscMalloc1(leafCounter, &leafLocationsNew);CHKERRQ(ierr); 1641 1642 counter = 0; 1643 leafCounter = 0; 1644 for(i=0; i<pEnd-pStart; i++) { 1645 if(toBalance[i]){ 1646 if(isLeaf[i]) { 1647 counter++; 1648 } 1649 if(newOwners[i] >= 0 && newOwners[i] != rank) { 1650 leafsNew[leafCounter] = i; 1651 leafLocationsNew[leafCounter].rank = newOwners[i]; 1652 leafLocationsNew[leafCounter].index = newNumbers[i]; 1653 leafCounter++; 1654 } 1655 } else { 1656 if(isLeaf[i]) { 1657 leafsNew[leafCounter] = i; 1658 leafLocationsNew[leafCounter].rank = iremote[counter].rank; 1659 leafLocationsNew[leafCounter].index = iremote[counter].index; 1660 leafCounter++; 1661 counter++; 1662 } 1663 } 1664 } 1665 1666 ierr = PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1667 1668 1669 ierr = PetscFree(toBalance);CHKERRQ(ierr); 1670 ierr = PetscFree(isLeaf);CHKERRQ(ierr); 1671 ierr = PetscFree(numLocalVerticesAllProcesses);CHKERRQ(ierr); 1672 ierr = PetscFree(cumSumVertices);CHKERRQ(ierr); 1673 ierr = PetscFree(locationsOfLeafs);CHKERRQ(ierr); 1674 ierr = PetscFree(rankOnLeafs);CHKERRQ(ierr); 1675 ierr = PetscFree(remoteLocalPointOfLeafs);CHKERRQ(ierr); 1676 ierr = PetscFree(points);CHKERRQ(ierr); 1677 ierr = PetscFree(leafGlobalNumbers);CHKERRQ(ierr); 1678 ierr = PetscFree(newOwners);CHKERRQ(ierr); 1679 ierr = PetscFree(newNumbers);CHKERRQ(ierr); 1680 PetscFunctionReturn(0); 1681 } 1682 1683 /*@ 1684 DMPlexSetPartitionBalance - Should distribution of the DM attempt to balance the shared point partition? 1685 1686 Input Parameters: 1687 + dm - The DMPlex object 1688 - flg - Balance the partition? 1689 1690 Level: intermediate 1691 1692 .seealso: DMPlexDistribute(), DMPlexGetPartitionBalance() 1693 @*/ 1694 PetscErrorCode DMPlexSetPartitionBalance(DM dm, PetscBool flg) 1695 { 1696 DM_Plex *mesh = (DM_Plex *)dm->data; 1697 1698 PetscFunctionBegin; 1699 mesh->partitionBalance = flg; 1700 PetscFunctionReturn(0); 1701 } 1702 1703 /*@ 1704 DMPlexGetPartitionBalance - Does distribution of the DM attempt to balance the shared point partition? 1705 1706 Input Parameter: 1707 + dm - The DMPlex object 1708 1709 Output Parameter: 1710 + flg - Balance the partition? 1711 1712 Level: intermediate 1713 1714 .seealso: DMPlexDistribute(), DMPlexSetPartitionBalance() 1715 @*/ 1716 PetscErrorCode DMPlexGetPartitionBalance(DM dm, PetscBool *flg) 1717 { 1718 DM_Plex *mesh = (DM_Plex *)dm->data; 1719 1720 PetscFunctionBegin; 1721 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1722 PetscValidIntPointer(flg, 2); 1723 *flg = mesh->partitionBalance; 1724 PetscFunctionReturn(0); 1725 } 1726 1727 /*@C 1728 DMPlexDerivePointSF - Build a point SF from an SF describing a point migration 1729 1730 Input Parameter: 1731 + dm - The source DMPlex object 1732 . migrationSF - The star forest that describes the parallel point remapping 1733 . ownership - Flag causing a vote to determine point ownership 1734 1735 Output Parameter: 1736 - pointSF - The star forest describing the point overlap in the remapped DM 1737 1738 Level: developer 1739 1740 .seealso: DMPlexDistribute(), DMPlexDistributeOverlap() 1741 @*/ 1742 PetscErrorCode DMPlexCreatePointSF(DM dm, PetscSF migrationSF, PetscBool ownership, PetscSF *pointSF) 1743 { 1744 PetscMPIInt rank, size; 1745 PetscInt p, nroots, nleaves, idx, npointLeaves; 1746 PetscInt *pointLocal; 1747 const PetscInt *leaves; 1748 const PetscSFNode *roots; 1749 PetscSFNode *rootNodes, *leafNodes, *pointRemote; 1750 Vec shifts; 1751 const PetscInt numShifts = 13759; 1752 const PetscScalar *shift = NULL; 1753 const PetscBool shiftDebug = PETSC_FALSE; 1754 PetscBool balance; 1755 PetscErrorCode ierr; 1756 1757 PetscFunctionBegin; 1758 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1759 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1760 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 1761 1762 ierr = DMPlexGetPartitionBalance(dm, &balance);CHKERRQ(ierr); 1763 ierr = PetscSFGetGraph(migrationSF, &nroots, &nleaves, &leaves, &roots);CHKERRQ(ierr); 1764 ierr = PetscMalloc2(nroots, &rootNodes, nleaves, &leafNodes);CHKERRQ(ierr); 1765 if (ownership) { 1766 /* If balancing, we compute a random cyclic shift of the rank for each remote point. That way, the max will evenly distribute among ranks. */ 1767 if (balance) { 1768 PetscRandom r; 1769 1770 ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr); 1771 ierr = PetscRandomSetInterval(r, 0, 2467*size);CHKERRQ(ierr); 1772 ierr = VecCreate(PETSC_COMM_SELF, &shifts);CHKERRQ(ierr); 1773 ierr = VecSetSizes(shifts, numShifts, numShifts);CHKERRQ(ierr); 1774 ierr = VecSetType(shifts, VECSTANDARD);CHKERRQ(ierr); 1775 ierr = VecSetRandom(shifts, r);CHKERRQ(ierr); 1776 ierr = PetscRandomDestroy(&r);CHKERRQ(ierr); 1777 ierr = VecGetArrayRead(shifts, &shift);CHKERRQ(ierr); 1778 } 1779 1780 /* Point ownership vote: Process with highest rank owns shared points */ 1781 for (p = 0; p < nleaves; ++p) { 1782 if (shiftDebug) { 1783 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); 1784 } 1785 /* Either put in a bid or we know we own it */ 1786 leafNodes[p].rank = (rank + (shift ? (PetscInt) PetscRealPart(shift[roots[p].index%numShifts]) : 0))%size; 1787 leafNodes[p].index = p; 1788 } 1789 for (p = 0; p < nroots; p++) { 1790 /* Root must not participate in the reduction, flag so that MAXLOC does not use */ 1791 rootNodes[p].rank = -3; 1792 rootNodes[p].index = -3; 1793 } 1794 ierr = PetscSFReduceBegin(migrationSF, MPIU_2INT, leafNodes, rootNodes, MPI_MAXLOC);CHKERRQ(ierr); 1795 ierr = PetscSFReduceEnd(migrationSF, MPIU_2INT, leafNodes, rootNodes, MPI_MAXLOC);CHKERRQ(ierr); 1796 if (balance) { 1797 /* We've voted, now we need to get the rank. When we're balancing the partition, the "rank" in rootNotes is not 1798 * the rank but rather (rank + random)%size. So we do another reduction, voting the same way, but sending the 1799 * rank instead of the index. */ 1800 PetscSFNode *rootRanks = NULL; 1801 ierr = PetscMalloc1(nroots, &rootRanks);CHKERRQ(ierr); 1802 for (p = 0; p < nroots; p++) { 1803 rootRanks[p].rank = -3; 1804 rootRanks[p].index = -3; 1805 } 1806 for (p = 0; p < nleaves; p++) leafNodes[p].index = rank; 1807 ierr = PetscSFReduceBegin(migrationSF, MPIU_2INT, leafNodes, rootRanks, MPI_MAXLOC);CHKERRQ(ierr); 1808 ierr = PetscSFReduceEnd(migrationSF, MPIU_2INT, leafNodes, rootRanks, MPI_MAXLOC);CHKERRQ(ierr); 1809 for (p = 0; p < nroots; p++) rootNodes[p].rank = rootRanks[p].index; 1810 ierr = PetscFree(rootRanks);CHKERRQ(ierr); 1811 } 1812 } else { 1813 for (p = 0; p < nroots; p++) { 1814 rootNodes[p].index = -1; 1815 rootNodes[p].rank = rank; 1816 }; 1817 for (p = 0; p < nleaves; p++) { 1818 /* Write new local id into old location */ 1819 if (roots[p].rank == rank) { 1820 rootNodes[roots[p].index].index = leaves ? leaves[p] : p; 1821 } 1822 } 1823 } 1824 ierr = PetscSFBcastBegin(migrationSF, MPIU_2INT, rootNodes, leafNodes);CHKERRQ(ierr); 1825 ierr = PetscSFBcastEnd(migrationSF, MPIU_2INT, rootNodes, leafNodes);CHKERRQ(ierr); 1826 1827 for (npointLeaves = 0, p = 0; p < nleaves; p++) { 1828 if (leafNodes[p].rank != rank) npointLeaves++; 1829 } 1830 ierr = PetscMalloc1(npointLeaves, &pointLocal);CHKERRQ(ierr); 1831 ierr = PetscMalloc1(npointLeaves, &pointRemote);CHKERRQ(ierr); 1832 for (idx = 0, p = 0; p < nleaves; p++) { 1833 if (leafNodes[p].rank != rank) { 1834 pointLocal[idx] = p; 1835 pointRemote[idx] = leafNodes[p]; 1836 idx++; 1837 } 1838 } 1839 if (shift) { 1840 ierr = VecRestoreArrayRead(shifts, &shift);CHKERRQ(ierr); 1841 ierr = VecDestroy(&shifts);CHKERRQ(ierr); 1842 } 1843 if (shiftDebug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), PETSC_STDOUT);CHKERRQ(ierr);} 1844 ierr = PetscSFCreate(PetscObjectComm((PetscObject) dm), pointSF);CHKERRQ(ierr); 1845 ierr = PetscSFSetFromOptions(*pointSF);CHKERRQ(ierr); 1846 ierr = PetscSFSetGraph(*pointSF, nleaves, npointLeaves, pointLocal, PETSC_OWN_POINTER, pointRemote, PETSC_OWN_POINTER);CHKERRQ(ierr); 1847 ierr = PetscFree2(rootNodes, leafNodes);CHKERRQ(ierr); 1848 PetscFunctionReturn(0); 1849 } 1850 1851 /*@C 1852 DMPlexMigrate - Migrates internal DM data over the supplied star forest 1853 1854 Collective on DM and PetscSF 1855 1856 Input Parameter: 1857 + dm - The source DMPlex object 1858 . sf - The star forest communication context describing the migration pattern 1859 1860 Output Parameter: 1861 - targetDM - The target DMPlex object 1862 1863 Level: intermediate 1864 1865 .seealso: DMPlexDistribute(), DMPlexDistributeOverlap() 1866 @*/ 1867 PetscErrorCode DMPlexMigrate(DM dm, PetscSF sf, DM targetDM) 1868 { 1869 MPI_Comm comm; 1870 PetscInt dim, cdim, nroots; 1871 PetscSF sfPoint; 1872 ISLocalToGlobalMapping ltogMigration; 1873 ISLocalToGlobalMapping ltogOriginal = NULL; 1874 PetscBool flg; 1875 PetscErrorCode ierr; 1876 1877 PetscFunctionBegin; 1878 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1879 ierr = PetscLogEventBegin(DMPLEX_Migrate, dm, 0, 0, 0);CHKERRQ(ierr); 1880 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 1881 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1882 ierr = DMSetDimension(targetDM, dim);CHKERRQ(ierr); 1883 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1884 ierr = DMSetCoordinateDim(targetDM, cdim);CHKERRQ(ierr); 1885 1886 /* Check for a one-to-all distribution pattern */ 1887 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1888 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1889 if (nroots >= 0) { 1890 IS isOriginal; 1891 PetscInt n, size, nleaves; 1892 PetscInt *numbering_orig, *numbering_new; 1893 1894 /* Get the original point numbering */ 1895 ierr = DMPlexCreatePointNumbering(dm, &isOriginal);CHKERRQ(ierr); 1896 ierr = ISLocalToGlobalMappingCreateIS(isOriginal, <ogOriginal);CHKERRQ(ierr); 1897 ierr = ISLocalToGlobalMappingGetSize(ltogOriginal, &size);CHKERRQ(ierr); 1898 ierr = ISLocalToGlobalMappingGetBlockIndices(ltogOriginal, (const PetscInt**)&numbering_orig);CHKERRQ(ierr); 1899 /* Convert to positive global numbers */ 1900 for (n=0; n<size; n++) {if (numbering_orig[n] < 0) numbering_orig[n] = -(numbering_orig[n]+1);} 1901 /* Derive the new local-to-global mapping from the old one */ 1902 ierr = PetscSFGetGraph(sf, NULL, &nleaves, NULL, NULL);CHKERRQ(ierr); 1903 ierr = PetscMalloc1(nleaves, &numbering_new);CHKERRQ(ierr); 1904 ierr = PetscSFBcastBegin(sf, MPIU_INT, (PetscInt *) numbering_orig, numbering_new);CHKERRQ(ierr); 1905 ierr = PetscSFBcastEnd(sf, MPIU_INT, (PetscInt *) numbering_orig, numbering_new);CHKERRQ(ierr); 1906 ierr = ISLocalToGlobalMappingCreate(comm, 1, nleaves, (const PetscInt*) numbering_new, PETSC_OWN_POINTER, <ogMigration);CHKERRQ(ierr); 1907 ierr = ISLocalToGlobalMappingRestoreIndices(ltogOriginal, (const PetscInt**)&numbering_orig);CHKERRQ(ierr); 1908 ierr = ISDestroy(&isOriginal);CHKERRQ(ierr); 1909 } else { 1910 /* One-to-all distribution pattern: We can derive LToG from SF */ 1911 ierr = ISLocalToGlobalMappingCreateSF(sf, 0, <ogMigration);CHKERRQ(ierr); 1912 } 1913 ierr = PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);CHKERRQ(ierr); 1914 if (flg) { 1915 ierr = PetscPrintf(comm, "Point renumbering for DM migration:\n");CHKERRQ(ierr); 1916 ierr = ISLocalToGlobalMappingView(ltogMigration, NULL);CHKERRQ(ierr); 1917 } 1918 /* Migrate DM data to target DM */ 1919 ierr = DMPlexDistributeCones(dm, sf, ltogOriginal, ltogMigration, targetDM);CHKERRQ(ierr); 1920 ierr = DMPlexDistributeLabels(dm, sf, targetDM);CHKERRQ(ierr); 1921 ierr = DMPlexDistributeCoordinates(dm, sf, targetDM);CHKERRQ(ierr); 1922 ierr = DMPlexDistributeSetupHybrid(dm, sf, ltogMigration, targetDM);CHKERRQ(ierr); 1923 ierr = DMPlexDistributeSetupTree(dm, sf, ltogOriginal, ltogMigration, targetDM);CHKERRQ(ierr); 1924 ierr = ISLocalToGlobalMappingDestroy(<ogOriginal);CHKERRQ(ierr); 1925 ierr = ISLocalToGlobalMappingDestroy(<ogMigration);CHKERRQ(ierr); 1926 ierr = PetscLogEventEnd(DMPLEX_Migrate, dm, 0, 0, 0);CHKERRQ(ierr); 1927 PetscFunctionReturn(0); 1928 } 1929 1930 /*@C 1931 DMPlexDistribute - Distributes the mesh and any associated sections. 1932 1933 Collective on DM 1934 1935 Input Parameter: 1936 + dm - The original DMPlex object 1937 - overlap - The overlap of partitions, 0 is the default 1938 1939 Output Parameter: 1940 + sf - The PetscSF used for point distribution, or NULL if not needed 1941 - dmParallel - The distributed DMPlex object 1942 1943 Note: If the mesh was not distributed, the output dmParallel will be NULL. 1944 1945 The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function 1946 representation on the mesh. 1947 1948 Level: intermediate 1949 1950 .keywords: mesh, elements 1951 .seealso: DMPlexCreate(), DMPlexDistributeByFace(), DMSetAdjacency() 1952 @*/ 1953 PetscErrorCode DMPlexDistribute(DM dm, PetscInt overlap, PetscSF *sf, DM *dmParallel) 1954 { 1955 MPI_Comm comm; 1956 PetscPartitioner partitioner; 1957 IS cellPart; 1958 PetscSection cellPartSection; 1959 DM dmCoord; 1960 DMLabel lblPartition, lblMigration; 1961 PetscSF sfProcess, sfMigration, sfStratified, sfPoint; 1962 PetscBool flg, balance; 1963 PetscMPIInt rank, size, p; 1964 PetscErrorCode ierr; 1965 1966 PetscFunctionBegin; 1967 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1968 PetscValidLogicalCollectiveInt(dm, overlap, 2); 1969 if (sf) PetscValidPointer(sf,3); 1970 PetscValidPointer(dmParallel,4); 1971 1972 if (sf) *sf = NULL; 1973 *dmParallel = NULL; 1974 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 1975 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1976 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1977 if (size == 1) PetscFunctionReturn(0); 1978 1979 ierr = PetscLogEventBegin(DMPLEX_Distribute,dm,0,0,0);CHKERRQ(ierr); 1980 /* Create cell partition */ 1981 ierr = PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0);CHKERRQ(ierr); 1982 ierr = PetscSectionCreate(comm, &cellPartSection);CHKERRQ(ierr); 1983 ierr = DMPlexGetPartitioner(dm, &partitioner);CHKERRQ(ierr); 1984 ierr = PetscPartitionerPartition(partitioner, dm, cellPartSection, &cellPart);CHKERRQ(ierr); 1985 { 1986 /* Convert partition to DMLabel */ 1987 PetscInt proc, pStart, pEnd, npoints, poffset; 1988 const PetscInt *points; 1989 ierr = DMLabelCreate(PETSC_COMM_SELF, "Point Partition", &lblPartition);CHKERRQ(ierr); 1990 ierr = ISGetIndices(cellPart, &points);CHKERRQ(ierr); 1991 ierr = PetscSectionGetChart(cellPartSection, &pStart, &pEnd);CHKERRQ(ierr); 1992 for (proc = pStart; proc < pEnd; proc++) { 1993 ierr = PetscSectionGetDof(cellPartSection, proc, &npoints);CHKERRQ(ierr); 1994 ierr = PetscSectionGetOffset(cellPartSection, proc, &poffset);CHKERRQ(ierr); 1995 for (p = poffset; p < poffset+npoints; p++) { 1996 ierr = DMLabelSetValue(lblPartition, points[p], proc);CHKERRQ(ierr); 1997 } 1998 } 1999 ierr = ISRestoreIndices(cellPart, &points);CHKERRQ(ierr); 2000 } 2001 ierr = DMPlexPartitionLabelClosure(dm, lblPartition);CHKERRQ(ierr); 2002 { 2003 /* Build a global process SF */ 2004 PetscSFNode *remoteProc; 2005 ierr = PetscMalloc1(size, &remoteProc);CHKERRQ(ierr); 2006 for (p = 0; p < size; ++p) { 2007 remoteProc[p].rank = p; 2008 remoteProc[p].index = rank; 2009 } 2010 ierr = PetscSFCreate(comm, &sfProcess);CHKERRQ(ierr); 2011 ierr = PetscObjectSetName((PetscObject) sfProcess, "Process SF");CHKERRQ(ierr); 2012 ierr = PetscSFSetGraph(sfProcess, size, size, NULL, PETSC_OWN_POINTER, remoteProc, PETSC_OWN_POINTER);CHKERRQ(ierr); 2013 } 2014 ierr = DMLabelCreate(PETSC_COMM_SELF, "Point migration", &lblMigration);CHKERRQ(ierr); 2015 ierr = DMPlexPartitionLabelInvert(dm, lblPartition, sfProcess, lblMigration);CHKERRQ(ierr); 2016 ierr = DMPlexPartitionLabelCreateSF(dm, lblMigration, &sfMigration);CHKERRQ(ierr); 2017 /* Stratify the SF in case we are migrating an already parallel plex */ 2018 ierr = DMPlexStratifyMigrationSF(dm, sfMigration, &sfStratified);CHKERRQ(ierr); 2019 ierr = PetscSFDestroy(&sfMigration);CHKERRQ(ierr); 2020 sfMigration = sfStratified; 2021 ierr = PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0);CHKERRQ(ierr); 2022 ierr = PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);CHKERRQ(ierr); 2023 if (flg) { 2024 ierr = DMLabelView(lblPartition, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 2025 ierr = PetscSFView(sfMigration, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 2026 } 2027 2028 /* Create non-overlapping parallel DM and migrate internal data */ 2029 ierr = DMPlexCreate(comm, dmParallel);CHKERRQ(ierr); 2030 ierr = PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");CHKERRQ(ierr); 2031 ierr = DMPlexMigrate(dm, sfMigration, *dmParallel);CHKERRQ(ierr); 2032 2033 /* Build the point SF without overlap */ 2034 ierr = DMPlexGetPartitionBalance(dm, &balance);CHKERRQ(ierr); 2035 ierr = DMPlexSetPartitionBalance(*dmParallel, balance);CHKERRQ(ierr); 2036 ierr = DMPlexCreatePointSF(*dmParallel, sfMigration, PETSC_TRUE, &sfPoint);CHKERRQ(ierr); 2037 ierr = DMSetPointSF(*dmParallel, sfPoint);CHKERRQ(ierr); 2038 ierr = DMGetCoordinateDM(*dmParallel, &dmCoord);CHKERRQ(ierr); 2039 if (dmCoord) {ierr = DMSetPointSF(dmCoord, sfPoint);CHKERRQ(ierr);} 2040 if (flg) {ierr = PetscSFView(sfPoint, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 2041 2042 if (overlap > 0) { 2043 DM dmOverlap; 2044 PetscInt nroots, nleaves; 2045 PetscSFNode *newRemote; 2046 const PetscSFNode *oldRemote; 2047 PetscSF sfOverlap, sfOverlapPoint; 2048 /* Add the partition overlap to the distributed DM */ 2049 ierr = DMPlexDistributeOverlap(*dmParallel, overlap, &sfOverlap, &dmOverlap);CHKERRQ(ierr); 2050 ierr = DMDestroy(dmParallel);CHKERRQ(ierr); 2051 *dmParallel = dmOverlap; 2052 if (flg) { 2053 ierr = PetscPrintf(comm, "Overlap Migration SF:\n");CHKERRQ(ierr); 2054 ierr = PetscSFView(sfOverlap, NULL);CHKERRQ(ierr); 2055 } 2056 2057 /* Re-map the migration SF to establish the full migration pattern */ 2058 ierr = PetscSFGetGraph(sfMigration, &nroots, NULL, NULL, &oldRemote);CHKERRQ(ierr); 2059 ierr = PetscSFGetGraph(sfOverlap, NULL, &nleaves, NULL, NULL);CHKERRQ(ierr); 2060 ierr = PetscMalloc1(nleaves, &newRemote);CHKERRQ(ierr); 2061 ierr = PetscSFBcastBegin(sfOverlap, MPIU_2INT, oldRemote, newRemote);CHKERRQ(ierr); 2062 ierr = PetscSFBcastEnd(sfOverlap, MPIU_2INT, oldRemote, newRemote);CHKERRQ(ierr); 2063 ierr = PetscSFCreate(comm, &sfOverlapPoint);CHKERRQ(ierr); 2064 ierr = PetscSFSetGraph(sfOverlapPoint, nroots, nleaves, NULL, PETSC_OWN_POINTER, newRemote, PETSC_OWN_POINTER);CHKERRQ(ierr); 2065 ierr = PetscSFDestroy(&sfOverlap);CHKERRQ(ierr); 2066 ierr = PetscSFDestroy(&sfMigration);CHKERRQ(ierr); 2067 sfMigration = sfOverlapPoint; 2068 } 2069 /* Cleanup Partition */ 2070 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr); 2071 ierr = DMLabelDestroy(&lblPartition);CHKERRQ(ierr); 2072 ierr = DMLabelDestroy(&lblMigration);CHKERRQ(ierr); 2073 ierr = PetscSectionDestroy(&cellPartSection);CHKERRQ(ierr); 2074 ierr = ISDestroy(&cellPart);CHKERRQ(ierr); 2075 /* Copy BC */ 2076 ierr = DMCopyBoundary(dm, *dmParallel);CHKERRQ(ierr); 2077 /* Create sfNatural */ 2078 if (dm->useNatural) { 2079 PetscSection section; 2080 2081 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 2082 ierr = DMPlexCreateGlobalToNaturalSF(*dmParallel, section, sfMigration, &(*dmParallel)->sfNatural);CHKERRQ(ierr); 2083 ierr = DMSetUseNatural(*dmParallel, PETSC_TRUE);CHKERRQ(ierr); 2084 } 2085 /* Cleanup */ 2086 if (sf) {*sf = sfMigration;} 2087 else {ierr = PetscSFDestroy(&sfMigration);CHKERRQ(ierr);} 2088 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 2089 ierr = PetscLogEventEnd(DMPLEX_Distribute,dm,0,0,0);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 /*@C 2094 DMPlexDistributeOverlap - Add partition overlap to a distributed non-overlapping DM. 2095 2096 Collective on DM 2097 2098 Input Parameter: 2099 + dm - The non-overlapping distrbuted DMPlex object 2100 - overlap - The overlap of partitions 2101 2102 Output Parameter: 2103 + sf - The PetscSF used for point distribution 2104 - dmOverlap - The overlapping distributed DMPlex object, or NULL 2105 2106 Note: If the mesh was not distributed, the return value is NULL. 2107 2108 The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function 2109 representation on the mesh. 2110 2111 Level: intermediate 2112 2113 .keywords: mesh, elements 2114 .seealso: DMPlexCreate(), DMPlexDistributeByFace(), DMSetAdjacency() 2115 @*/ 2116 PetscErrorCode DMPlexDistributeOverlap(DM dm, PetscInt overlap, PetscSF *sf, DM *dmOverlap) 2117 { 2118 MPI_Comm comm; 2119 PetscMPIInt size, rank; 2120 PetscSection rootSection, leafSection; 2121 IS rootrank, leafrank; 2122 DM dmCoord; 2123 DMLabel lblOverlap; 2124 PetscSF sfOverlap, sfStratified, sfPoint; 2125 PetscErrorCode ierr; 2126 2127 PetscFunctionBegin; 2128 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2129 if (sf) PetscValidPointer(sf, 3); 2130 PetscValidPointer(dmOverlap, 4); 2131 2132 if (sf) *sf = NULL; 2133 *dmOverlap = NULL; 2134 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2135 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 2136 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2137 if (size == 1) PetscFunctionReturn(0); 2138 2139 ierr = PetscLogEventBegin(DMPLEX_DistributeOverlap, dm, 0, 0, 0);CHKERRQ(ierr); 2140 /* Compute point overlap with neighbouring processes on the distributed DM */ 2141 ierr = PetscLogEventBegin(DMPLEX_Partition,dm,0,0,0);CHKERRQ(ierr); 2142 ierr = PetscSectionCreate(comm, &rootSection);CHKERRQ(ierr); 2143 ierr = PetscSectionCreate(comm, &leafSection);CHKERRQ(ierr); 2144 ierr = DMPlexDistributeOwnership(dm, rootSection, &rootrank, leafSection, &leafrank);CHKERRQ(ierr); 2145 ierr = DMPlexCreateOverlap(dm, overlap, rootSection, rootrank, leafSection, leafrank, &lblOverlap);CHKERRQ(ierr); 2146 /* Convert overlap label to stratified migration SF */ 2147 ierr = DMPlexPartitionLabelCreateSF(dm, lblOverlap, &sfOverlap);CHKERRQ(ierr); 2148 ierr = DMPlexStratifyMigrationSF(dm, sfOverlap, &sfStratified);CHKERRQ(ierr); 2149 ierr = PetscSFDestroy(&sfOverlap);CHKERRQ(ierr); 2150 sfOverlap = sfStratified; 2151 ierr = PetscObjectSetName((PetscObject) sfOverlap, "Overlap SF");CHKERRQ(ierr); 2152 ierr = PetscSFSetFromOptions(sfOverlap);CHKERRQ(ierr); 2153 2154 ierr = PetscSectionDestroy(&rootSection);CHKERRQ(ierr); 2155 ierr = PetscSectionDestroy(&leafSection);CHKERRQ(ierr); 2156 ierr = ISDestroy(&rootrank);CHKERRQ(ierr); 2157 ierr = ISDestroy(&leafrank);CHKERRQ(ierr); 2158 ierr = PetscLogEventEnd(DMPLEX_Partition,dm,0,0,0);CHKERRQ(ierr); 2159 2160 /* Build the overlapping DM */ 2161 ierr = DMPlexCreate(comm, dmOverlap);CHKERRQ(ierr); 2162 ierr = PetscObjectSetName((PetscObject) *dmOverlap, "Parallel Mesh");CHKERRQ(ierr); 2163 ierr = DMPlexMigrate(dm, sfOverlap, *dmOverlap);CHKERRQ(ierr); 2164 /* Build the new point SF */ 2165 ierr = DMPlexCreatePointSF(*dmOverlap, sfOverlap, PETSC_FALSE, &sfPoint);CHKERRQ(ierr); 2166 ierr = DMSetPointSF(*dmOverlap, sfPoint);CHKERRQ(ierr); 2167 ierr = DMGetCoordinateDM(*dmOverlap, &dmCoord);CHKERRQ(ierr); 2168 if (dmCoord) {ierr = DMSetPointSF(dmCoord, sfPoint);CHKERRQ(ierr);} 2169 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 2170 /* Cleanup overlap partition */ 2171 ierr = DMLabelDestroy(&lblOverlap);CHKERRQ(ierr); 2172 if (sf) *sf = sfOverlap; 2173 else {ierr = PetscSFDestroy(&sfOverlap);CHKERRQ(ierr);} 2174 ierr = PetscLogEventEnd(DMPLEX_DistributeOverlap, dm, 0, 0, 0);CHKERRQ(ierr); 2175 PetscFunctionReturn(0); 2176 } 2177 2178 /*@C 2179 DMPlexGetGatherDM - Get a copy of the DMPlex that gathers all points on the 2180 root process of the original's communicator. 2181 2182 Collective on DM 2183 2184 Input Parameters: 2185 . dm - the original DMPlex object 2186 2187 Output Parameters: 2188 + sf - the PetscSF used for point distribution (optional) 2189 - gatherMesh - the gathered DM object, or NULL 2190 2191 Level: intermediate 2192 2193 .keywords: mesh 2194 .seealso: DMPlexDistribute(), DMPlexGetRedundantDM() 2195 @*/ 2196 PetscErrorCode DMPlexGetGatherDM(DM dm, PetscSF *sf, DM *gatherMesh) 2197 { 2198 MPI_Comm comm; 2199 PetscMPIInt size; 2200 PetscPartitioner oldPart, gatherPart; 2201 PetscErrorCode ierr; 2202 2203 PetscFunctionBegin; 2204 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2205 PetscValidPointer(gatherMesh,2); 2206 *gatherMesh = NULL; 2207 if (sf) *sf = NULL; 2208 comm = PetscObjectComm((PetscObject)dm); 2209 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2210 if (size == 1) PetscFunctionReturn(0); 2211 ierr = DMPlexGetPartitioner(dm,&oldPart);CHKERRQ(ierr); 2212 ierr = PetscObjectReference((PetscObject)oldPart);CHKERRQ(ierr); 2213 ierr = PetscPartitionerCreate(comm,&gatherPart);CHKERRQ(ierr); 2214 ierr = PetscPartitionerSetType(gatherPart,PETSCPARTITIONERGATHER);CHKERRQ(ierr); 2215 ierr = DMPlexSetPartitioner(dm,gatherPart);CHKERRQ(ierr); 2216 ierr = DMPlexDistribute(dm,0,sf,gatherMesh);CHKERRQ(ierr); 2217 2218 ierr = DMPlexSetPartitioner(dm,oldPart);CHKERRQ(ierr); 2219 ierr = PetscPartitionerDestroy(&gatherPart);CHKERRQ(ierr); 2220 ierr = PetscPartitionerDestroy(&oldPart);CHKERRQ(ierr); 2221 PetscFunctionReturn(0); 2222 } 2223 2224 /*@C 2225 DMPlexGetRedundantDM - Get a copy of the DMPlex that is completely copied on each process. 2226 2227 Collective on DM 2228 2229 Input Parameters: 2230 . dm - the original DMPlex object 2231 2232 Output Parameters: 2233 + sf - the PetscSF used for point distribution (optional) 2234 - redundantMesh - the redundant DM object, or NULL 2235 2236 Level: intermediate 2237 2238 .keywords: mesh 2239 .seealso: DMPlexDistribute(), DMPlexGetGatherDM() 2240 @*/ 2241 PetscErrorCode DMPlexGetRedundantDM(DM dm, PetscSF *sf, DM *redundantMesh) 2242 { 2243 MPI_Comm comm; 2244 PetscMPIInt size, rank; 2245 PetscInt pStart, pEnd, p; 2246 PetscInt numPoints = -1; 2247 PetscSF migrationSF, sfPoint, gatherSF; 2248 DM gatherDM, dmCoord; 2249 PetscSFNode *points; 2250 PetscErrorCode ierr; 2251 2252 PetscFunctionBegin; 2253 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2254 PetscValidPointer(redundantMesh,2); 2255 *redundantMesh = NULL; 2256 comm = PetscObjectComm((PetscObject)dm); 2257 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2258 if (size == 1) { 2259 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 2260 *redundantMesh = dm; 2261 if (sf) *sf = NULL; 2262 PetscFunctionReturn(0); 2263 } 2264 ierr = DMPlexGetGatherDM(dm,&gatherSF,&gatherDM);CHKERRQ(ierr); 2265 if (!gatherDM) PetscFunctionReturn(0); 2266 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2267 ierr = DMPlexGetChart(gatherDM,&pStart,&pEnd);CHKERRQ(ierr); 2268 numPoints = pEnd - pStart; 2269 ierr = MPI_Bcast(&numPoints,1,MPIU_INT,0,comm);CHKERRQ(ierr); 2270 ierr = PetscMalloc1(numPoints,&points);CHKERRQ(ierr); 2271 ierr = PetscSFCreate(comm,&migrationSF);CHKERRQ(ierr); 2272 for (p = 0; p < numPoints; p++) { 2273 points[p].index = p; 2274 points[p].rank = 0; 2275 } 2276 ierr = PetscSFSetGraph(migrationSF,pEnd-pStart,numPoints,NULL,PETSC_OWN_POINTER,points,PETSC_OWN_POINTER);CHKERRQ(ierr); 2277 ierr = DMPlexCreate(comm, redundantMesh);CHKERRQ(ierr); 2278 ierr = PetscObjectSetName((PetscObject) *redundantMesh, "Redundant Mesh");CHKERRQ(ierr); 2279 ierr = DMPlexMigrate(gatherDM, migrationSF, *redundantMesh);CHKERRQ(ierr); 2280 ierr = DMPlexCreatePointSF(*redundantMesh, migrationSF, PETSC_FALSE, &sfPoint);CHKERRQ(ierr); 2281 ierr = DMSetPointSF(*redundantMesh, sfPoint);CHKERRQ(ierr); 2282 ierr = DMGetCoordinateDM(*redundantMesh, &dmCoord);CHKERRQ(ierr); 2283 if (dmCoord) {ierr = DMSetPointSF(dmCoord, sfPoint);CHKERRQ(ierr);} 2284 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 2285 if (sf) { 2286 PetscSF tsf; 2287 2288 ierr = PetscSFCompose(gatherSF,migrationSF,&tsf);CHKERRQ(ierr); 2289 ierr = DMPlexStratifyMigrationSF(dm, tsf, sf);CHKERRQ(ierr); 2290 ierr = PetscSFDestroy(&tsf);CHKERRQ(ierr); 2291 } 2292 ierr = PetscSFDestroy(&migrationSF);CHKERRQ(ierr); 2293 ierr = PetscSFDestroy(&gatherSF);CHKERRQ(ierr); 2294 ierr = DMDestroy(&gatherDM);CHKERRQ(ierr); 2295 PetscFunctionReturn(0); 2296 } 2297