Actual source code: plexlabel.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petscsf.h>

  6: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
  7: {

 11:   PetscNew(label);
 12:   PetscStrallocpy(name, &(*label)->name);

 14:   (*label)->refct          = 1;
 15:   (*label)->numStrata      = 0;
 16:   (*label)->stratumValues  = NULL;
 17:   (*label)->arrayValid     = NULL;
 18:   (*label)->stratumSizes   = NULL;
 19:   (*label)->points         = NULL;
 20:   (*label)->ht             = NULL;
 21:   (*label)->pStart         = -1;
 22:   (*label)->pEnd           = -1;
 23:   (*label)->bt             = NULL;
 24:   return(0);
 25: }

 29: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 30: {
 31:   PetscInt       off;

 34:   if (label->arrayValid[v]) return 0;
 35:   if (v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %d in DMLabelMakeValid_Private\n", v);
 37:   PetscHashISize(label->ht[v], label->stratumSizes[v]);

 39:   PetscMalloc1(label->stratumSizes[v], &label->points[v]);
 40:   off = 0;
 41:   PetscHashIGetKeys(label->ht[v], &off, &(label->points[v][0]));
 42:   if (off != label->stratumSizes[v]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of contributed points %d from value %d should be %d", off, label->stratumValues[v], label->stratumSizes[v]);
 43:   PetscHashIClear(label->ht[v]);
 44:   PetscSortInt(label->stratumSizes[v], label->points[v]);
 45:   if (label->bt) {
 46:     PetscInt p;

 48:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 49:       const PetscInt point = label->points[v][p];

 51:       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
 52:       PetscBTSet(label->bt, point - label->pStart);
 53:     }
 54:   }
 55:   label->arrayValid[v] = PETSC_TRUE;
 56:   return(0);
 57: }

 61: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
 62: {
 63:   PetscInt       v;

 67:   for (v = 0; v < label->numStrata; v++){
 68:     DMLabelMakeValid_Private(label, v);
 69:   }
 70:   return(0);
 71: }

 75: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
 76: {
 77:   PETSC_UNUSED PetscHashIIter ret, iter;
 78:   PetscInt                    p;

 82:   if (!label->arrayValid[v]) return(0);
 83:   for (p = 0; p < label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[v][p], ret, iter);
 84:   PetscFree(label->points[v]);
 85:   label->arrayValid[v] = PETSC_FALSE;
 86:   return(0);
 87: }

 91: static PetscErrorCode DMLabelAddStratum_Private(DMLabel label, PetscInt value)
 92: {
 93:   PetscInt    v, *tmpV, *tmpS, **tmpP;
 94:   PetscHashI *tmpH;
 95:   PetscBool  *tmpB;


100:   PetscMalloc1((label->numStrata+1), &tmpV);
101:   PetscMalloc1((label->numStrata+1), &tmpS);
102:   PetscMalloc1((label->numStrata+1), &tmpH);
103:   PetscMalloc1((label->numStrata+1), &tmpP);
104:   PetscMalloc1((label->numStrata+1), &tmpB);
105:   for (v = 0; v < label->numStrata; ++v) {
106:     tmpV[v] = label->stratumValues[v];
107:     tmpS[v] = label->stratumSizes[v];
108:     tmpH[v] = label->ht[v];
109:     tmpP[v] = label->points[v];
110:     tmpB[v] = label->arrayValid[v];
111:   }
112:   tmpV[v] = value;
113:   tmpS[v] = 0;
114:   PetscHashICreate(tmpH[v]);
115:   tmpP[v] = NULL;
116:   tmpB[v] = PETSC_TRUE;
117:   ++label->numStrata;
118:   PetscFree(label->stratumValues);
119:   PetscFree(label->stratumSizes);
120:   PetscFree(label->ht);
121:   PetscFree(label->points);
122:   PetscFree(label->arrayValid);
123:   label->stratumValues = tmpV;
124:   label->stratumSizes  = tmpS;
125:   label->ht            = tmpH;
126:   label->points        = tmpP;
127:   label->arrayValid    = tmpB;

129:   return(0);
130: }

134: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
135: {
138:   *name = label->name;
139:   return(0);
140: }

144: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
145: {
146:   PetscInt       v;
147:   PetscMPIInt    rank;

151:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
152:   if (label) {
153:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
154:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%d, %d)\n", label->pStart, label->pEnd);}
155:     for (v = 0; v < label->numStrata; ++v) {
156:       const PetscInt value = label->stratumValues[v];
157:       PetscInt       p;

159:       for (p = 0; p < label->stratumSizes[v]; ++p) {
160:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D (%D)\n", rank, label->points[v][p], value);
161:       }
162:     }
163:   }
164:   PetscViewerFlush(viewer);
165:   return(0);
166: }

170: /*@C
171:   DMLabelView - View the label

173:   Input Parameters:
174: + label - The DMLabel
175: - viewer - The PetscViewer

177:   Level: intermediate

179: .seealso: DMLabelCreate(), DMLabelDestroy()
180: @*/
181: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
182: {
183:   PetscBool      iascii;

188:   if (label) {DMLabelMakeAllValid_Private(label);}
189:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
190:   if (iascii) {
191:     DMLabelView_Ascii(label, viewer);
192:   } else SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not supported by this mesh object", ((PetscObject) viewer)->type_name);
193:   return(0);
194: }

198: PetscErrorCode DMLabelDestroy(DMLabel *label)
199: {
200:   PetscInt       v;

204:   if (!(*label)) return(0);
205:   if (--(*label)->refct > 0) return(0);
206:   PetscFree((*label)->name);
207:   PetscFree((*label)->stratumValues);
208:   PetscFree((*label)->stratumSizes);
209:   for (v = 0; v < (*label)->numStrata; ++v) {PetscFree((*label)->points[v]);}
210:   PetscFree((*label)->points);
211:   PetscFree((*label)->arrayValid);
212:   if ((*label)->ht) {
213:     for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
214:     PetscFree((*label)->ht);
215:   }
216:   PetscBTDestroy(&(*label)->bt);
217:   PetscFree(*label);
218:   return(0);
219: }

223: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
224: {
225:   PetscInt       v, q;

229:   DMLabelMakeAllValid_Private(label);
230:   PetscNew(labelnew);
231:   PetscStrallocpy(label->name, &(*labelnew)->name);

233:   (*labelnew)->refct      = 1;
234:   (*labelnew)->numStrata  = label->numStrata;
235:   if (label->numStrata) {
236:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
237:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
238:     PetscMalloc1(label->numStrata, &(*labelnew)->ht);
239:     PetscMalloc1(label->numStrata, &(*labelnew)->points);
240:     PetscMalloc1(label->numStrata, &(*labelnew)->arrayValid);
241:     /* Could eliminate unused space here */
242:     for (v = 0; v < label->numStrata; ++v) {
243:       PetscMalloc1(label->stratumSizes[v], &(*labelnew)->points[v]);
244:       PetscHashICreate((*labelnew)->ht[v]);
245:       (*labelnew)->arrayValid[v]     = PETSC_TRUE;
246:       (*labelnew)->stratumValues[v]  = label->stratumValues[v];
247:       (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
248:       for (q = 0; q < label->stratumSizes[v]; ++q) {
249:         (*labelnew)->points[v][q] = label->points[v][q];
250:       }
251:     }
252:   }
253:   (*labelnew)->pStart = -1;
254:   (*labelnew)->pEnd   = -1;
255:   (*labelnew)->bt     = NULL;
256:   return(0);
257: }

261: /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
262: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
263: {
264:   PetscInt       v;

268:   DMLabelMakeAllValid_Private(label);
269:   if (label->bt) {PetscBTDestroy(&label->bt);}
270:   label->pStart = pStart;
271:   label->pEnd   = pEnd;
272:   PetscBTCreate(pEnd - pStart, &label->bt);
273:   PetscBTMemzero(pEnd - pStart, label->bt);
274:   for (v = 0; v < label->numStrata; ++v) {
275:     PetscInt i;

277:     for (i = 0; i < label->stratumSizes[v]; ++i) {
278:       const PetscInt point = label->points[v][i];

280:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, pStart, pEnd);
281:       PetscBTSet(label->bt, point - pStart);
282:     }
283:   }
284:   return(0);
285: }

289: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
290: {

294:   label->pStart = -1;
295:   label->pEnd   = -1;
296:   if (label->bt) {PetscBTDestroy(&label->bt);}
297:   return(0);
298: }

302: /*@
303:   DMLabelHasValue - Determine whether a label assigns the value to any point

305:   Input Parameters:
306: + label - the DMLabel
307: - value - the value

309:   Output Parameter:
310: . contains - Flag indicating whether the label maps this value to any point

312:   Level: developer

314: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
315: @*/
316: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
317: {
318:   PetscInt v;

322:   for (v = 0; v < label->numStrata; ++v) {
323:     if (value == label->stratumValues[v]) break;
324:   }
325:   *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
326:   return(0);
327: }

331: /*@
332:   DMLabelHasPoint - Determine whether a label assigns a value to a point

334:   Input Parameters:
335: + label - the DMLabel
336: - point - the point

338:   Output Parameter:
339: . contains - Flag indicating whether the label maps this point to a value

341:   Note: The user must call DMLabelCreateIndex() before this function.

343:   Level: developer

345: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
346: @*/
347: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
348: {

353:   DMLabelMakeAllValid_Private(label);
354: #if defined(PETSC_USE_DEBUG)
355:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
356:   if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
357: #endif
358:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
359:   return(0);
360: }

364: /*@
365:   DMLabelStratumHasPoint - Return true if the stratum contains a point

367:   Input Parameters:
368: + label - the DMLabel
369: . value - the stratum value
370: - point - the point

372:   Output Parameter:
373: . contains - true if the stratum contains the point

375:   Level: intermediate

377: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
378: @*/
379: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
380: {
381:   PetscInt       v;

386:   *contains = PETSC_FALSE;
387:   for (v = 0; v < label->numStrata; ++v) {
388:     if (label->stratumValues[v] == value) {
389:       if (label->arrayValid[v]) {
390:         PetscInt i;

392:         PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
393:         if (i >= 0) {
394:           *contains = PETSC_TRUE;
395:           break;
396:         }
397:       } else {
398:         PetscBool has;

400:         PetscHashIHasKey(label->ht[v], point, has);
401:         if (has) {
402:           *contains = PETSC_TRUE;
403:           break;
404:         }
405:       }
406:     }
407:   }
408:   return(0);
409: }

413: /*@
414:   DMLabelGetValue - Return the value a label assigns to a point, or -1

416:   Input Parameters:
417: + label - the DMLabel
418: - point - the point

420:   Output Parameter:
421: . value - The point value, or -1

423:   Level: intermediate

425: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
426: @*/
427: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
428: {
429:   PetscInt       v;

434:   *value = -1;
435:   for (v = 0; v < label->numStrata; ++v) {
436:     if (label->arrayValid[v]) {
437:       PetscInt i;

439:       PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
440:       if (i >= 0) {
441:         *value = label->stratumValues[v];
442:         break;
443:       }
444:     } else {
445:       PetscBool has;

447:       PetscHashIHasKey(label->ht[v], point, has);
448:       if (has) {
449:         *value = label->stratumValues[v];
450:         break;
451:       }
452:     }
453:   }
454:   return(0);
455: }

459: /*@
460:   DMLabelSetValue - Set the value a label assigns to a point

462:   Input Parameters:
463: + label - the DMLabel
464: . point - the point
465: - value - The point value

467:   Level: intermediate

469: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue()
470: @*/
471: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
472: {
473:   PETSC_UNUSED PetscHashIIter iter, ret;
474:   PetscInt                    v;
475:   PetscErrorCode              ierr;

478:   /* Find, or add, label value */
479:   for (v = 0; v < label->numStrata; ++v) {
480:     if (label->stratumValues[v] == value) break;
481:   }
482:   /* Create new table */
483:   if (v >= label->numStrata) {
484:     DMLabelAddStratum_Private(label, value);
485:   }
486:   DMLabelMakeInvalid_Private(label, v);
487:   /* Set key */
488:   PetscHashIPut(label->ht[v], point, ret, iter);
489:   return(0);
490: }

494: /*@
495:   DMLabelClearValue - Clear the value a label assigns to a point

497:   Input Parameters:
498: + label - the DMLabel
499: . point - the point
500: - value - The point value

502:   Level: intermediate

504: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
505: @*/
506: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
507: {
508:   PetscInt       v, p;

512:   /* Find label value */
513:   for (v = 0; v < label->numStrata; ++v) {
514:     if (label->stratumValues[v] == value) break;
515:   }
516:   if (v >= label->numStrata) return(0);
517:   if (label->arrayValid[v]) {
518:     /* Check whether point exists */
519:     PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &p);
520:     if (p >= 0) {
521:       PetscMemmove(&label->points[v][p], &label->points[v][p+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
522:       --label->stratumSizes[v];
523:       if (label->bt) {
524:         if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
525:         PetscBTClear(label->bt, point - label->pStart);
526:       }
527:     }
528:   } else {
529:     PetscHashIDelKey(label->ht[v], point);
530:   }
531:   return(0);
532: }

536: /*@
537:   DMLabelInsertIS - Set all points in the IS to a value

539:   Input Parameters:
540: + label - the DMLabel
541: . is    - the point IS
542: - value - The point value

544:   Level: intermediate

546: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
547: @*/
548: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
549: {
550:   const PetscInt *points;
551:   PetscInt        n, p;
552:   PetscErrorCode  ierr;

556:   ISGetLocalSize(is, &n);
557:   ISGetIndices(is, &points);
558:   for (p = 0; p < n; ++p) {DMLabelSetValue(label, points[p], value);}
559:   ISRestoreIndices(is, &points);
560:   return(0);
561: }

565: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
566: {
569:   *numValues = label->numStrata;
570:   return(0);
571: }

575: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
576: {

581:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
582:   return(0);
583: }

587: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
588: {
589:   PetscInt       v;

594:   *size = 0;
595:   for (v = 0; v < label->numStrata; ++v) {
596:     if (label->stratumValues[v] == value) {
597:       DMLabelMakeValid_Private(label, v);
598:       *size = label->stratumSizes[v];
599:       break;
600:     }
601:   }
602:   return(0);
603: }

607: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
608: {
609:   PetscInt       v;

615:   for (v = 0; v < label->numStrata; ++v) {
616:     if (label->stratumValues[v] != value) continue;
617:     DMLabelMakeValid_Private(label, v);
618:     if (start) *start = label->points[v][0];
619:     if (end)   *end   = label->points[v][label->stratumSizes[v]-1]+1;
620:     break;
621:   }
622:   return(0);
623: }

627: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
628: {
629:   PetscInt       v;

634:   *points = NULL;
635:   for (v = 0; v < label->numStrata; ++v) {
636:     if (label->stratumValues[v] == value) {
637:       DMLabelMakeValid_Private(label, v);
638:       if (label->arrayValid[v]) {
639:         ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[v][0], PETSC_COPY_VALUES, points);
640:         PetscObjectSetName((PetscObject) *points, "indices");
641:       } else {
642:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
643:       }
644:       break;
645:     }
646:   }
647:   return(0);
648: }

652: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
653: {
654:   PetscInt       v;

658:   for (v = 0; v < label->numStrata; ++v) {
659:     if (label->stratumValues[v] == value) break;
660:   }
661:   if (v >= label->numStrata) return(0);
662:   if (label->bt) {
663:     PetscInt i;

665:     for (i = 0; i < label->stratumSizes[v]; ++i) {
666:       const PetscInt point = label->points[v][i];

668:       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, label->pStart, label->pEnd);
669:       PetscBTClear(label->bt, point - label->pStart);
670:     }
671:   }
672:   if (label->arrayValid[v]) {
673:     label->stratumSizes[v] = 0;
674:   } else {
675:     PetscHashIClear(label->ht[v]);
676:   }
677:   return(0);
678: }

682: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
683: {
684:   PetscInt       v;

688:   DMLabelMakeAllValid_Private(label);
689:   label->pStart = start;
690:   label->pEnd   = end;
691:   if (label->bt) {PetscBTDestroy(&label->bt);}
692:   /* Could squish offsets, but would only make sense if I reallocate the storage */
693:   for (v = 0; v < label->numStrata; ++v) {
694:     PetscInt off, q;

696:     for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
697:       const PetscInt point = label->points[v][q];

699:       if ((point < start) || (point >= end)) continue;
700:       label->points[v][off++] = point;
701:     }
702:     label->stratumSizes[v] = off;
703:   }
704:   DMLabelCreateIndex(label, start, end);
705:   return(0);
706: }

710: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
711: {
712:   const PetscInt *perm;
713:   PetscInt        numValues, numPoints, v, q;
714:   PetscErrorCode  ierr;

717:   DMLabelMakeAllValid_Private(label);
718:   DMLabelDuplicate(label, labelNew);
719:   DMLabelGetNumValues(*labelNew, &numValues);
720:   ISGetLocalSize(permutation, &numPoints);
721:   ISGetIndices(permutation, &perm);
722:   for (v = 0; v < numValues; ++v) {
723:     const PetscInt size   = (*labelNew)->stratumSizes[v];

725:     for (q = 0; q < size; ++q) {
726:       const PetscInt point = (*labelNew)->points[v][q];

728:       if ((point < 0) || (point >= numPoints)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [0, %d) for the remapping", point, numPoints);
729:       (*labelNew)->points[v][q] = perm[point];
730:     }
731:     PetscSortInt(size, &(*labelNew)->points[v][0]);
732:   }
733:   ISRestoreIndices(permutation, &perm);
734:   if (label->bt) {
735:     PetscBTDestroy(&label->bt);
736:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
737:   }
738:   return(0);
739: }

743: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
744: {
745:   MPI_Comm       comm;
746:   PetscSection   rootSection, leafSection;
747:   PetscSF        labelSF;
748:   PetscInt       p, pStart, pEnd, l, lStart, lEnd, s, nroots, nleaves, size, dof, offset, stratum;
749:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx, *leafStrata, *strataIdx;
750:   char          *name;
751:   PetscInt       nameSize;
752:   size_t         len = 0;
753:   PetscMPIInt    rank, numProcs;

757:   if (label) {DMLabelMakeAllValid_Private(label);}
758:   PetscObjectGetComm((PetscObject)sf, &comm);
759:   MPI_Comm_rank(comm, &rank);
760:   MPI_Comm_size(comm, &numProcs);
761:   /* Bcast name */
762:   if (!rank) {PetscStrlen(label->name, &len);}
763:   nameSize = len;
764:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
765:   PetscMalloc1(nameSize+1, &name);
766:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
767:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
768:   DMLabelCreate(name, labelNew);
769:   PetscFree(name);
770:   /* Bcast numStrata */
771:   if (!rank) (*labelNew)->numStrata = label->numStrata;
772:   MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
773:   /* Bcast stratumValues */
774:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
775:   if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
776:   MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);
777:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->arrayValid);
778:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->arrayValid[s] = PETSC_TRUE;

780:   /* Build a section detailing strata-per-point, distribute and build SF
781:      from that and then send our points. */
782:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
783:   PetscSectionCreate(comm, &rootSection);
784:   PetscSectionSetChart(rootSection, 0, nroots);
785:   if (label) {
786:     for (s = 0; s < label->numStrata; ++s) {
787:       lStart = 0;
788:       lEnd = label->stratumSizes[s];
789:       for (l=lStart; l<lEnd; l++) {
790:         PetscSectionGetDof(rootSection, label->points[s][l], &dof);
791:         PetscSectionSetDof(rootSection, label->points[s][l], dof+1);
792:       }
793:     }
794:   }
795:   PetscSectionSetUp(rootSection);

797:   /* Create a point-wise array of point strata */
798:   PetscSectionGetStorageSize(rootSection, &size);
799:   PetscMalloc1(size, &rootStrata);
800:   PetscCalloc1(nroots, &rootIdx);
801:   if (label) {
802:     for (s = 0; s < label->numStrata; ++s) {
803:       lStart = 0;
804:       lEnd = label->stratumSizes[s];
805:       for (l=lStart; l<lEnd; l++) {
806:         p = label->points[s][l];
807:         PetscSectionGetOffset(rootSection, p, &offset);
808:         rootStrata[offset+rootIdx[p]++] = s;
809:       }
810:     }
811:   }

813:   /* Build SF that maps label points to remote processes */
814:   PetscSectionCreate(comm, &leafSection);
815:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, leafSection);
816:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, leafSection, &labelSF);

818:   /* Send the strata for each point over the derived SF */
819:   PetscSectionGetStorageSize(leafSection, &size);
820:   PetscMalloc1(size, &leafStrata);
821:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, leafStrata);
822:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, leafStrata);

824:   /* Rebuild the point strata on the receiver */
825:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
826:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
827:   for (p=pStart; p<pEnd; p++) {
828:     PetscSectionGetDof(leafSection, p, &dof);
829:     PetscSectionGetOffset(leafSection, p, &offset);
830:     for (s=0; s<dof; s++) {
831:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
832:     }
833:   }
834:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
835:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
836:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
837:     PetscHashICreate((*labelNew)->ht[s]);
838:     PetscMalloc1((*labelNew)->stratumSizes[s], &(*labelNew)->points[s]);
839:   }

841:   /* Insert points into new strata */
842:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
843:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
844:   for (p=pStart; p<pEnd; p++) {
845:     PetscSectionGetDof(leafSection, p, &dof);
846:     PetscSectionGetOffset(leafSection, p, &offset);
847:     for (s=0; s<dof; s++) {
848:       stratum = leafStrata[offset+s];
849:       (*labelNew)->points[stratum][strataIdx[stratum]++] = p;
850:     }
851:   }
852:   PetscFree(rootStrata);
853:   PetscFree(leafStrata);
854:   PetscFree(rootIdx);
855:   PetscFree(strataIdx);
856:   PetscSectionDestroy(&rootSection);
857:   PetscSectionDestroy(&leafSection);
858:   PetscSFDestroy(&labelSF);
859:   return(0);
860: }


865: /*@C
866:   DMPlexCreateLabel - Create a label of the given name if it does not already exist

868:   Not Collective

870:   Input Parameters:
871: + dm   - The DMPlex object
872: - name - The label name

874:   Level: intermediate

876: .keywords: mesh
877: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
878: @*/
879: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
880: {
881:   DM_Plex        *mesh = (DM_Plex*) dm->data;
882:   PlexLabel      next  = mesh->labels;
883:   PetscBool      flg   = PETSC_FALSE;

889:   while (next) {
890:     PetscStrcmp(name, next->label->name, &flg);
891:     if (flg) break;
892:     next = next->next;
893:   }
894:   if (!flg) {
895:     PlexLabel tmpLabel;

897:     PetscCalloc1(1, &tmpLabel);
898:     DMLabelCreate(name, &tmpLabel->label);
899:     tmpLabel->output = PETSC_TRUE;
900:     tmpLabel->next   = mesh->labels;
901:     mesh->labels     = tmpLabel;
902:   }
903:   return(0);
904: }

908: /*@C
909:   DMPlexGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default

911:   Not Collective

913:   Input Parameters:
914: + dm   - The DMPlex object
915: . name - The label name
916: - point - The mesh point

918:   Output Parameter:
919: . value - The label value for this point, or -1 if the point is not in the label

921:   Level: beginner

923: .keywords: mesh
924: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
925: @*/
926: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
927: {
928:   DMLabel        label;

934:   DMPlexGetLabel(dm, name, &label);
935:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
936:   DMLabelGetValue(label, point, value);
937:   return(0);
938: }

942: /*@C
943:   DMPlexSetLabelValue - Add a point to a Sieve Label with given value

945:   Not Collective

947:   Input Parameters:
948: + dm   - The DMPlex object
949: . name - The label name
950: . point - The mesh point
951: - value - The label value for this point

953:   Output Parameter:

955:   Level: beginner

957: .keywords: mesh
958: .seealso: DMLabelSetValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
959: @*/
960: PetscErrorCode DMPlexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
961: {
962:   DMLabel        label;

968:   DMPlexGetLabel(dm, name, &label);
969:   if (!label) {
970:     DMPlexCreateLabel(dm, name);
971:     DMPlexGetLabel(dm, name, &label);
972:   }
973:   DMLabelSetValue(label, point, value);
974:   return(0);
975: }

979: /*@C
980:   DMPlexClearLabelValue - Remove a point from a Sieve Label with given value

982:   Not Collective

984:   Input Parameters:
985: + dm   - The DMPlex object
986: . name - The label name
987: . point - The mesh point
988: - value - The label value for this point

990:   Output Parameter:

992:   Level: beginner

994: .keywords: mesh
995: .seealso: DMLabelClearValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
996: @*/
997: PetscErrorCode DMPlexClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
998: {
999:   DMLabel        label;

1005:   DMPlexGetLabel(dm, name, &label);
1006:   if (!label) return(0);
1007:   DMLabelClearValue(label, point, value);
1008:   return(0);
1009: }

1013: /*@C
1014:   DMPlexGetLabelSize - Get the number of different integer ids in a Label

1016:   Not Collective

1018:   Input Parameters:
1019: + dm   - The DMPlex object
1020: - name - The label name

1022:   Output Parameter:
1023: . size - The number of different integer ids, or 0 if the label does not exist

1025:   Level: beginner

1027: .keywords: mesh
1028: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
1029: @*/
1030: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
1031: {
1032:   DMLabel        label;

1039:   DMPlexGetLabel(dm, name, &label);
1040:   *size = 0;
1041:   if (!label) return(0);
1042:   DMLabelGetNumValues(label, size);
1043:   return(0);
1044: }

1048: /*@C
1049:   DMPlexGetLabelIdIS - Get the integer ids in a label

1051:   Not Collective

1053:   Input Parameters:
1054: + mesh - The DMPlex object
1055: - name - The label name

1057:   Output Parameter:
1058: . ids - The integer ids, or NULL if the label does not exist

1060:   Level: beginner

1062: .keywords: mesh
1063: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
1064: @*/
1065: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
1066: {
1067:   DMLabel        label;

1074:   DMPlexGetLabel(dm, name, &label);
1075:   *ids = NULL;
1076:   if (!label) return(0);
1077:   DMLabelGetValueIS(label, ids);
1078:   return(0);
1079: }

1083: /*@C
1084:   DMPlexGetStratumSize - Get the number of points in a label stratum

1086:   Not Collective

1088:   Input Parameters:
1089: + dm - The DMPlex object
1090: . name - The label name
1091: - value - The stratum value

1093:   Output Parameter:
1094: . size - The stratum size

1096:   Level: beginner

1098: .keywords: mesh
1099: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
1100: @*/
1101: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1102: {
1103:   DMLabel        label;

1110:   DMPlexGetLabel(dm, name, &label);
1111:   *size = 0;
1112:   if (!label) return(0);
1113:   DMLabelGetStratumSize(label, value, size);
1114:   return(0);
1115: }

1119: /*@C
1120:   DMPlexGetStratumIS - Get the points in a label stratum

1122:   Not Collective

1124:   Input Parameters:
1125: + dm - The DMPlex object
1126: . name - The label name
1127: - value - The stratum value

1129:   Output Parameter:
1130: . points - The stratum points, or NULL if the label does not exist or does not have that value

1132:   Level: beginner

1134: .keywords: mesh
1135: .seealso: DMLabelGetStratumIS(), DMPlexGetStratumSize()
1136: @*/
1137: PetscErrorCode DMPlexGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
1138: {
1139:   DMLabel        label;

1146:   DMPlexGetLabel(dm, name, &label);
1147:   *points = NULL;
1148:   if (!label) return(0);
1149:   DMLabelGetStratumIS(label, value, points);
1150:   return(0);
1151: }

1155: /*@C
1156:   DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label

1158:   Not Collective

1160:   Input Parameters:
1161: + dm   - The DMPlex object
1162: . name - The label name
1163: - value - The label value for this point

1165:   Output Parameter:

1167:   Level: beginner

1169: .keywords: mesh
1170: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1171: @*/
1172: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1173: {
1174:   DMLabel        label;

1180:   DMPlexGetLabel(dm, name, &label);
1181:   if (!label) return(0);
1182:   DMLabelClearStratum(label, value);
1183:   return(0);
1184: }

1188: /*@
1189:   DMPlexGetNumLabels - Return the number of labels defined by the mesh

1191:   Not Collective

1193:   Input Parameter:
1194: . dm   - The DMPlex object

1196:   Output Parameter:
1197: . numLabels - the number of Labels

1199:   Level: intermediate

1201: .keywords: mesh
1202: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1203: @*/
1204: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1205: {
1206:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1207:   PlexLabel next = mesh->labels;
1208:   PetscInt  n    = 0;

1213:   while (next) {++n; next = next->next;}
1214:   *numLabels = n;
1215:   return(0);
1216: }

1220: /*@C
1221:   DMPlexGetLabelName - Return the name of nth label

1223:   Not Collective

1225:   Input Parameters:
1226: + dm - The DMPlex object
1227: - n  - the label number

1229:   Output Parameter:
1230: . name - the label name

1232:   Level: intermediate

1234: .keywords: mesh
1235: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1236: @*/
1237: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1238: {
1239:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1240:   PlexLabel next = mesh->labels;
1241:   PetscInt  l    = 0;

1246:   while (next) {
1247:     if (l == n) {
1248:       *name = next->label->name;
1249:       return(0);
1250:     }
1251:     ++l;
1252:     next = next->next;
1253:   }
1254:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1255: }

1259: /*@C
1260:   DMPlexHasLabel - Determine whether the mesh has a label of a given name

1262:   Not Collective

1264:   Input Parameters:
1265: + dm   - The DMPlex object
1266: - name - The label name

1268:   Output Parameter:
1269: . hasLabel - PETSC_TRUE if the label is present

1271:   Level: intermediate

1273: .keywords: mesh
1274: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1275: @*/
1276: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1277: {
1278:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1279:   PlexLabel      next = mesh->labels;

1286:   *hasLabel = PETSC_FALSE;
1287:   while (next) {
1288:     PetscStrcmp(name, next->label->name, hasLabel);
1289:     if (*hasLabel) break;
1290:     next = next->next;
1291:   }
1292:   return(0);
1293: }

1297: /*@C
1298:   DMPlexGetLabel - Return the label of a given name, or NULL

1300:   Not Collective

1302:   Input Parameters:
1303: + dm   - The DMPlex object
1304: - name - The label name

1306:   Output Parameter:
1307: . label - The DMLabel, or NULL if the label is absent

1309:   Level: intermediate

1311: .keywords: mesh
1312: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1313: @*/
1314: PetscErrorCode DMPlexGetLabel(DM dm, const char name[], DMLabel *label)
1315: {
1316:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1317:   PlexLabel      next = mesh->labels;
1318:   PetscBool      hasLabel;

1325:   *label = NULL;
1326:   while (next) {
1327:     PetscStrcmp(name, next->label->name, &hasLabel);
1328:     if (hasLabel) {
1329:       *label = next->label;
1330:       break;
1331:     }
1332:     next = next->next;
1333:   }
1334:   return(0);
1335: }

1339: /*@C
1340:   DMPlexGetLabelByNum - Return the nth label

1342:   Not Collective

1344:   Input Parameters:
1345: + dm - The DMPlex object
1346: - n  - the label number

1348:   Output Parameter:
1349: . label - the label

1351:   Level: intermediate

1353: .keywords: mesh
1354: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1355: @*/
1356: PetscErrorCode DMPlexGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
1357: {
1358:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1359:   PlexLabel next = mesh->labels;
1360:   PetscInt  l    = 0;

1365:   while (next) {
1366:     if (l == n) {
1367:       *label = next->label;
1368:       return(0);
1369:     }
1370:     ++l;
1371:     next = next->next;
1372:   }
1373:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1374: }

1378: /*@C
1379:   DMPlexAddLabel - Add the label to this mesh

1381:   Not Collective

1383:   Input Parameters:
1384: + dm   - The DMPlex object
1385: - label - The DMLabel

1387:   Level: developer

1389: .keywords: mesh
1390: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1391: @*/
1392: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1393: {
1394:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1395:   PlexLabel      tmpLabel;
1396:   PetscBool      hasLabel;

1401:   DMPlexHasLabel(dm, label->name, &hasLabel);
1402:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1403:   PetscCalloc1(1, &tmpLabel);
1404:   tmpLabel->label  = label;
1405:   tmpLabel->output = PETSC_TRUE;
1406:   tmpLabel->next   = mesh->labels;
1407:   mesh->labels     = tmpLabel;
1408:   return(0);
1409: }

1413: /*@C
1414:   DMPlexRemoveLabel - Remove the label from this mesh

1416:   Not Collective

1418:   Input Parameters:
1419: + dm   - The DMPlex object
1420: - name - The label name

1422:   Output Parameter:
1423: . label - The DMLabel, or NULL if the label is absent

1425:   Level: developer

1427: .keywords: mesh
1428: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1429: @*/
1430: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1431: {
1432:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1433:   PlexLabel      next = mesh->labels;
1434:   PlexLabel      last = NULL;
1435:   PetscBool      hasLabel;

1440:   DMPlexHasLabel(dm, name, &hasLabel);
1441:   *label = NULL;
1442:   if (!hasLabel) return(0);
1443:   while (next) {
1444:     PetscStrcmp(name, next->label->name, &hasLabel);
1445:     if (hasLabel) {
1446:       if (last) last->next   = next->next;
1447:       else      mesh->labels = next->next;
1448:       next->next = NULL;
1449:       *label     = next->label;
1450:       PetscFree(next);
1451:       break;
1452:     }
1453:     last = next;
1454:     next = next->next;
1455:   }
1456:   return(0);
1457: }

1461: /*@C
1462:   DMPlexGetLabelOutput - Get the output flag for a given label

1464:   Not Collective

1466:   Input Parameters:
1467: + dm   - The DMPlex object
1468: - name - The label name

1470:   Output Parameter:
1471: . output - The flag for output

1473:   Level: developer

1475: .keywords: mesh
1476: .seealso: DMPlexSetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1477: @*/
1478: PetscErrorCode DMPlexGetLabelOutput(DM dm, const char name[], PetscBool *output)
1479: {
1480:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1481:   PlexLabel      next = mesh->labels;

1488:   while (next) {
1489:     PetscBool flg;

1491:     PetscStrcmp(name, next->label->name, &flg);
1492:     if (flg) {*output = next->output; return(0);}
1493:     next = next->next;
1494:   }
1495:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1496: }

1500: /*@C
1501:   DMPlexSetLabelOutput - Set the output flag for a given label

1503:   Not Collective

1505:   Input Parameters:
1506: + dm     - The DMPlex object
1507: . name   - The label name
1508: - output - The flag for output

1510:   Level: developer

1512: .keywords: mesh
1513: .seealso: DMPlexGetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1514: @*/
1515: PetscErrorCode DMPlexSetLabelOutput(DM dm, const char name[], PetscBool output)
1516: {
1517:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1518:   PlexLabel      next = mesh->labels;

1524:   while (next) {
1525:     PetscBool flg;

1527:     PetscStrcmp(name, next->label->name, &flg);
1528:     if (flg) {next->output = output; return(0);}
1529:     next = next->next;
1530:   }
1531:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1532: }