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