Actual source code: plexdistribute.c

petsc-master 2018-05-21
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>

  4: /*@C
  5:   DMPlexSetAdjacencyUser - Define adjacency in the mesh using a user-provided callback

  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

 12:   Level: advanced

 14:   Notes:
 15:      The caller of DMPlexGetAdjacency may need to arrange that a large enough array is available for the adjacency.

 17:      Any setting here overrides other configuration of DMPlex adjacency determination.

 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;

 27:   mesh->useradjacency = user;
 28:   mesh->useradjacencyctx = ctx;
 29:   return(0);
 30: }

 32: /*@C
 33:   DMPlexGetAdjacencyUser - get the user-defined adjacency callback

 35:   Input Parameter:
 36: . dm      - The DM object

 38:   Output Parameters:
 39: - user    - The user callback
 40: - ctx     - context for callback evaluation

 42:   Level: advanced

 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;

 52:   if (user) *user = mesh->useradjacency;
 53:   if (ctx) *ctx = mesh->useradjacencyctx;
 54:   return(0);
 55: }

 57: /*@
 58:   DMPlexSetAdjacencyUseCone - Define adjacency in the mesh using either the cone or the support first

 60:   Input Parameters:
 61: + dm      - The DM object
 62: - useCone - Flag to use the cone first

 64:   Level: intermediate

 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

 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;

 81:   DMGetDS(dm, &prob);
 82:   PetscDSGetNumFields(prob, &Nf);
 83:   if (!Nf) {
 84:     PetscDSGetAdjacency(prob, PETSC_DEFAULT, NULL, &useClosure);
 85:     PetscDSSetAdjacency(prob, PETSC_DEFAULT, useCone, useClosure);
 86:   } else {
 87:     PetscDSGetAdjacency(prob, 0, NULL, &useClosure);
 88:     PetscDSSetAdjacency(prob, 0, useCone, useClosure);
 89:   }
 90:   return(0);
 91: }

 93: /*@
 94:   DMPlexGetAdjacencyUseCone - Query whether adjacency in the mesh uses the cone or the support first

 96:   Input Parameter:
 97: . dm      - The DM object

 99:   Output Parameter:
100: . useCone - Flag to use the cone first

102:   Level: intermediate

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

109: .seealso: DMPlexSetAdjacencyUseCone(), DMPlexSetAdjacencyUseClosure(), DMPlexGetAdjacencyUseClosure(), DMPlexDistribute(), DMPlexPreallocateOperator()
110: @*/
111: PetscErrorCode DMPlexGetAdjacencyUseCone(DM dm, PetscBool *useCone)
112: {
113:   PetscDS        prob;
114:   PetscInt       Nf;

118:   DMGetDS(dm, &prob);
119:   PetscDSGetNumFields(prob, &Nf);
120:   if (!Nf) {
121:     PetscDSGetAdjacency(prob, PETSC_DEFAULT, useCone, NULL);
122:   } else {
123:     PetscDSGetAdjacency(prob, 0, useCone, NULL);
124:   }
125:   return(0);
126: }

128: /*@
129:   DMPlexSetAdjacencyUseClosure - Define adjacency in the mesh using the transitive closure

131:   Input Parameters:
132: + dm      - The DM object
133: - useClosure - Flag to use the closure

135:   Level: intermediate

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

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;

152:   DMGetDS(dm, &prob);
153:   PetscDSGetNumFields(prob, &Nf);
154:   if (!Nf) {
155:     PetscDSGetAdjacency(prob, PETSC_DEFAULT, &useCone, NULL);
156:     PetscDSSetAdjacency(prob, PETSC_DEFAULT, useCone, useClosure);
157:   } else {
158:     PetscDSGetAdjacency(prob, 0, &useCone, NULL);
159:     PetscDSSetAdjacency(prob, 0, useCone, useClosure);
160:   }
161:   return(0);
162: }

164: /*@
165:   DMPlexGetAdjacencyUseClosure - Query whether adjacency in the mesh uses the transitive closure

167:   Input Parameter:
168: . dm      - The DM object

170:   Output Parameter:
171: . useClosure - Flag to use the closure

173:   Level: intermediate

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

180: .seealso: DMPlexSetAdjacencyUseClosure(), DMPlexSetAdjacencyUseCone(), DMPlexGetAdjacencyUseCone(), DMPlexDistribute(), DMPlexPreallocateOperator()
181: @*/
182: PetscErrorCode DMPlexGetAdjacencyUseClosure(DM dm, PetscBool *useClosure)
183: {
184:   PetscDS        prob;
185:   PetscInt       Nf;

189:   DMGetDS(dm, &prob);
190:   PetscDSGetNumFields(prob, &Nf);
191:   if (!Nf) {
192:     PetscDSGetAdjacency(prob, PETSC_DEFAULT, NULL, useClosure);
193:   } else {
194:     PetscDSGetAdjacency(prob, 0, NULL, useClosure);
195:   }
196:   return(0);
197: }

199: /*@
200:   DMPlexSetAdjacencyUseAnchors - Define adjacency in the mesh using the point-to-point constraints.

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.

206:   Level: intermediate

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;

216:   mesh->useAnchors = useAnchors;
217:   return(0);
218: }

220: /*@
221:   DMPlexGetAdjacencyUseAnchors - Query whether adjacency in the mesh uses the point-to-point constraints.

223:   Input Parameter:
224: . dm      - The DM object

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.

229:   Level: intermediate

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;

240:   *useAnchors = mesh->useAnchors;
241:   return(0);
242: }

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;

251:   DMPlexGetConeSize(dm, p, &coneSize);
252:   DMPlexGetCone(dm, p, &cone);
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;

258:     DMPlexGetSupportSize(dm, point, &supportSize);
259:     DMPlexGetSupport(dm, point, &support);
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:   return(0);
269: }

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;

278:   DMPlexGetSupportSize(dm, p, &supportSize);
279:   DMPlexGetSupport(dm, p, &support);
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;

285:     DMPlexGetConeSize(dm, point, &coneSize);
286:     DMPlexGetCone(dm, point, &cone);
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:   return(0);
296: }

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;

305:   DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, &star);
306:   for (s = 0; s < starSize*2; s += 2) {
307:     const PetscInt *closure = NULL;
308:     PetscInt        closureSize, c, q;

310:     DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);
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:     DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);
318:   }
319:   DMPlexRestoreTransitiveClosure(dm, p, useClosure, &starSize, &star);
320:   *adjSize = numAdj;
321:   return(0);
322: }

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;

337:   if (useAnchors) {
338:     DMPlexGetAnchors(dm,&aSec,&aIS);
339:     if (aSec) {
340:       PetscSectionGetMaxDof(aSec,&maxAnchors);
341:       maxAnchors = PetscMax(1,maxAnchors);
342:       PetscSectionGetChart(aSec,&aStart,&aEnd);
343:       ISGetIndices(aIS,&anchors);
344:     }
345:   }
346:   if (!*adj) {
347:     PetscInt depth, coneSeries, supportSeries, maxC, maxS, pStart, pEnd;

349:     DMPlexGetChart(dm, &pStart,&pEnd);
350:     DMPlexGetDepth(dm, &depth);
351:     DMPlexGetMaxSizes(dm, &maxC, &maxS);
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:     PetscMalloc1(asiz,adj);
358:   }
359:   if (*adjSize < 0) *adjSize = asiz;
360:   maxAdjSize = *adjSize;
361:   if (mesh->useradjacency) {
362:     mesh->useradjacency(dm, p, adjSize, *adj, mesh->useradjacencyctx);
363:   } else if (useTransitiveClosure) {
364:     DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, *adj);
365:   } else if (useCone) {
366:     DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, *adj);
367:   } else {
368:     DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, *adj);
369:   }
370:   if (useAnchors && aSec) {
371:     PetscInt origSize = *adjSize;
372:     PetscInt numAdj = origSize;
373:     PetscInt i = 0, j;
374:     PetscInt *orig = *adj;

376:     while (i < origSize) {
377:       PetscInt p = orig[i];
378:       PetscInt aDof = 0;

380:       if (p >= aStart && p < aEnd) {
381:         PetscSectionGetDof(aSec,p,&aDof);
382:       }
383:       if (aDof) {
384:         PetscInt aOff;
385:         PetscInt s, q;

387:         for (j = i + 1; j < numAdj; j++) {
388:           orig[j - 1] = orig[j];
389:         }
390:         origSize--;
391:         numAdj--;
392:         PetscSectionGetOffset(aSec,p,&aOff);
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:     ISRestoreIndices(aIS,&anchors);
406:   }
407:   return(0);
408: }

410: /*@
411:   DMPlexGetAdjacency - Return all points adjacent to the given point

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

419:   Output Parameters:
420: + adjSize - The number of adjacent points
421: - adj - The adjacent points

423:   Level: advanced

425:   Notes:
426:     The user must PetscFree the adj array if it was not passed in.

428: .seealso: DMPlexSetAdjacencyUseCone(), DMPlexSetAdjacencyUseClosure(), DMPlexDistribute(), DMCreateMatrix(), DMPlexPreallocateOperator()
429: @*/
430: PetscErrorCode DMPlexGetAdjacency(DM dm, PetscInt p, PetscInt *adjSize, PetscInt *adj[])
431: {
432:   PetscBool      useCone, useClosure, useAnchors;

439:   DMPlexGetAdjacencyUseCone(dm, &useCone);
440:   DMPlexGetAdjacencyUseClosure(dm, &useClosure);
441:   DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
442:   DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, adjSize, adj);
443:   return(0);
444: }

446: /*@
447:   DMPlexCreateTwoSidedProcessSF - Create an SF which just has process connectivity

449:   Collective on DM

451:   Input Parameters:
452: + dm      - The DM
453: - sfPoint - The PetscSF which encodes point connectivity

455:   Output Parameters:
456: + processRanks - A list of process neighbors, or NULL
457: - sfProcess    - An SF encoding the two-sided process connectivity, or NULL

459:   Level: developer

461: .seealso: PetscSFCreate(), DMPlexCreateProcessSF()
462: @*/
463: PetscErrorCode DMPlexCreateTwoSidedProcessSF(DM dm, PetscSF sfPoint, PetscSection rootRankSection, IS rootRanks, PetscSection leafRankSection, IS leafRanks, IS *processRanks, PetscSF *sfProcess)
464: {
465:   const PetscSFNode *remotePoints;
466:   PetscInt          *localPointsNew;
467:   PetscSFNode       *remotePointsNew;
468:   const PetscInt    *nranks;
469:   PetscInt          *ranksNew;
470:   PetscBT            neighbors;
471:   PetscInt           pStart, pEnd, p, numLeaves, l, numNeighbors, n;
472:   PetscMPIInt        size, proc, rank;
473:   PetscErrorCode     ierr;

480:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
481:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
482:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, NULL, &remotePoints);
483:   PetscBTCreate(size, &neighbors);
484:   PetscBTMemzero(size, neighbors);
485:   /* Compute root-to-leaf process connectivity */
486:   PetscSectionGetChart(rootRankSection, &pStart, &pEnd);
487:   ISGetIndices(rootRanks, &nranks);
488:   for (p = pStart; p < pEnd; ++p) {
489:     PetscInt ndof, noff, n;

491:     PetscSectionGetDof(rootRankSection, p, &ndof);
492:     PetscSectionGetOffset(rootRankSection, p, &noff);
493:     for (n = 0; n < ndof; ++n) {PetscBTSet(neighbors, nranks[noff+n]);}
494:   }
495:   ISRestoreIndices(rootRanks, &nranks);
496:   /* Compute leaf-to-neighbor process connectivity */
497:   PetscSectionGetChart(leafRankSection, &pStart, &pEnd);
498:   ISGetIndices(leafRanks, &nranks);
499:   for (p = pStart; p < pEnd; ++p) {
500:     PetscInt ndof, noff, n;

502:     PetscSectionGetDof(leafRankSection, p, &ndof);
503:     PetscSectionGetOffset(leafRankSection, p, &noff);
504:     for (n = 0; n < ndof; ++n) {PetscBTSet(neighbors, nranks[noff+n]);}
505:   }
506:   ISRestoreIndices(leafRanks, &nranks);
507:   /* Compute leaf-to-root process connectivity */
508:   for (l = 0; l < numLeaves; ++l) {PetscBTSet(neighbors, remotePoints[l].rank);}
509:   /* Calculate edges */
510:   PetscBTClear(neighbors, rank);
511:   for(proc = 0, numNeighbors = 0; proc < size; ++proc) {if (PetscBTLookup(neighbors, proc)) ++numNeighbors;}
512:   PetscMalloc1(numNeighbors, &ranksNew);
513:   PetscMalloc1(numNeighbors, &localPointsNew);
514:   PetscMalloc1(numNeighbors, &remotePointsNew);
515:   for(proc = 0, n = 0; proc < size; ++proc) {
516:     if (PetscBTLookup(neighbors, proc)) {
517:       ranksNew[n]              = proc;
518:       localPointsNew[n]        = proc;
519:       remotePointsNew[n].index = rank;
520:       remotePointsNew[n].rank  = proc;
521:       ++n;
522:     }
523:   }
524:   PetscBTDestroy(&neighbors);
525:   if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numNeighbors, ranksNew, PETSC_OWN_POINTER, processRanks);}
526:   else              {PetscFree(ranksNew);}
527:   if (sfProcess) {
528:     PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
529:     PetscObjectSetName((PetscObject) *sfProcess, "Two-Sided Process SF");
530:     PetscSFSetFromOptions(*sfProcess);
531:     PetscSFSetGraph(*sfProcess, size, numNeighbors, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
532:   }
533:   return(0);
534: }

536: /*@
537:   DMPlexDistributeOwnership - Compute owner information for shared points. This basically gets two-sided for an SF.

539:   Collective on DM

541:   Input Parameter:
542: . dm - The DM

544:   Output Parameters:
545: + rootSection - The number of leaves for a given root point
546: . rootrank    - The rank of each edge into the root point
547: . leafSection - The number of processes sharing a given leaf point
548: - leafrank    - The rank of each process sharing a leaf point

550:   Level: developer

552: .seealso: DMPlexCreateOverlap()
553: @*/
554: PetscErrorCode DMPlexDistributeOwnership(DM dm, PetscSection rootSection, IS *rootrank, PetscSection leafSection, IS *leafrank)
555: {
556:   MPI_Comm        comm;
557:   PetscSF         sfPoint;
558:   const PetscInt *rootdegree;
559:   PetscInt       *myrank, *remoterank;
560:   PetscInt        pStart, pEnd, p, nedges;
561:   PetscMPIInt     rank;
562:   PetscErrorCode  ierr;

565:   PetscObjectGetComm((PetscObject) dm, &comm);
566:   MPI_Comm_rank(comm, &rank);
567:   DMPlexGetChart(dm, &pStart, &pEnd);
568:   DMGetPointSF(dm, &sfPoint);
569:   /* Compute number of leaves for each root */
570:   PetscObjectSetName((PetscObject) rootSection, "Root Section");
571:   PetscSectionSetChart(rootSection, pStart, pEnd);
572:   PetscSFComputeDegreeBegin(sfPoint, &rootdegree);
573:   PetscSFComputeDegreeEnd(sfPoint, &rootdegree);
574:   for (p = pStart; p < pEnd; ++p) {PetscSectionSetDof(rootSection, p, rootdegree[p-pStart]);}
575:   PetscSectionSetUp(rootSection);
576:   /* Gather rank of each leaf to root */
577:   PetscSectionGetStorageSize(rootSection, &nedges);
578:   PetscMalloc1(pEnd-pStart, &myrank);
579:   PetscMalloc1(nedges,  &remoterank);
580:   for (p = 0; p < pEnd-pStart; ++p) myrank[p] = rank;
581:   PetscSFGatherBegin(sfPoint, MPIU_INT, myrank, remoterank);
582:   PetscSFGatherEnd(sfPoint, MPIU_INT, myrank, remoterank);
583:   PetscFree(myrank);
584:   ISCreateGeneral(comm, nedges, remoterank, PETSC_OWN_POINTER, rootrank);
585:   /* Distribute remote ranks to leaves */
586:   PetscObjectSetName((PetscObject) leafSection, "Leaf Section");
587:   DMPlexDistributeFieldIS(dm, sfPoint, rootSection, *rootrank, leafSection, leafrank);
588:   return(0);
589: }

591: /*@C
592:   DMPlexCreateOverlap - Compute owner information for shared points. This basically gets two-sided for an SF.

594:   Collective on DM

596:   Input Parameters:
597: + dm          - The DM
598: . levels      - Number of overlap levels
599: . rootSection - The number of leaves for a given root point
600: . rootrank    - The rank of each edge into the root point
601: . leafSection - The number of processes sharing a given leaf point
602: - leafrank    - The rank of each process sharing a leaf point

604:   Output Parameters:
605: + ovLabel     - DMLabel containing remote overlap contributions as point/rank pairings

607:   Level: developer

609: .seealso: DMPlexDistributeOwnership(), DMPlexDistribute()
610: @*/
611: PetscErrorCode DMPlexCreateOverlap(DM dm, PetscInt levels, PetscSection rootSection, IS rootrank, PetscSection leafSection, IS leafrank, DMLabel *ovLabel)
612: {
613:   MPI_Comm           comm;
614:   DMLabel            ovAdjByRank; /* A DMLabel containing all points adjacent to shared points, separated by rank (value in label) */
615:   PetscSF            sfPoint, sfProc;
616:   const PetscSFNode *remote;
617:   const PetscInt    *local;
618:   const PetscInt    *nrank, *rrank;
619:   PetscInt          *adj = NULL;
620:   PetscInt           pStart, pEnd, p, sStart, sEnd, nleaves, l;
621:   PetscMPIInt        rank, size;
622:   PetscBool          useCone, useClosure, flg;
623:   PetscErrorCode     ierr;

626:   PetscObjectGetComm((PetscObject) dm, &comm);
627:   MPI_Comm_size(comm, &size);
628:   MPI_Comm_rank(comm, &rank);
629:   DMGetPointSF(dm, &sfPoint);
630:   DMPlexGetChart(dm, &pStart, &pEnd);
631:   PetscSectionGetChart(leafSection, &sStart, &sEnd);
632:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
633:   DMLabelCreate("Overlap adjacency", &ovAdjByRank);
634:   /* Handle leaves: shared with the root point */
635:   for (l = 0; l < nleaves; ++l) {
636:     PetscInt adjSize = PETSC_DETERMINE, a;

638:     DMPlexGetAdjacency(dm, local[l], &adjSize, &adj);
639:     for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remote[l].rank);}
640:   }
641:   ISGetIndices(rootrank, &rrank);
642:   ISGetIndices(leafrank, &nrank);
643:   /* Handle roots */
644:   for (p = pStart; p < pEnd; ++p) {
645:     PetscInt adjSize = PETSC_DETERMINE, neighbors = 0, noff, n, a;

647:     if ((p >= sStart) && (p < sEnd)) {
648:       /* Some leaves share a root with other leaves on different processes */
649:       PetscSectionGetDof(leafSection, p, &neighbors);
650:       if (neighbors) {
651:         PetscSectionGetOffset(leafSection, p, &noff);
652:         DMPlexGetAdjacency(dm, p, &adjSize, &adj);
653:         for (n = 0; n < neighbors; ++n) {
654:           const PetscInt remoteRank = nrank[noff+n];

656:           if (remoteRank == rank) continue;
657:           for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remoteRank);}
658:         }
659:       }
660:     }
661:     /* Roots are shared with leaves */
662:     PetscSectionGetDof(rootSection, p, &neighbors);
663:     if (!neighbors) continue;
664:     PetscSectionGetOffset(rootSection, p, &noff);
665:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
666:     for (n = 0; n < neighbors; ++n) {
667:       const PetscInt remoteRank = rrank[noff+n];

669:       if (remoteRank == rank) continue;
670:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(ovAdjByRank, adj[a], remoteRank);}
671:     }
672:   }
673:   PetscFree(adj);
674:   ISRestoreIndices(rootrank, &rrank);
675:   ISRestoreIndices(leafrank, &nrank);
676:   /* Add additional overlap levels */
677:   for (l = 1; l < levels; l++) {
678:     /* Propagate point donations over SF to capture remote connections */
679:     DMPlexPartitionLabelPropagate(dm, ovAdjByRank);
680:     /* Add next level of point donations to the label */
681:     DMPlexPartitionLabelAdjacency(dm, ovAdjByRank);
682:   }
683:   /* We require the closure in the overlap */
684:   DMPlexGetAdjacencyUseCone(dm, &useCone);
685:   DMPlexGetAdjacencyUseClosure(dm, &useClosure);
686:   if (useCone || !useClosure) {
687:     DMPlexPartitionLabelClosure(dm, ovAdjByRank);
688:   }
689:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-overlap_view", &flg);
690:   if (flg) {
691:     DMLabelView(ovAdjByRank, PETSC_VIEWER_STDOUT_WORLD);
692:   }
693:   /* Make global process SF and invert sender to receiver label */
694:   {
695:     /* Build a global process SF */
696:     PetscSFNode *remoteProc;
697:     PetscMalloc1(size, &remoteProc);
698:     for (p = 0; p < size; ++p) {
699:       remoteProc[p].rank  = p;
700:       remoteProc[p].index = rank;
701:     }
702:     PetscSFCreate(comm, &sfProc);
703:     PetscObjectSetName((PetscObject) sfProc, "Process SF");
704:     PetscSFSetGraph(sfProc, size, size, NULL, PETSC_OWN_POINTER, remoteProc, PETSC_OWN_POINTER);
705:   }
706:   DMLabelCreate("Overlap label", ovLabel);
707:   DMPlexPartitionLabelInvert(dm, ovAdjByRank, sfProc, *ovLabel);
708:   /* Add owned points, except for shared local points */
709:   for (p = pStart; p < pEnd; ++p) {DMLabelSetValue(*ovLabel, p, rank);}
710:   for (l = 0; l < nleaves; ++l) {
711:     DMLabelClearValue(*ovLabel, local[l], rank);
712:     DMLabelSetValue(*ovLabel, remote[l].index, remote[l].rank);
713:   }
714:   /* Clean up */
715:   DMLabelDestroy(&ovAdjByRank);
716:   PetscSFDestroy(&sfProc);
717:   return(0);
718: }

720: /*@C
721:   DMPlexCreateOverlapMigrationSF - Create an SF describing the new mesh distribution to make the overlap described by the input SF

723:   Collective on DM

725:   Input Parameters:
726: + dm          - The DM
727: - overlapSF   - The SF mapping ghost points in overlap to owner points on other processes

729:   Output Parameters:
730: + migrationSF - An SF that maps original points in old locations to points in new locations

732:   Level: developer

734: .seealso: DMPlexCreateOverlap(), DMPlexDistribute()
735: @*/
736: PetscErrorCode DMPlexCreateOverlapMigrationSF(DM dm, PetscSF overlapSF, PetscSF *migrationSF)
737: {
738:   MPI_Comm           comm;
739:   PetscMPIInt        rank, size;
740:   PetscInt           d, dim, p, pStart, pEnd, nroots, nleaves, newLeaves, point, numSharedPoints;
741:   PetscInt          *pointDepths, *remoteDepths, *ilocal;
742:   PetscInt          *depthRecv, *depthShift, *depthIdx;
743:   PetscSFNode       *iremote;
744:   PetscSF            pointSF;
745:   const PetscInt    *sharedLocal;
746:   const PetscSFNode *overlapRemote, *sharedRemote;
747:   PetscErrorCode     ierr;

751:   PetscObjectGetComm((PetscObject)dm, &comm);
752:   MPI_Comm_rank(comm, &rank);
753:   MPI_Comm_size(comm, &size);
754:   DMGetDimension(dm, &dim);

756:   /* Before building the migration SF we need to know the new stratum offsets */
757:   PetscSFGetGraph(overlapSF, &nroots, &nleaves, NULL, &overlapRemote);
758:   PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths);
759:   for (d=0; d<dim+1; d++) {
760:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
761:     for (p=pStart; p<pEnd; p++) pointDepths[p] = d;
762:   }
763:   for (p=0; p<nleaves; p++) remoteDepths[p] = -1;
764:   PetscSFBcastBegin(overlapSF, MPIU_INT, pointDepths, remoteDepths);
765:   PetscSFBcastEnd(overlapSF, MPIU_INT, pointDepths, remoteDepths);

767:   /* Count recevied points in each stratum and compute the internal strata shift */
768:   PetscMalloc3(dim+1, &depthRecv, dim+1, &depthShift, dim+1, &depthIdx);
769:   for (d=0; d<dim+1; d++) depthRecv[d]=0;
770:   for (p=0; p<nleaves; p++) depthRecv[remoteDepths[p]]++;
771:   depthShift[dim] = 0;
772:   for (d=0; d<dim; d++) depthShift[d] = depthRecv[dim];
773:   for (d=1; d<dim; d++) depthShift[d] += depthRecv[0];
774:   for (d=dim-2; d>0; d--) depthShift[d] += depthRecv[d+1];
775:   for (d=0; d<dim+1; d++) {
776:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
777:     depthIdx[d] = pStart + depthShift[d];
778:   }

780:   /* Form the overlap SF build an SF that describes the full overlap migration SF */
781:   DMPlexGetChart(dm, &pStart, &pEnd);
782:   newLeaves = pEnd - pStart + nleaves;
783:   PetscMalloc1(newLeaves, &ilocal);
784:   PetscMalloc1(newLeaves, &iremote);
785:   /* First map local points to themselves */
786:   for (d=0; d<dim+1; d++) {
787:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
788:     for (p=pStart; p<pEnd; p++) {
789:       point = p + depthShift[d];
790:       ilocal[point] = point;
791:       iremote[point].index = p;
792:       iremote[point].rank = rank;
793:       depthIdx[d]++;
794:     }
795:   }

797:   /* Add in the remote roots for currently shared points */
798:   DMGetPointSF(dm, &pointSF);
799:   PetscSFGetGraph(pointSF, NULL, &numSharedPoints, &sharedLocal, &sharedRemote);
800:   for (d=0; d<dim+1; d++) {
801:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
802:     for (p=0; p<numSharedPoints; p++) {
803:       if (pStart <= sharedLocal[p] && sharedLocal[p] < pEnd) {
804:         point = sharedLocal[p] + depthShift[d];
805:         iremote[point].index = sharedRemote[p].index;
806:         iremote[point].rank = sharedRemote[p].rank;
807:       }
808:     }
809:   }

811:   /* Now add the incoming overlap points */
812:   for (p=0; p<nleaves; p++) {
813:     point = depthIdx[remoteDepths[p]];
814:     ilocal[point] = point;
815:     iremote[point].index = overlapRemote[p].index;
816:     iremote[point].rank = overlapRemote[p].rank;
817:     depthIdx[remoteDepths[p]]++;
818:   }
819:   PetscFree2(pointDepths,remoteDepths);

821:   PetscSFCreate(comm, migrationSF);
822:   PetscObjectSetName((PetscObject) *migrationSF, "Overlap Migration SF");
823:   PetscSFSetFromOptions(*migrationSF);
824:   DMPlexGetChart(dm, &pStart, &pEnd);
825:   PetscSFSetGraph(*migrationSF, pEnd-pStart, newLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);

827:   PetscFree3(depthRecv, depthShift, depthIdx);
828:   return(0);
829: }

831: /*@
832:   DMPlexStratifyMigrationSF - Rearrange the leaves of a migration sf for stratification.

834:   Input Parameter:
835: + dm          - The DM
836: - sf          - A star forest with non-ordered leaves, usually defining a DM point migration

838:   Output Parameter:
839: . migrationSF - A star forest with added leaf indirection that ensures the resulting DM is stratified

841:   Level: developer

843: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexDistributeOverlap()
844: @*/
845: PetscErrorCode DMPlexStratifyMigrationSF(DM dm, PetscSF sf, PetscSF *migrationSF)
846: {
847:   MPI_Comm           comm;
848:   PetscMPIInt        rank, size;
849:   PetscInt           d, ldepth, depth, p, pStart, pEnd, nroots, nleaves;
850:   PetscInt          *pointDepths, *remoteDepths, *ilocal;
851:   PetscInt          *depthRecv, *depthShift, *depthIdx;
852:   PetscInt           hybEnd[4];
853:   const PetscSFNode *iremote;
854:   PetscErrorCode     ierr;

858:   PetscObjectGetComm((PetscObject) dm, &comm);
859:   MPI_Comm_rank(comm, &rank);
860:   MPI_Comm_size(comm, &size);
861:   DMPlexGetDepth(dm, &ldepth);
862:   MPIU_Allreduce(&ldepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
863:   if ((ldepth >= 0) && (depth != ldepth)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", ldepth, depth);

865:   /* Before building the migration SF we need to know the new stratum offsets */
866:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, &iremote);
867:   PetscMalloc2(nroots, &pointDepths, nleaves, &remoteDepths);
868:   DMPlexGetHybridBounds(dm,&hybEnd[depth],&hybEnd[depth-1],&hybEnd[1],&hybEnd[0]);
869:   for (d = 0; d < depth+1; ++d) {
870:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
871:     for (p = pStart; p < pEnd; ++p) {
872:       if (hybEnd[d] >= 0 && p >= hybEnd[d]) { /* put in a separate value for hybrid points */
873:         pointDepths[p] = 2 * d;
874:       } else {
875:         pointDepths[p] = 2 * d + 1;
876:       }
877:     }
878:   }
879:   for (p = 0; p < nleaves; ++p) remoteDepths[p] = -1;
880:   PetscSFBcastBegin(sf, MPIU_INT, pointDepths, remoteDepths);
881:   PetscSFBcastEnd(sf, MPIU_INT, pointDepths, remoteDepths);
882:   /* Count recevied points in each stratum and compute the internal strata shift */
883:   PetscMalloc3(2*(depth+1), &depthRecv, 2*(depth+1), &depthShift, 2*(depth+1), &depthIdx);
884:   for (d = 0; d < 2*(depth+1); ++d) depthRecv[d] = 0;
885:   for (p = 0; p < nleaves; ++p) depthRecv[remoteDepths[p]]++;
886:   depthShift[2*depth+1] = 0;
887:   for (d = 0; d < 2*depth+1; ++d) depthShift[d] = depthRecv[2 * depth + 1];
888:   for (d = 0; d < 2*depth; ++d) depthShift[d] += depthRecv[2 * depth];
889:   depthShift[0] += depthRecv[1];
890:   for (d = 2; d < 2*depth; ++d) depthShift[d] += depthRecv[1];
891:   for (d = 2; d < 2*depth; ++d) depthShift[d] += depthRecv[0];
892:   for (d = 2 * depth-1; d > 2; --d) {
893:     PetscInt e;

895:     for (e = d -1; e > 1; --e) depthShift[e] += depthRecv[d];
896:   }
897:   for (d = 0; d < 2*(depth+1); ++d) {depthIdx[d] = 0;}
898:   /* Derive a new local permutation based on stratified indices */
899:   PetscMalloc1(nleaves, &ilocal);
900:   for (p = 0; p < nleaves; ++p) {
901:     const PetscInt dep = remoteDepths[p];

903:     ilocal[p] = depthShift[dep] + depthIdx[dep];
904:     depthIdx[dep]++;
905:   }
906:   PetscSFCreate(comm, migrationSF);
907:   PetscObjectSetName((PetscObject) *migrationSF, "Migration SF");
908:   PetscSFSetGraph(*migrationSF, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_COPY_VALUES);
909:   PetscFree2(pointDepths,remoteDepths);
910:   PetscFree3(depthRecv, depthShift, depthIdx);
911:   return(0);
912: }

914: /*@
915:   DMPlexDistributeField - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

917:   Collective on DM

919:   Input Parameters:
920: + dm - The DMPlex object
921: . pointSF - The PetscSF describing the communication pattern
922: . originalSection - The PetscSection for existing data layout
923: - originalVec - The existing data

925:   Output Parameters:
926: + newSection - The PetscSF describing the new data layout
927: - newVec - The new data

929:   Level: developer

931: .seealso: DMPlexDistribute(), DMPlexDistributeFieldIS(), DMPlexDistributeData()
932: @*/
933: PetscErrorCode DMPlexDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
934: {
935:   PetscSF        fieldSF;
936:   PetscInt      *remoteOffsets, fieldSize;
937:   PetscScalar   *originalValues, *newValues;

941:   PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0);
942:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

944:   PetscSectionGetStorageSize(newSection, &fieldSize);
945:   VecSetSizes(newVec, fieldSize, PETSC_DETERMINE);
946:   VecSetType(newVec,dm->vectype);

948:   VecGetArray(originalVec, &originalValues);
949:   VecGetArray(newVec, &newValues);
950:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
951:   PetscFree(remoteOffsets);
952:   PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues);
953:   PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues);
954:   PetscSFDestroy(&fieldSF);
955:   VecRestoreArray(newVec, &newValues);
956:   VecRestoreArray(originalVec, &originalValues);
957:   PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0);
958:   return(0);
959: }

961: /*@
962:   DMPlexDistributeFieldIS - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

964:   Collective on DM

966:   Input Parameters:
967: + dm - The DMPlex object
968: . pointSF - The PetscSF describing the communication pattern
969: . originalSection - The PetscSection for existing data layout
970: - originalIS - The existing data

972:   Output Parameters:
973: + newSection - The PetscSF describing the new data layout
974: - newIS - The new data

976:   Level: developer

978: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexDistributeData()
979: @*/
980: PetscErrorCode DMPlexDistributeFieldIS(DM dm, PetscSF pointSF, PetscSection originalSection, IS originalIS, PetscSection newSection, IS *newIS)
981: {
982:   PetscSF         fieldSF;
983:   PetscInt       *newValues, *remoteOffsets, fieldSize;
984:   const PetscInt *originalValues;
985:   PetscErrorCode  ierr;

988:   PetscLogEventBegin(DMPLEX_DistributeField,dm,0,0,0);
989:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

991:   PetscSectionGetStorageSize(newSection, &fieldSize);
992:   PetscMalloc1(fieldSize, &newValues);

994:   ISGetIndices(originalIS, &originalValues);
995:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
996:   PetscFree(remoteOffsets);
997:   PetscSFBcastBegin(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues);
998:   PetscSFBcastEnd(fieldSF, MPIU_INT, (PetscInt *) originalValues, newValues);
999:   PetscSFDestroy(&fieldSF);
1000:   ISRestoreIndices(originalIS, &originalValues);
1001:   ISCreateGeneral(PetscObjectComm((PetscObject) pointSF), fieldSize, newValues, PETSC_OWN_POINTER, newIS);
1002:   PetscLogEventEnd(DMPLEX_DistributeField,dm,0,0,0);
1003:   return(0);
1004: }

1006: /*@
1007:   DMPlexDistributeData - Distribute field data to match a given PetscSF, usually the SF from mesh distribution

1009:   Collective on DM

1011:   Input Parameters:
1012: + dm - The DMPlex object
1013: . pointSF - The PetscSF describing the communication pattern
1014: . originalSection - The PetscSection for existing data layout
1015: . datatype - The type of data
1016: - originalData - The existing data

1018:   Output Parameters:
1019: + newSection - The PetscSection describing the new data layout
1020: - newData - The new data

1022:   Level: developer

1024: .seealso: DMPlexDistribute(), DMPlexDistributeField()
1025: @*/
1026: PetscErrorCode DMPlexDistributeData(DM dm, PetscSF pointSF, PetscSection originalSection, MPI_Datatype datatype, void *originalData, PetscSection newSection, void **newData)
1027: {
1028:   PetscSF        fieldSF;
1029:   PetscInt      *remoteOffsets, fieldSize;
1030:   PetscMPIInt    dataSize;

1034:   PetscLogEventBegin(DMPLEX_DistributeData,dm,0,0,0);
1035:   PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);

1037:   PetscSectionGetStorageSize(newSection, &fieldSize);
1038:   MPI_Type_size(datatype, &dataSize);
1039:   PetscMalloc(fieldSize * dataSize, newData);

1041:   PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
1042:   PetscFree(remoteOffsets);
1043:   PetscSFBcastBegin(fieldSF, datatype, originalData, *newData);
1044:   PetscSFBcastEnd(fieldSF, datatype, originalData, *newData);
1045:   PetscSFDestroy(&fieldSF);
1046:   PetscLogEventEnd(DMPLEX_DistributeData,dm,0,0,0);
1047:   return(0);
1048: }

1050: static PetscErrorCode DMPlexDistributeCones(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
1051: {
1052:   DM_Plex               *pmesh = (DM_Plex*) (dmParallel)->data;
1053:   MPI_Comm               comm;
1054:   PetscSF                coneSF;
1055:   PetscSection           originalConeSection, newConeSection;
1056:   PetscInt              *remoteOffsets, *cones, *globCones, *newCones, newConesSize;
1057:   PetscBool              flg;
1058:   PetscErrorCode         ierr;

1063:   PetscLogEventBegin(DMPLEX_DistributeCones,dm,0,0,0);

1065:   /* Distribute cone section */
1066:   PetscObjectGetComm((PetscObject)dm, &comm);
1067:   DMPlexGetConeSection(dm, &originalConeSection);
1068:   DMPlexGetConeSection(dmParallel, &newConeSection);
1069:   PetscSFDistributeSection(migrationSF, originalConeSection, &remoteOffsets, newConeSection);
1070:   DMSetUp(dmParallel);
1071:   {
1072:     PetscInt pStart, pEnd, p;

1074:     PetscSectionGetChart(newConeSection, &pStart, &pEnd);
1075:     for (p = pStart; p < pEnd; ++p) {
1076:       PetscInt coneSize;
1077:       PetscSectionGetDof(newConeSection, p, &coneSize);
1078:       pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize);
1079:     }
1080:   }
1081:   /* Communicate and renumber cones */
1082:   PetscSFCreateSectionSF(migrationSF, originalConeSection, remoteOffsets, newConeSection, &coneSF);
1083:   PetscFree(remoteOffsets);
1084:   DMPlexGetCones(dm, &cones);
1085:   if (original) {
1086:     PetscInt numCones;

1088:     PetscSectionGetStorageSize(originalConeSection,&numCones); PetscMalloc1(numCones,&globCones);
1089:     ISLocalToGlobalMappingApplyBlock(original, numCones, cones, globCones);
1090:   }
1091:   else {
1092:     globCones = cones;
1093:   }
1094:   DMPlexGetCones(dmParallel, &newCones);
1095:   PetscSFBcastBegin(coneSF, MPIU_INT, globCones, newCones);
1096:   PetscSFBcastEnd(coneSF, MPIU_INT, globCones, newCones);
1097:   if (original) {
1098:     PetscFree(globCones);
1099:   }
1100:   PetscSectionGetStorageSize(newConeSection, &newConesSize);
1101:   ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, newConesSize, newCones, NULL, newCones);
1102: #if PETSC_USE_DEBUG
1103:   {
1104:     PetscInt  p;
1105:     PetscBool valid = PETSC_TRUE;
1106:     for (p = 0; p < newConesSize; ++p) {
1107:       if (newCones[p] < 0) {valid = PETSC_FALSE; PetscPrintf(PETSC_COMM_SELF, "Point %d not in overlap SF\n", p);}
1108:     }
1109:     if (!valid) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
1110:   }
1111: #endif
1112:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-cones_view", &flg);
1113:   if (flg) {
1114:     PetscPrintf(comm, "Serial Cone Section:\n");
1115:     PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_WORLD);
1116:     PetscPrintf(comm, "Parallel Cone Section:\n");
1117:     PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_WORLD);
1118:     PetscSFView(coneSF, NULL);
1119:   }
1120:   DMPlexGetConeOrientations(dm, &cones);
1121:   DMPlexGetConeOrientations(dmParallel, &newCones);
1122:   PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones);
1123:   PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones);
1124:   PetscSFDestroy(&coneSF);
1125:   PetscLogEventEnd(DMPLEX_DistributeCones,dm,0,0,0);
1126:   /* Create supports and stratify DMPlex */
1127:   {
1128:     PetscInt pStart, pEnd;

1130:     PetscSectionGetChart(pmesh->coneSection, &pStart, &pEnd);
1131:     PetscSectionSetChart(pmesh->supportSection, pStart, pEnd);
1132:   }
1133:   DMPlexSymmetrize(dmParallel);
1134:   DMPlexStratify(dmParallel);
1135:   {
1136:     PetscBool useCone, useClosure, useAnchors;

1138:     DMPlexGetAdjacencyUseCone(dm, &useCone);
1139:     DMPlexGetAdjacencyUseClosure(dm, &useClosure);
1140:     DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
1141:     DMPlexSetAdjacencyUseCone(dmParallel, useCone);
1142:     DMPlexSetAdjacencyUseClosure(dmParallel, useClosure);
1143:     DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors);
1144:   }
1145:   return(0);
1146: }

1148: static PetscErrorCode DMPlexDistributeCoordinates(DM dm, PetscSF migrationSF, DM dmParallel)
1149: {
1150:   MPI_Comm         comm;
1151:   PetscSection     originalCoordSection, newCoordSection;
1152:   Vec              originalCoordinates, newCoordinates;
1153:   PetscInt         bs;
1154:   PetscBool        isper;
1155:   const char      *name;
1156:   const PetscReal *maxCell, *L;
1157:   const DMBoundaryType *bd;
1158:   PetscErrorCode   ierr;


1164:   PetscObjectGetComm((PetscObject)dm, &comm);
1165:   DMGetCoordinateSection(dm, &originalCoordSection);
1166:   DMGetCoordinateSection(dmParallel, &newCoordSection);
1167:   DMGetCoordinatesLocal(dm, &originalCoordinates);
1168:   if (originalCoordinates) {
1169:     VecCreate(PETSC_COMM_SELF, &newCoordinates);
1170:     PetscObjectGetName((PetscObject) originalCoordinates, &name);
1171:     PetscObjectSetName((PetscObject) newCoordinates, name);

1173:     DMPlexDistributeField(dm, migrationSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);
1174:     DMSetCoordinatesLocal(dmParallel, newCoordinates);
1175:     VecGetBlockSize(originalCoordinates, &bs);
1176:     VecSetBlockSize(newCoordinates, bs);
1177:     VecDestroy(&newCoordinates);
1178:   }
1179:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
1180:   DMSetPeriodicity(dmParallel, isper, maxCell, L, bd);
1181:   return(0);
1182: }

1184: /* Here we are assuming that process 0 always has everything */
1185: static PetscErrorCode DMPlexDistributeLabels(DM dm, PetscSF migrationSF, DM dmParallel)
1186: {
1187:   DM_Plex         *mesh = (DM_Plex*) dm->data;
1188:   MPI_Comm         comm;
1189:   DMLabel          depthLabel;
1190:   PetscMPIInt      rank;
1191:   PetscInt         depth, d, numLabels, numLocalLabels, l;
1192:   PetscBool        hasLabels = PETSC_FALSE, lsendDepth, sendDepth;
1193:   PetscObjectState depthState = -1;
1194:   PetscErrorCode   ierr;

1199:   PetscLogEventBegin(DMPLEX_DistributeLabels,dm,0,0,0);
1200:   PetscObjectGetComm((PetscObject)dm, &comm);
1201:   MPI_Comm_rank(comm, &rank);

1203:   /* If the user has changed the depth label, communicate it instead */
1204:   DMPlexGetDepth(dm, &depth);
1205:   DMPlexGetDepthLabel(dm, &depthLabel);
1206:   if (depthLabel) {DMLabelGetState(depthLabel, &depthState);}
1207:   lsendDepth = mesh->depthState != depthState ? PETSC_TRUE : PETSC_FALSE;
1208:   MPIU_Allreduce(&lsendDepth, &sendDepth, 1, MPIU_BOOL, MPI_LOR, comm);
1209:   if (sendDepth) {
1210:     DMRemoveLabel(dmParallel, "depth", &depthLabel);
1211:     DMLabelDestroy(&depthLabel);
1212:   }
1213:   /* Everyone must have either the same number of labels, or none */
1214:   DMGetNumLabels(dm, &numLocalLabels);
1215:   numLabels = numLocalLabels;
1216:   MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm);
1217:   if (numLabels == numLocalLabels) hasLabels = PETSC_TRUE;
1218:   for (l = numLabels-1; l >= 0; --l) {
1219:     DMLabel     label = NULL, labelNew = NULL;
1220:     PetscBool   isDepth, lisOutput = PETSC_TRUE, isOutput;

1222:     if (hasLabels) {DMGetLabelByNum(dm, l, &label);}
1223:     /* Skip "depth" because it is recreated */
1224:     if (hasLabels) {PetscStrcmp(label->name, "depth", &isDepth);}
1225:     MPI_Bcast(&isDepth, 1, MPIU_BOOL, 0, comm);
1226:     if (isDepth && !sendDepth) continue;
1227:     DMLabelDistribute(label, migrationSF, &labelNew);
1228:     if (isDepth) {
1229:       /* Put in any missing strata which can occur if users are managing the depth label themselves */
1230:       PetscInt gdepth;

1232:       MPIU_Allreduce(&depth, &gdepth, 1, MPIU_INT, MPI_MAX, comm);
1233:       if ((depth >= 0) && (gdepth != depth)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent Plex depth %d != %d", depth, gdepth);
1234:       for (d = 0; d <= gdepth; ++d) {
1235:         PetscBool has;

1237:         DMLabelHasStratum(labelNew, d, &has);
1238:         if (!has) {DMLabelAddStratum(labelNew, d);}
1239:       }
1240:     }
1241:     DMAddLabel(dmParallel, labelNew);
1242:     /* Put the output flag in the new label */
1243:     if (hasLabels) {DMGetLabelOutput(dm, label->name, &lisOutput);}
1244:     MPIU_Allreduce(&lisOutput, &isOutput, 1, MPIU_BOOL, MPI_LAND, comm);
1245:     DMSetLabelOutput(dmParallel, labelNew->name, isOutput);
1246:   }
1247:   PetscLogEventEnd(DMPLEX_DistributeLabels,dm,0,0,0);
1248:   return(0);
1249: }

1251: static PetscErrorCode DMPlexDistributeSetupHybrid(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping renumbering, DM dmParallel)
1252: {
1253:   DM_Plex        *mesh  = (DM_Plex*) dm->data;
1254:   DM_Plex        *pmesh = (DM_Plex*) (dmParallel)->data;
1255:   PetscBool      *isHybrid, *isHybridParallel;
1256:   PetscInt        dim, depth, d;
1257:   PetscInt        pStart, pEnd, pStartP, pEndP;
1258:   PetscErrorCode  ierr;


1264:   DMGetDimension(dm, &dim);
1265:   DMPlexGetDepth(dm, &depth);
1266:   DMPlexGetChart(dm,&pStart,&pEnd);
1267:   DMPlexGetChart(dmParallel,&pStartP,&pEndP);
1268:   PetscCalloc2(pEnd-pStart,&isHybrid,pEndP-pStartP,&isHybridParallel);
1269:   for (d = 0; d <= depth; d++) {
1270:     PetscInt hybridMax = (depth == 1 && d == 1) ? mesh->hybridPointMax[dim] : mesh->hybridPointMax[d];

1272:     if (hybridMax >= 0) {
1273:       PetscInt sStart, sEnd, p;

1275:       DMPlexGetDepthStratum(dm,d,&sStart,&sEnd);
1276:       for (p = hybridMax; p < sEnd; p++) isHybrid[p-pStart] = PETSC_TRUE;
1277:     }
1278:   }
1279:   PetscSFBcastBegin(migrationSF,MPIU_BOOL,isHybrid,isHybridParallel);
1280:   PetscSFBcastEnd(migrationSF,MPIU_BOOL,isHybrid,isHybridParallel);
1281:   for (d = 0; d <= dim; d++) pmesh->hybridPointMax[d] = -1;
1282:   for (d = 0; d <= depth; d++) {
1283:     PetscInt sStart, sEnd, p, dd;

1285:     DMPlexGetDepthStratum(dmParallel,d,&sStart,&sEnd);
1286:     dd = (depth == 1 && d == 1) ? dim : d;
1287:     for (p = sStart; p < sEnd; p++) {
1288:       if (isHybridParallel[p-pStartP]) {
1289:         pmesh->hybridPointMax[dd] = p;
1290:         break;
1291:       }
1292:     }
1293:   }
1294:   PetscFree2(isHybrid,isHybridParallel);
1295:   return(0);
1296: }

1298: static PetscErrorCode DMPlexDistributeSetupTree(DM dm, PetscSF migrationSF, ISLocalToGlobalMapping original, ISLocalToGlobalMapping renumbering, DM dmParallel)
1299: {
1300:   DM_Plex        *mesh  = (DM_Plex*) dm->data;
1301:   DM_Plex        *pmesh = (DM_Plex*) (dmParallel)->data;
1302:   MPI_Comm        comm;
1303:   DM              refTree;
1304:   PetscSection    origParentSection, newParentSection;
1305:   PetscInt        *origParents, *origChildIDs;
1306:   PetscBool       flg;
1307:   PetscErrorCode  ierr;

1312:   PetscObjectGetComm((PetscObject)dm, &comm);

1314:   /* Set up tree */
1315:   DMPlexGetReferenceTree(dm,&refTree);
1316:   DMPlexSetReferenceTree(dmParallel,refTree);
1317:   DMPlexGetTree(dm,&origParentSection,&origParents,&origChildIDs,NULL,NULL);
1318:   if (origParentSection) {
1319:     PetscInt        pStart, pEnd;
1320:     PetscInt        *newParents, *newChildIDs, *globParents;
1321:     PetscInt        *remoteOffsetsParents, newParentSize;
1322:     PetscSF         parentSF;

1324:     DMPlexGetChart(dmParallel, &pStart, &pEnd);
1325:     PetscSectionCreate(PetscObjectComm((PetscObject)dmParallel),&newParentSection);
1326:     PetscSectionSetChart(newParentSection,pStart,pEnd);
1327:     PetscSFDistributeSection(migrationSF, origParentSection, &remoteOffsetsParents, newParentSection);
1328:     PetscSFCreateSectionSF(migrationSF, origParentSection, remoteOffsetsParents, newParentSection, &parentSF);
1329:     PetscFree(remoteOffsetsParents);
1330:     PetscSectionGetStorageSize(newParentSection,&newParentSize);
1331:     PetscMalloc2(newParentSize,&newParents,newParentSize,&newChildIDs);
1332:     if (original) {
1333:       PetscInt numParents;

1335:       PetscSectionGetStorageSize(origParentSection,&numParents);
1336:       PetscMalloc1(numParents,&globParents);
1337:       ISLocalToGlobalMappingApplyBlock(original, numParents, origParents, globParents);
1338:     }
1339:     else {
1340:       globParents = origParents;
1341:     }
1342:     PetscSFBcastBegin(parentSF, MPIU_INT, globParents, newParents);
1343:     PetscSFBcastEnd(parentSF, MPIU_INT, globParents, newParents);
1344:     if (original) {
1345:       PetscFree(globParents);
1346:     }
1347:     PetscSFBcastBegin(parentSF, MPIU_INT, origChildIDs, newChildIDs);
1348:     PetscSFBcastEnd(parentSF, MPIU_INT, origChildIDs, newChildIDs);
1349:     ISGlobalToLocalMappingApplyBlock(renumbering,IS_GTOLM_MASK, newParentSize, newParents, NULL, newParents);
1350: #if PETSC_USE_DEBUG
1351:     {
1352:       PetscInt  p;
1353:       PetscBool valid = PETSC_TRUE;
1354:       for (p = 0; p < newParentSize; ++p) {
1355:         if (newParents[p] < 0) {valid = PETSC_FALSE; PetscPrintf(PETSC_COMM_SELF, "Point %d not in overlap SF\n", p);}
1356:       }
1357:       if (!valid) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid global to local map");
1358:     }
1359: #endif
1360:     PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-parents_view", &flg);
1361:     if (flg) {
1362:       PetscPrintf(comm, "Serial Parent Section: \n");
1363:       PetscSectionView(origParentSection, PETSC_VIEWER_STDOUT_WORLD);
1364:       PetscPrintf(comm, "Parallel Parent Section: \n");
1365:       PetscSectionView(newParentSection, PETSC_VIEWER_STDOUT_WORLD);
1366:       PetscSFView(parentSF, NULL);
1367:     }
1368:     DMPlexSetTree(dmParallel,newParentSection,newParents,newChildIDs);
1369:     PetscSectionDestroy(&newParentSection);
1370:     PetscFree2(newParents,newChildIDs);
1371:     PetscSFDestroy(&parentSF);
1372:   }
1373:   pmesh->useAnchors = mesh->useAnchors;
1374:   return(0);
1375: }

1377: PETSC_UNUSED static PetscErrorCode DMPlexDistributeSF(DM dm, PetscSF migrationSF, DM dmParallel)
1378: {
1379:   PetscMPIInt            rank, size;
1380:   MPI_Comm               comm;
1381:   PetscErrorCode         ierr;


1387:   /* Create point SF for parallel mesh */
1388:   PetscLogEventBegin(DMPLEX_DistributeSF,dm,0,0,0);
1389:   PetscObjectGetComm((PetscObject)dm, &comm);
1390:   MPI_Comm_rank(comm, &rank);
1391:   MPI_Comm_size(comm, &size);
1392:   {
1393:     const PetscInt *leaves;
1394:     PetscSFNode    *remotePoints, *rowners, *lowners;
1395:     PetscInt        numRoots, numLeaves, numGhostPoints = 0, p, gp, *ghostPoints;
1396:     PetscInt        pStart, pEnd;

1398:     DMPlexGetChart(dmParallel, &pStart, &pEnd);
1399:     PetscSFGetGraph(migrationSF, &numRoots, &numLeaves, &leaves, NULL);
1400:     PetscMalloc2(numRoots,&rowners,numLeaves,&lowners);
1401:     for (p=0; p<numRoots; p++) {
1402:       rowners[p].rank  = -1;
1403:       rowners[p].index = -1;
1404:     }
1405:     PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners);
1406:     PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners);
1407:     for (p = 0; p < numLeaves; ++p) {
1408:       if (lowners[p].rank < 0 || lowners[p].rank == rank) { /* Either put in a bid or we know we own it */
1409:         lowners[p].rank  = rank;
1410:         lowners[p].index = leaves ? leaves[p] : p;
1411:       } else if (lowners[p].rank >= 0) { /* Point already claimed so flag so that MAXLOC does not listen to us */
1412:         lowners[p].rank  = -2;
1413:         lowners[p].index = -2;
1414:       }
1415:     }
1416:     for (p=0; p<numRoots; p++) { /* Root must not participate in the rediction, flag so that MAXLOC does not use */
1417:       rowners[p].rank  = -3;
1418:       rowners[p].index = -3;
1419:     }
1420:     PetscSFReduceBegin(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC);
1421:     PetscSFReduceEnd(migrationSF, MPIU_2INT, lowners, rowners, MPI_MAXLOC);
1422:     PetscSFBcastBegin(migrationSF, MPIU_2INT, rowners, lowners);
1423:     PetscSFBcastEnd(migrationSF, MPIU_2INT, rowners, lowners);
1424:     for (p = 0; p < numLeaves; ++p) {
1425:       if (lowners[p].rank < 0 || lowners[p].index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cell partition corrupt: point not claimed");
1426:       if (lowners[p].rank != rank) ++numGhostPoints;
1427:     }
1428:     PetscMalloc1(numGhostPoints, &ghostPoints);
1429:     PetscMalloc1(numGhostPoints, &remotePoints);
1430:     for (p = 0, gp = 0; p < numLeaves; ++p) {
1431:       if (lowners[p].rank != rank) {
1432:         ghostPoints[gp]        = leaves ? leaves[p] : p;
1433:         remotePoints[gp].rank  = lowners[p].rank;
1434:         remotePoints[gp].index = lowners[p].index;
1435:         ++gp;
1436:       }
1437:     }
1438:     PetscFree2(rowners,lowners);
1439:     PetscSFSetGraph((dmParallel)->sf, pEnd - pStart, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1440:     PetscSFSetFromOptions((dmParallel)->sf);
1441:   }
1442:   {
1443:     PetscBool useCone, useClosure, useAnchors;

1445:     DMPlexGetAdjacencyUseCone(dm, &useCone);
1446:     DMPlexGetAdjacencyUseClosure(dm, &useClosure);
1447:     DMPlexGetAdjacencyUseAnchors(dm, &useAnchors);
1448:     DMPlexSetAdjacencyUseCone(dmParallel, useCone);
1449:     DMPlexSetAdjacencyUseClosure(dmParallel, useClosure);
1450:     DMPlexSetAdjacencyUseAnchors(dmParallel, useAnchors);
1451:   }
1452:   PetscLogEventEnd(DMPLEX_DistributeSF,dm,0,0,0);
1453:   return(0);
1454: }

1456: /*@C
1457:   DMPlexDerivePointSF - Build a point SF from an SF describing a point migration

1459:   Input Parameter:
1460: + dm          - The source DMPlex object
1461: . migrationSF - The star forest that describes the parallel point remapping
1462: . ownership   - Flag causing a vote to determine point ownership

1464:   Output Parameter:
1465: - pointSF     - The star forest describing the point overlap in the remapped DM

1467:   Level: developer

1469: .seealso: DMPlexDistribute(), DMPlexDistributeOverlap()
1470: @*/
1471: PetscErrorCode DMPlexCreatePointSF(DM dm, PetscSF migrationSF, PetscBool ownership, PetscSF *pointSF)
1472: {
1473:   PetscMPIInt        rank;
1474:   PetscInt           p, nroots, nleaves, idx, npointLeaves;
1475:   PetscInt          *pointLocal;
1476:   const PetscInt    *leaves;
1477:   const PetscSFNode *roots;
1478:   PetscSFNode       *rootNodes, *leafNodes, *pointRemote;
1479:   PetscErrorCode     ierr;

1483:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);

1485:   PetscSFGetGraph(migrationSF, &nroots, &nleaves, &leaves, &roots);
1486:   PetscMalloc2(nroots, &rootNodes, nleaves, &leafNodes);
1487:   if (ownership) {
1488:     /* Point ownership vote: Process with highest rank ownes shared points */
1489:     for (p = 0; p < nleaves; ++p) {
1490:       /* Either put in a bid or we know we own it */
1491:       leafNodes[p].rank  = rank;
1492:       leafNodes[p].index = p;
1493:     }
1494:     for (p = 0; p < nroots; p++) {
1495:       /* Root must not participate in the reduction, flag so that MAXLOC does not use */
1496:       rootNodes[p].rank  = -3;
1497:       rootNodes[p].index = -3;
1498:     }
1499:     PetscSFReduceBegin(migrationSF, MPIU_2INT, leafNodes, rootNodes, MPI_MAXLOC);
1500:     PetscSFReduceEnd(migrationSF, MPIU_2INT, leafNodes, rootNodes, MPI_MAXLOC);
1501:   } else {
1502:     for (p = 0; p < nroots; p++) {
1503:       rootNodes[p].index = -1;
1504:       rootNodes[p].rank = rank;
1505:     };
1506:     for (p = 0; p < nleaves; p++) {
1507:       /* Write new local id into old location */
1508:       if (roots[p].rank == rank) {
1509:         rootNodes[roots[p].index].index = leaves ? leaves[p] : p;
1510:       }
1511:     }
1512:   }
1513:   PetscSFBcastBegin(migrationSF, MPIU_2INT, rootNodes, leafNodes);
1514:   PetscSFBcastEnd(migrationSF, MPIU_2INT, rootNodes, leafNodes);

1516:   for (npointLeaves = 0, p = 0; p < nleaves; p++) {if (leafNodes[p].rank != rank) npointLeaves++;}
1517:   PetscMalloc1(npointLeaves, &pointLocal);
1518:   PetscMalloc1(npointLeaves, &pointRemote);
1519:   for (idx = 0, p = 0; p < nleaves; p++) {
1520:     if (leafNodes[p].rank != rank) {
1521:       pointLocal[idx] = p;
1522:       pointRemote[idx] = leafNodes[p];
1523:       idx++;
1524:     }
1525:   }
1526:   PetscSFCreate(PetscObjectComm((PetscObject) dm), pointSF);
1527:   PetscSFSetFromOptions(*pointSF);
1528:   PetscSFSetGraph(*pointSF, nleaves, npointLeaves, pointLocal, PETSC_OWN_POINTER, pointRemote, PETSC_OWN_POINTER);
1529:   PetscFree2(rootNodes, leafNodes);
1530:   return(0);
1531: }

1533: /*@C
1534:   DMPlexMigrate  - Migrates internal DM data over the supplied star forest

1536:   Input Parameter:
1537: + dm       - The source DMPlex object
1538: . sf       - The star forest communication context describing the migration pattern

1540:   Output Parameter:
1541: - targetDM - The target DMPlex object

1543:   Level: intermediate

1545: .seealso: DMPlexDistribute(), DMPlexDistributeOverlap()
1546: @*/
1547: PetscErrorCode DMPlexMigrate(DM dm, PetscSF sf, DM targetDM)
1548: {
1549:   MPI_Comm               comm;
1550:   PetscInt               dim, nroots;
1551:   PetscSF                sfPoint;
1552:   ISLocalToGlobalMapping ltogMigration;
1553:   ISLocalToGlobalMapping ltogOriginal = NULL;
1554:   PetscBool              flg;
1555:   PetscErrorCode         ierr;

1559:   PetscLogEventBegin(DMPLEX_Migrate, dm, 0, 0, 0);
1560:   PetscObjectGetComm((PetscObject) dm, &comm);
1561:   DMGetDimension(dm, &dim);
1562:   DMSetDimension(targetDM, dim);

1564:   /* Check for a one-to-all distribution pattern */
1565:   DMGetPointSF(dm, &sfPoint);
1566:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
1567:   if (nroots >= 0) {
1568:     IS                     isOriginal;
1569:     PetscInt               n, size, nleaves;
1570:     PetscInt              *numbering_orig, *numbering_new;
1571:     /* Get the original point numbering */
1572:     DMPlexCreatePointNumbering(dm, &isOriginal);
1573:     ISLocalToGlobalMappingCreateIS(isOriginal, &ltogOriginal);
1574:     ISLocalToGlobalMappingGetSize(ltogOriginal, &size);
1575:     ISLocalToGlobalMappingGetBlockIndices(ltogOriginal, (const PetscInt**)&numbering_orig);
1576:     /* Convert to positive global numbers */
1577:     for (n=0; n<size; n++) {if (numbering_orig[n] < 0) numbering_orig[n] = -(numbering_orig[n]+1);}
1578:     /* Derive the new local-to-global mapping from the old one */
1579:     PetscSFGetGraph(sf, NULL, &nleaves, NULL, NULL);
1580:     PetscMalloc1(nleaves, &numbering_new);
1581:     PetscSFBcastBegin(sf, MPIU_INT, (PetscInt *) numbering_orig, numbering_new);
1582:     PetscSFBcastEnd(sf, MPIU_INT, (PetscInt *) numbering_orig, numbering_new);
1583:     ISLocalToGlobalMappingCreate(comm, 1, nleaves, (const PetscInt*) numbering_new, PETSC_OWN_POINTER, &ltogMigration);
1584:     ISLocalToGlobalMappingRestoreIndices(ltogOriginal, (const PetscInt**)&numbering_orig);
1585:     ISDestroy(&isOriginal);
1586:   } else {
1587:     /* One-to-all distribution pattern: We can derive LToG from SF */
1588:     ISLocalToGlobalMappingCreateSF(sf, 0, &ltogMigration);
1589:   }
1590:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);
1591:   if (flg) {
1592:     PetscPrintf(comm, "Point renumbering for DM migration:\n");
1593:     ISLocalToGlobalMappingView(ltogMigration, NULL);
1594:   }
1595:   /* Migrate DM data to target DM */
1596:   DMPlexDistributeCones(dm, sf, ltogOriginal, ltogMigration, targetDM);
1597:   DMPlexDistributeLabels(dm, sf, targetDM);
1598:   DMPlexDistributeCoordinates(dm, sf, targetDM);
1599:   DMPlexDistributeSetupHybrid(dm, sf, ltogMigration, targetDM);
1600:   DMPlexDistributeSetupTree(dm, sf, ltogOriginal, ltogMigration, targetDM);
1601:   ISLocalToGlobalMappingDestroy(&ltogOriginal);
1602:   ISLocalToGlobalMappingDestroy(&ltogMigration);
1603:   PetscLogEventEnd(DMPLEX_Migrate, dm, 0, 0, 0);
1604:   return(0);
1605: }

1607: /*@C
1608:   DMPlexDistribute - Distributes the mesh and any associated sections.

1610:   Not Collective

1612:   Input Parameter:
1613: + dm  - The original DMPlex object
1614: - overlap - The overlap of partitions, 0 is the default

1616:   Output Parameter:
1617: + sf - The PetscSF used for point distribution
1618: - parallelMesh - The distributed DMPlex object, or NULL

1620:   Note: If the mesh was not distributed, the return value is NULL.

1622:   The user can control the definition of adjacency for the mesh using DMPlexSetAdjacencyUseCone() and
1623:   DMPlexSetAdjacencyUseClosure(). They should choose the combination appropriate for the function
1624:   representation on the mesh.

1626:   Level: intermediate

1628: .keywords: mesh, elements
1629: .seealso: DMPlexCreate(), DMPlexDistributeByFace(), DMPlexSetAdjacencyUseCone(), DMPlexSetAdjacencyUseClosure()
1630: @*/
1631: PetscErrorCode DMPlexDistribute(DM dm, PetscInt overlap, PetscSF *sf, DM *dmParallel)
1632: {
1633:   MPI_Comm               comm;
1634:   PetscPartitioner       partitioner;
1635:   IS                     cellPart;
1636:   PetscSection           cellPartSection;
1637:   DM                     dmCoord;
1638:   DMLabel                lblPartition, lblMigration;
1639:   PetscSF                sfProcess, sfMigration, sfStratified, sfPoint;
1640:   PetscBool              flg;
1641:   PetscMPIInt            rank, size, p;
1642:   PetscErrorCode         ierr;


1649:   PetscLogEventBegin(DMPLEX_Distribute,dm,0,0,0);
1650:   PetscObjectGetComm((PetscObject)dm,&comm);
1651:   MPI_Comm_rank(comm, &rank);
1652:   MPI_Comm_size(comm, &size);

1654:   if (sf) *sf = NULL;
1655:   *dmParallel = NULL;
1656:   if (size == 1) {
1657:     PetscLogEventEnd(DMPLEX_Distribute,dm,0,0,0);
1658:     return(0);
1659:   }

1661:   /* Create cell partition */
1662:   PetscLogEventBegin(PETSCPARTITIONER_Partition,dm,0,0,0);
1663:   PetscSectionCreate(comm, &cellPartSection);
1664:   DMPlexGetPartitioner(dm, &partitioner);
1665:   PetscPartitionerPartition(partitioner, dm, cellPartSection, &cellPart);
1666:   {
1667:     /* Convert partition to DMLabel */
1668:     PetscInt proc, pStart, pEnd, npoints, poffset;
1669:     const PetscInt *points;
1670:     DMLabelCreate("Point Partition", &lblPartition);
1671:     ISGetIndices(cellPart, &points);
1672:     PetscSectionGetChart(cellPartSection, &pStart, &pEnd);
1673:     for (proc = pStart; proc < pEnd; proc++) {
1674:       PetscSectionGetDof(cellPartSection, proc, &npoints);
1675:       PetscSectionGetOffset(cellPartSection, proc, &poffset);
1676:       for (p = poffset; p < poffset+npoints; p++) {
1677:         DMLabelSetValue(lblPartition, points[p], proc);
1678:       }
1679:     }
1680:     ISRestoreIndices(cellPart, &points);
1681:   }
1682:   DMPlexPartitionLabelClosure(dm, lblPartition);
1683:   {
1684:     /* Build a global process SF */
1685:     PetscSFNode *remoteProc;
1686:     PetscMalloc1(size, &remoteProc);
1687:     for (p = 0; p < size; ++p) {
1688:       remoteProc[p].rank  = p;
1689:       remoteProc[p].index = rank;
1690:     }
1691:     PetscSFCreate(comm, &sfProcess);
1692:     PetscObjectSetName((PetscObject) sfProcess, "Process SF");
1693:     PetscSFSetGraph(sfProcess, size, size, NULL, PETSC_OWN_POINTER, remoteProc, PETSC_OWN_POINTER);
1694:   }
1695:   DMLabelCreate("Point migration", &lblMigration);
1696:   DMPlexPartitionLabelInvert(dm, lblPartition, sfProcess, lblMigration);
1697:   DMPlexPartitionLabelCreateSF(dm, lblMigration, &sfMigration);
1698:   /* Stratify the SF in case we are migrating an already parallel plex */
1699:   DMPlexStratifyMigrationSF(dm, sfMigration, &sfStratified);
1700:   PetscSFDestroy(&sfMigration);
1701:   sfMigration = sfStratified;
1702:   PetscLogEventEnd(PETSCPARTITIONER_Partition,dm,0,0,0);
1703:   PetscOptionsHasName(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-partition_view", &flg);
1704:   if (flg) {
1705:     DMLabelView(lblPartition, PETSC_VIEWER_STDOUT_WORLD);
1706:     PetscSFView(sfMigration, PETSC_VIEWER_STDOUT_WORLD);
1707:   }

1709:   /* Create non-overlapping parallel DM and migrate internal data */
1710:   DMPlexCreate(comm, dmParallel);
1711:   PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");
1712:   DMPlexMigrate(dm, sfMigration, *dmParallel);

1714:   /* Build the point SF without overlap */
1715:   DMPlexCreatePointSF(*dmParallel, sfMigration, PETSC_TRUE, &sfPoint);
1716:   DMSetPointSF(*dmParallel, sfPoint);
1717:   DMGetCoordinateDM(*dmParallel, &dmCoord);
1718:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
1719:   if (flg) {PetscSFView(sfPoint, PETSC_VIEWER_STDOUT_WORLD);}

1721:   if (overlap > 0) {
1722:     DM                 dmOverlap;
1723:     PetscInt           nroots, nleaves;
1724:     PetscSFNode       *newRemote;
1725:     const PetscSFNode *oldRemote;
1726:     PetscSF            sfOverlap, sfOverlapPoint;
1727:     /* Add the partition overlap to the distributed DM */
1728:     DMPlexDistributeOverlap(*dmParallel, overlap, &sfOverlap, &dmOverlap);
1729:     DMDestroy(dmParallel);
1730:     *dmParallel = dmOverlap;
1731:     if (flg) {
1732:       PetscPrintf(comm, "Overlap Migration SF:\n");
1733:       PetscSFView(sfOverlap, NULL);
1734:     }

1736:     /* Re-map the migration SF to establish the full migration pattern */
1737:     PetscSFGetGraph(sfMigration, &nroots, NULL, NULL, &oldRemote);
1738:     PetscSFGetGraph(sfOverlap, NULL, &nleaves, NULL, NULL);
1739:     PetscMalloc1(nleaves, &newRemote);
1740:     PetscSFBcastBegin(sfOverlap, MPIU_2INT, oldRemote, newRemote);
1741:     PetscSFBcastEnd(sfOverlap, MPIU_2INT, oldRemote, newRemote);
1742:     PetscSFCreate(comm, &sfOverlapPoint);
1743:     PetscSFSetGraph(sfOverlapPoint, nroots, nleaves, NULL, PETSC_OWN_POINTER, newRemote, PETSC_OWN_POINTER);
1744:     PetscSFDestroy(&sfOverlap);
1745:     PetscSFDestroy(&sfMigration);
1746:     sfMigration = sfOverlapPoint;
1747:   }
1748:   /* Cleanup Partition */
1749:   PetscSFDestroy(&sfProcess);
1750:   DMLabelDestroy(&lblPartition);
1751:   DMLabelDestroy(&lblMigration);
1752:   PetscSectionDestroy(&cellPartSection);
1753:   ISDestroy(&cellPart);
1754:   /* Copy BC */
1755:   DMCopyBoundary(dm, *dmParallel);
1756:   /* Create sfNatural */
1757:   if (dm->useNatural) {
1758:     PetscSection section;

1760:     DMGetDefaultSection(dm, &section);
1761:     DMPlexCreateGlobalToNaturalSF(*dmParallel, section, sfMigration, &(*dmParallel)->sfNatural);
1762:     DMSetUseNatural(*dmParallel, PETSC_TRUE);
1763:   }
1764:   /* Cleanup */
1765:   if (sf) {*sf = sfMigration;}
1766:   else    {PetscSFDestroy(&sfMigration);}
1767:   PetscSFDestroy(&sfPoint);
1768:   PetscLogEventEnd(DMPLEX_Distribute,dm,0,0,0);
1769:   return(0);
1770: }

1772: /*@C
1773:   DMPlexDistributeOverlap - Add partition overlap to a distributed non-overlapping DM.

1775:   Not Collective

1777:   Input Parameter:
1778: + dm  - The non-overlapping distrbuted DMPlex object
1779: - overlap - The overlap of partitions, 0 is the default

1781:   Output Parameter:
1782: + sf - The PetscSF used for point distribution
1783: - dmOverlap - The overlapping distributed DMPlex object, or NULL

1785:   Note: If the mesh was not distributed, the return value is NULL.

1787:   The user can control the definition of adjacency for the mesh using DMPlexGetAdjacencyUseCone() and
1788:   DMPlexSetAdjacencyUseClosure(). They should choose the combination appropriate for the function
1789:   representation on the mesh.

1791:   Level: intermediate

1793: .keywords: mesh, elements
1794: .seealso: DMPlexCreate(), DMPlexDistributeByFace(), DMPlexSetAdjacencyUseCone(), DMPlexSetAdjacencyUseClosure()
1795: @*/
1796: PetscErrorCode DMPlexDistributeOverlap(DM dm, PetscInt overlap, PetscSF *sf, DM *dmOverlap)
1797: {
1798:   MPI_Comm               comm;
1799:   PetscMPIInt            size, rank;
1800:   PetscSection           rootSection, leafSection;
1801:   IS                     rootrank, leafrank;
1802:   DM                     dmCoord;
1803:   DMLabel                lblOverlap;
1804:   PetscSF                sfOverlap, sfStratified, sfPoint;
1805:   PetscErrorCode         ierr;


1812:   PetscObjectGetComm((PetscObject)dm,&comm);
1813:   MPI_Comm_size(comm, &size);
1814:   MPI_Comm_rank(comm, &rank);
1815:   if (size == 1) {*dmOverlap = NULL; return(0);}
1816:   PetscLogEventBegin(DMPLEX_DistributeOverlap, dm, 0, 0, 0);

1818:   /* Compute point overlap with neighbouring processes on the distributed DM */
1819:   PetscLogEventBegin(PETSCPARTITIONER_Partition,dm,0,0,0);
1820:   PetscSectionCreate(comm, &rootSection);
1821:   PetscSectionCreate(comm, &leafSection);
1822:   DMPlexDistributeOwnership(dm, rootSection, &rootrank, leafSection, &leafrank);
1823:   DMPlexCreateOverlap(dm, overlap, rootSection, rootrank, leafSection, leafrank, &lblOverlap);
1824:   /* Convert overlap label to stratified migration SF */
1825:   DMPlexPartitionLabelCreateSF(dm, lblOverlap, &sfOverlap);
1826:   DMPlexStratifyMigrationSF(dm, sfOverlap, &sfStratified);
1827:   PetscSFDestroy(&sfOverlap);
1828:   sfOverlap = sfStratified;
1829:   PetscObjectSetName((PetscObject) sfOverlap, "Overlap SF");
1830:   PetscSFSetFromOptions(sfOverlap);

1832:   PetscSectionDestroy(&rootSection);
1833:   PetscSectionDestroy(&leafSection);
1834:   ISDestroy(&rootrank);
1835:   ISDestroy(&leafrank);
1836:   PetscLogEventEnd(PETSCPARTITIONER_Partition,dm,0,0,0);

1838:   /* Build the overlapping DM */
1839:   DMPlexCreate(comm, dmOverlap);
1840:   PetscObjectSetName((PetscObject) *dmOverlap, "Parallel Mesh");
1841:   DMPlexMigrate(dm, sfOverlap, *dmOverlap);
1842:   /* Build the new point SF */
1843:   DMPlexCreatePointSF(*dmOverlap, sfOverlap, PETSC_FALSE, &sfPoint);
1844:   DMSetPointSF(*dmOverlap, sfPoint);
1845:   DMGetCoordinateDM(*dmOverlap, &dmCoord);
1846:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
1847:   PetscSFDestroy(&sfPoint);
1848:   /* Cleanup overlap partition */
1849:   DMLabelDestroy(&lblOverlap);
1850:   if (sf) *sf = sfOverlap;
1851:   else    {PetscSFDestroy(&sfOverlap);}
1852:   PetscLogEventEnd(DMPLEX_DistributeOverlap, dm, 0, 0, 0);
1853:   return(0);
1854: }

1856: /*@C
1857:   DMPlexGetGatherDM - Get a copy of the DMPlex that gathers all points on the
1858:   root process of the original's communicator.

1860:   Input Parameters:
1861: . dm - the original DMPlex object

1863:   Output Parameters:
1864: . gatherMesh - the gathered DM object, or NULL

1866:   Level: intermediate

1868: .keywords: mesh
1869: .seealso: DMPlexDistribute(), DMPlexGetRedundantDM()
1870: @*/
1871: PetscErrorCode DMPlexGetGatherDM(DM dm, DM * gatherMesh)
1872: {
1873:   MPI_Comm       comm;
1874:   PetscMPIInt    size;
1875:   PetscPartitioner oldPart, gatherPart;

1880:   comm = PetscObjectComm((PetscObject)dm);
1881:   MPI_Comm_size(comm,&size);
1882:   *gatherMesh = NULL;
1883:   if (size == 1) return(0);
1884:   DMPlexGetPartitioner(dm,&oldPart);
1885:   PetscObjectReference((PetscObject)oldPart);
1886:   PetscPartitionerCreate(comm,&gatherPart);
1887:   PetscPartitionerSetType(gatherPart,PETSCPARTITIONERGATHER);
1888:   DMPlexSetPartitioner(dm,gatherPart);
1889:   DMPlexDistribute(dm,0,NULL,gatherMesh);
1890:   DMPlexSetPartitioner(dm,oldPart);
1891:   PetscPartitionerDestroy(&gatherPart);
1892:   PetscPartitionerDestroy(&oldPart);
1893:   return(0);
1894: }

1896: /*@C
1897:   DMPlexGetRedundantDM - Get a copy of the DMPlex that is completely copied on each process.

1899:   Input Parameters:
1900: . dm - the original DMPlex object

1902:   Output Parameters:
1903: . redundantMesh - the redundant DM object, or NULL

1905:   Level: intermediate

1907: .keywords: mesh
1908: .seealso: DMPlexDistribute(), DMPlexGetGatherDM()
1909: @*/
1910: PetscErrorCode DMPlexGetRedundantDM(DM dm, DM * redundantMesh)
1911: {
1912:   MPI_Comm       comm;
1913:   PetscMPIInt    size, rank;
1914:   PetscInt       pStart, pEnd, p;
1915:   PetscInt       numPoints = -1;
1916:   PetscSF        migrationSF, sfPoint;
1917:   DM             gatherDM, dmCoord;
1918:   PetscSFNode    *points;

1923:   comm = PetscObjectComm((PetscObject)dm);
1924:   MPI_Comm_size(comm,&size);
1925:   *redundantMesh = NULL;
1926:   if (size == 1) {
1927:     PetscObjectReference((PetscObject) dm);
1928:     *redundantMesh = dm;
1929:     return(0);
1930:   }
1931:   DMPlexGetGatherDM(dm,&gatherDM);
1932:   if (!gatherDM) return(0);
1933:   MPI_Comm_rank(comm,&rank);
1934:   DMPlexGetChart(gatherDM,&pStart,&pEnd);
1935:   numPoints = pEnd - pStart;
1936:   MPI_Bcast(&numPoints,1,MPIU_INT,0,comm);
1937:   PetscMalloc1(numPoints,&points);
1938:   PetscSFCreate(comm,&migrationSF);
1939:   for (p = 0; p < numPoints; p++) {
1940:     points[p].index = p;
1941:     points[p].rank  = 0;
1942:   }
1943:   PetscSFSetGraph(migrationSF,pEnd-pStart,numPoints,NULL,PETSC_OWN_POINTER,points,PETSC_OWN_POINTER);
1944:   DMPlexCreate(comm, redundantMesh);
1945:   PetscObjectSetName((PetscObject) *redundantMesh, "Redundant Mesh");
1946:   DMPlexMigrate(gatherDM, migrationSF, *redundantMesh);
1947:   DMPlexCreatePointSF(*redundantMesh, migrationSF, PETSC_FALSE, &sfPoint);
1948:   DMSetPointSF(*redundantMesh, sfPoint);
1949:   DMGetCoordinateDM(*redundantMesh, &dmCoord);
1950:   if (dmCoord) {DMSetPointSF(dmCoord, sfPoint);}
1951:   PetscSFDestroy(&sfPoint);
1952:   PetscSFDestroy(&migrationSF);
1953:   DMDestroy(&gatherDM);
1954:   return(0);
1955: }