Actual source code: plexlabel.c

petsc-master 2015-05-04
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: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
537: {
540:   *numValues = label->numStrata;
541:   return(0);
542: }

546: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
547: {

552:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
553:   return(0);
554: }

558: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
559: {
560:   PetscInt       v;

565:   *size = 0;
566:   for (v = 0; v < label->numStrata; ++v) {
567:     if (label->stratumValues[v] == value) {
568:       DMLabelMakeValid_Private(label, v);
569:       *size = label->stratumSizes[v];
570:       break;
571:     }
572:   }
573:   return(0);
574: }

578: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
579: {
580:   PetscInt       v;

586:   for (v = 0; v < label->numStrata; ++v) {
587:     if (label->stratumValues[v] != value) continue;
588:     DMLabelMakeValid_Private(label, v);
589:     if (start) *start = label->points[v][0];
590:     if (end)   *end   = label->points[v][label->stratumSizes[v]-1]+1;
591:     break;
592:   }
593:   return(0);
594: }

598: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
599: {
600:   PetscInt       v;

605:   *points = NULL;
606:   for (v = 0; v < label->numStrata; ++v) {
607:     if (label->stratumValues[v] == value) {
608:       DMLabelMakeValid_Private(label, v);
609:       if (label->arrayValid[v]) {
610:         ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[v][0], PETSC_COPY_VALUES, points);
611:         PetscObjectSetName((PetscObject) *points, "indices");
612:       } else {
613:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
614:       }
615:       break;
616:     }
617:   }
618:   return(0);
619: }

623: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
624: {
625:   PetscInt       v;

629:   for (v = 0; v < label->numStrata; ++v) {
630:     if (label->stratumValues[v] == value) break;
631:   }
632:   if (v >= label->numStrata) return(0);
633:   if (label->bt) {
634:     PetscInt i;

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

639:       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);
640:       PetscBTClear(label->bt, point - label->pStart);
641:     }
642:   }
643:   if (label->arrayValid[v]) {
644:     label->stratumSizes[v] = 0;
645:   } else {
646:     PetscHashIClear(label->ht[v]);
647:   }
648:   return(0);
649: }

653: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
654: {
655:   PetscInt       v;

659:   DMLabelMakeAllValid_Private(label);
660:   label->pStart = start;
661:   label->pEnd   = end;
662:   if (label->bt) {PetscBTDestroy(&label->bt);}
663:   /* Could squish offsets, but would only make sense if I reallocate the storage */
664:   for (v = 0; v < label->numStrata; ++v) {
665:     PetscInt off, q;

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

670:       if ((point < start) || (point >= end)) continue;
671:       label->points[v][off++] = point;
672:     }
673:     label->stratumSizes[v] = off;
674:   }
675:   DMLabelCreateIndex(label, start, end);
676:   return(0);
677: }

681: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
682: {
683:   const PetscInt *perm;
684:   PetscInt        numValues, numPoints, v, q;
685:   PetscErrorCode  ierr;

688:   DMLabelMakeAllValid_Private(label);
689:   DMLabelDuplicate(label, labelNew);
690:   DMLabelGetNumValues(*labelNew, &numValues);
691:   ISGetLocalSize(permutation, &numPoints);
692:   ISGetIndices(permutation, &perm);
693:   for (v = 0; v < numValues; ++v) {
694:     const PetscInt size   = (*labelNew)->stratumSizes[v];

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

699:       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);
700:       (*labelNew)->points[v][q] = perm[point];
701:     }
702:     PetscSortInt(size, &(*labelNew)->points[v][0]);
703:   }
704:   ISRestoreIndices(permutation, &perm);
705:   if (label->bt) {
706:     PetscBTDestroy(&label->bt);
707:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
708:   }
709:   return(0);
710: }

714: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
715: {
716:   MPI_Comm       comm;
717:   PetscSection   rootSection, leafSection;
718:   PetscSF        labelSF;
719:   PetscInt       p, pStart, pEnd, l, lStart, lEnd, s, nroots, nleaves, size, dof, offset, stratum;
720:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx, *leafStrata, *strataIdx;
721:   char          *name;
722:   PetscInt       nameSize;
723:   size_t         len = 0;
724:   PetscMPIInt    rank, numProcs;

728:   if (label) {DMLabelMakeAllValid_Private(label);}
729:   PetscObjectGetComm((PetscObject)sf, &comm);
730:   MPI_Comm_rank(comm, &rank);
731:   MPI_Comm_size(comm, &numProcs);
732:   /* Bcast name */
733:   if (!rank) {PetscStrlen(label->name, &len);}
734:   nameSize = len;
735:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
736:   PetscMalloc1(nameSize+1, &name);
737:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
738:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
739:   DMLabelCreate(name, labelNew);
740:   PetscFree(name);
741:   /* Bcast numStrata */
742:   if (!rank) (*labelNew)->numStrata = label->numStrata;
743:   MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
744:   /* Bcast stratumValues */
745:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
746:   if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
747:   MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);
748:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->arrayValid);
749:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->arrayValid[s] = PETSC_TRUE;

751:   /* Build a section detailing strata-per-point, distribute and build SF
752:      from that and then send our points. */
753:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
754:   PetscSectionCreate(comm, &rootSection);
755:   PetscSectionSetChart(rootSection, 0, nroots);
756:   if (label) {
757:     for (s = 0; s < label->numStrata; ++s) {
758:       lStart = 0;
759:       lEnd = label->stratumSizes[s];
760:       for (l=lStart; l<lEnd; l++) {
761:         PetscSectionGetDof(rootSection, label->points[s][l], &dof);
762:         PetscSectionSetDof(rootSection, label->points[s][l], dof+1);
763:       }
764:     }
765:   }
766:   PetscSectionSetUp(rootSection);

768:   /* Create a point-wise array of point strata */
769:   PetscSectionGetStorageSize(rootSection, &size);
770:   PetscMalloc1(size, &rootStrata);
771:   PetscCalloc1(nroots, &rootIdx);
772:   if (label) {
773:     for (s = 0; s < label->numStrata; ++s) {
774:       lStart = 0;
775:       lEnd = label->stratumSizes[s];
776:       for (l=lStart; l<lEnd; l++) {
777:         p = label->points[s][l];
778:         PetscSectionGetOffset(rootSection, p, &offset);
779:         rootStrata[offset+rootIdx[p]++] = s;
780:       }
781:     }
782:   }

784:   /* Build SF that maps label points to remote processes */
785:   PetscSectionCreate(comm, &leafSection);
786:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, leafSection);
787:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, leafSection, &labelSF);

789:   /* Send the strata for each point over the derived SF */
790:   PetscSectionGetStorageSize(leafSection, &size);
791:   PetscMalloc1(size, &leafStrata);
792:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, leafStrata);
793:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, leafStrata);

795:   /* Rebuild the point strata on the receiver */
796:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
797:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
798:   for (p=pStart; p<pEnd; p++) {
799:     PetscSectionGetDof(leafSection, p, &dof);
800:     PetscSectionGetOffset(leafSection, p, &offset);
801:     for (s=0; s<dof; s++) {
802:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
803:     }
804:   }
805:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
806:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
807:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
808:     PetscHashICreate((*labelNew)->ht[s]);
809:     PetscMalloc1((*labelNew)->stratumSizes[s], &(*labelNew)->points[s]);
810:   }

812:   /* Insert points into new strata */
813:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
814:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
815:   for (p=pStart; p<pEnd; p++) {
816:     PetscSectionGetDof(leafSection, p, &dof);
817:     PetscSectionGetOffset(leafSection, p, &offset);
818:     for (s=0; s<dof; s++) {
819:       stratum = leafStrata[offset+s];
820:       (*labelNew)->points[stratum][strataIdx[stratum]++] = p;
821:     }
822:   }
823:   PetscFree(rootStrata);
824:   PetscFree(leafStrata);
825:   PetscFree(rootIdx);
826:   PetscFree(strataIdx);
827:   PetscSectionDestroy(&rootSection);
828:   PetscSectionDestroy(&leafSection);
829:   PetscSFDestroy(&labelSF);
830:   return(0);
831: }


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

839:   Not Collective

841:   Input Parameters:
842: + dm   - The DMPlex object
843: - name - The label name

845:   Level: intermediate

847: .keywords: mesh
848: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
849: @*/
850: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
851: {
852:   DM_Plex        *mesh = (DM_Plex*) dm->data;
853:   PlexLabel      next  = mesh->labels;
854:   PetscBool      flg   = PETSC_FALSE;

860:   while (next) {
861:     PetscStrcmp(name, next->label->name, &flg);
862:     if (flg) break;
863:     next = next->next;
864:   }
865:   if (!flg) {
866:     PlexLabel tmpLabel;

868:     PetscCalloc1(1, &tmpLabel);
869:     DMLabelCreate(name, &tmpLabel->label);
870:     tmpLabel->output = PETSC_TRUE;
871:     tmpLabel->next   = mesh->labels;
872:     mesh->labels     = tmpLabel;
873:   }
874:   return(0);
875: }

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

882:   Not Collective

884:   Input Parameters:
885: + dm   - The DMPlex object
886: . name - The label name
887: - point - The mesh point

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

892:   Level: beginner

894: .keywords: mesh
895: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
896: @*/
897: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
898: {
899:   DMLabel        label;

905:   DMPlexGetLabel(dm, name, &label);
906:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
907:   DMLabelGetValue(label, point, value);
908:   return(0);
909: }

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

916:   Not Collective

918:   Input Parameters:
919: + dm   - The DMPlex object
920: . name - The label name
921: . point - The mesh point
922: - value - The label value for this point

924:   Output Parameter:

926:   Level: beginner

928: .keywords: mesh
929: .seealso: DMLabelSetValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
930: @*/
931: PetscErrorCode DMPlexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
932: {
933:   DMLabel        label;

939:   DMPlexGetLabel(dm, name, &label);
940:   if (!label) {
941:     DMPlexCreateLabel(dm, name);
942:     DMPlexGetLabel(dm, name, &label);
943:   }
944:   DMLabelSetValue(label, point, value);
945:   return(0);
946: }

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

953:   Not Collective

955:   Input Parameters:
956: + dm   - The DMPlex object
957: . name - The label name
958: . point - The mesh point
959: - value - The label value for this point

961:   Output Parameter:

963:   Level: beginner

965: .keywords: mesh
966: .seealso: DMLabelClearValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
967: @*/
968: PetscErrorCode DMPlexClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
969: {
970:   DMLabel        label;

976:   DMPlexGetLabel(dm, name, &label);
977:   if (!label) return(0);
978:   DMLabelClearValue(label, point, value);
979:   return(0);
980: }

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

987:   Not Collective

989:   Input Parameters:
990: + dm   - The DMPlex object
991: - name - The label name

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

996:   Level: beginner

998: .keywords: mesh
999: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
1000: @*/
1001: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
1002: {
1003:   DMLabel        label;

1010:   DMPlexGetLabel(dm, name, &label);
1011:   *size = 0;
1012:   if (!label) return(0);
1013:   DMLabelGetNumValues(label, size);
1014:   return(0);
1015: }

1019: /*@C
1020:   DMPlexGetLabelIdIS - Get the integer ids in a label

1022:   Not Collective

1024:   Input Parameters:
1025: + mesh - The DMPlex object
1026: - name - The label name

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

1031:   Level: beginner

1033: .keywords: mesh
1034: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
1035: @*/
1036: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
1037: {
1038:   DMLabel        label;

1045:   DMPlexGetLabel(dm, name, &label);
1046:   *ids = NULL;
1047:   if (!label) return(0);
1048:   DMLabelGetValueIS(label, ids);
1049:   return(0);
1050: }

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

1057:   Not Collective

1059:   Input Parameters:
1060: + dm - The DMPlex object
1061: . name - The label name
1062: - value - The stratum value

1064:   Output Parameter:
1065: . size - The stratum size

1067:   Level: beginner

1069: .keywords: mesh
1070: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
1071: @*/
1072: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1073: {
1074:   DMLabel        label;

1081:   DMPlexGetLabel(dm, name, &label);
1082:   *size = 0;
1083:   if (!label) return(0);
1084:   DMLabelGetStratumSize(label, value, size);
1085:   return(0);
1086: }

1090: /*@C
1091:   DMPlexGetStratumIS - Get the points in a label stratum

1093:   Not Collective

1095:   Input Parameters:
1096: + dm - The DMPlex object
1097: . name - The label name
1098: - value - The stratum value

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

1103:   Level: beginner

1105: .keywords: mesh
1106: .seealso: DMLabelGetStratumIS(), DMPlexGetStratumSize()
1107: @*/
1108: PetscErrorCode DMPlexGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
1109: {
1110:   DMLabel        label;

1117:   DMPlexGetLabel(dm, name, &label);
1118:   *points = NULL;
1119:   if (!label) return(0);
1120:   DMLabelGetStratumIS(label, value, points);
1121:   return(0);
1122: }

1126: /*@C
1127:   DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label

1129:   Not Collective

1131:   Input Parameters:
1132: + dm   - The DMPlex object
1133: . name - The label name
1134: - value - The label value for this point

1136:   Output Parameter:

1138:   Level: beginner

1140: .keywords: mesh
1141: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1142: @*/
1143: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1144: {
1145:   DMLabel        label;

1151:   DMPlexGetLabel(dm, name, &label);
1152:   if (!label) return(0);
1153:   DMLabelClearStratum(label, value);
1154:   return(0);
1155: }

1159: /*@
1160:   DMPlexGetNumLabels - Return the number of labels defined by the mesh

1162:   Not Collective

1164:   Input Parameter:
1165: . dm   - The DMPlex object

1167:   Output Parameter:
1168: . numLabels - the number of Labels

1170:   Level: intermediate

1172: .keywords: mesh
1173: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1174: @*/
1175: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1176: {
1177:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1178:   PlexLabel next = mesh->labels;
1179:   PetscInt  n    = 0;

1184:   while (next) {++n; next = next->next;}
1185:   *numLabels = n;
1186:   return(0);
1187: }

1191: /*@C
1192:   DMPlexGetLabelName - Return the name of nth label

1194:   Not Collective

1196:   Input Parameters:
1197: + dm - The DMPlex object
1198: - n  - the label number

1200:   Output Parameter:
1201: . name - the label name

1203:   Level: intermediate

1205: .keywords: mesh
1206: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1207: @*/
1208: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1209: {
1210:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1211:   PlexLabel next = mesh->labels;
1212:   PetscInt  l    = 0;

1217:   while (next) {
1218:     if (l == n) {
1219:       *name = next->label->name;
1220:       return(0);
1221:     }
1222:     ++l;
1223:     next = next->next;
1224:   }
1225:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1226: }

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

1233:   Not Collective

1235:   Input Parameters:
1236: + dm   - The DMPlex object
1237: - name - The label name

1239:   Output Parameter:
1240: . hasLabel - PETSC_TRUE if the label is present

1242:   Level: intermediate

1244: .keywords: mesh
1245: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1246: @*/
1247: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1248: {
1249:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1250:   PlexLabel      next = mesh->labels;

1257:   *hasLabel = PETSC_FALSE;
1258:   while (next) {
1259:     PetscStrcmp(name, next->label->name, hasLabel);
1260:     if (*hasLabel) break;
1261:     next = next->next;
1262:   }
1263:   return(0);
1264: }

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

1271:   Not Collective

1273:   Input Parameters:
1274: + dm   - The DMPlex object
1275: - name - The label name

1277:   Output Parameter:
1278: . label - The DMLabel, or NULL if the label is absent

1280:   Level: intermediate

1282: .keywords: mesh
1283: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1284: @*/
1285: PetscErrorCode DMPlexGetLabel(DM dm, const char name[], DMLabel *label)
1286: {
1287:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1288:   PlexLabel      next = mesh->labels;
1289:   PetscBool      hasLabel;

1296:   *label = NULL;
1297:   while (next) {
1298:     PetscStrcmp(name, next->label->name, &hasLabel);
1299:     if (hasLabel) {
1300:       *label = next->label;
1301:       break;
1302:     }
1303:     next = next->next;
1304:   }
1305:   return(0);
1306: }

1310: /*@C
1311:   DMPlexGetLabelByNum - Return the nth label

1313:   Not Collective

1315:   Input Parameters:
1316: + dm - The DMPlex object
1317: - n  - the label number

1319:   Output Parameter:
1320: . label - the label

1322:   Level: intermediate

1324: .keywords: mesh
1325: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1326: @*/
1327: PetscErrorCode DMPlexGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
1328: {
1329:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1330:   PlexLabel next = mesh->labels;
1331:   PetscInt  l    = 0;

1336:   while (next) {
1337:     if (l == n) {
1338:       *label = next->label;
1339:       return(0);
1340:     }
1341:     ++l;
1342:     next = next->next;
1343:   }
1344:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1345: }

1349: /*@C
1350:   DMPlexAddLabel - Add the label to this mesh

1352:   Not Collective

1354:   Input Parameters:
1355: + dm   - The DMPlex object
1356: - label - The DMLabel

1358:   Level: developer

1360: .keywords: mesh
1361: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1362: @*/
1363: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1364: {
1365:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1366:   PlexLabel      tmpLabel;
1367:   PetscBool      hasLabel;

1372:   DMPlexHasLabel(dm, label->name, &hasLabel);
1373:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1374:   PetscCalloc1(1, &tmpLabel);
1375:   tmpLabel->label  = label;
1376:   tmpLabel->output = PETSC_TRUE;
1377:   tmpLabel->next   = mesh->labels;
1378:   mesh->labels     = tmpLabel;
1379:   return(0);
1380: }

1384: /*@C
1385:   DMPlexRemoveLabel - Remove the label from this mesh

1387:   Not Collective

1389:   Input Parameters:
1390: + dm   - The DMPlex object
1391: - name - The label name

1393:   Output Parameter:
1394: . label - The DMLabel, or NULL if the label is absent

1396:   Level: developer

1398: .keywords: mesh
1399: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1400: @*/
1401: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1402: {
1403:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1404:   PlexLabel      next = mesh->labels;
1405:   PlexLabel      last = NULL;
1406:   PetscBool      hasLabel;

1411:   DMPlexHasLabel(dm, name, &hasLabel);
1412:   *label = NULL;
1413:   if (!hasLabel) return(0);
1414:   while (next) {
1415:     PetscStrcmp(name, next->label->name, &hasLabel);
1416:     if (hasLabel) {
1417:       if (last) last->next   = next->next;
1418:       else      mesh->labels = next->next;
1419:       next->next = NULL;
1420:       *label     = next->label;
1421:       PetscFree(next);
1422:       break;
1423:     }
1424:     last = next;
1425:     next = next->next;
1426:   }
1427:   return(0);
1428: }

1432: /*@C
1433:   DMPlexGetLabelOutput - Get the output flag for a given label

1435:   Not Collective

1437:   Input Parameters:
1438: + dm   - The DMPlex object
1439: - name - The label name

1441:   Output Parameter:
1442: . output - The flag for output

1444:   Level: developer

1446: .keywords: mesh
1447: .seealso: DMPlexSetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1448: @*/
1449: PetscErrorCode DMPlexGetLabelOutput(DM dm, const char name[], PetscBool *output)
1450: {
1451:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1452:   PlexLabel      next = mesh->labels;

1459:   while (next) {
1460:     PetscBool flg;

1462:     PetscStrcmp(name, next->label->name, &flg);
1463:     if (flg) {*output = next->output; return(0);}
1464:     next = next->next;
1465:   }
1466:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1467: }

1471: /*@C
1472:   DMPlexSetLabelOutput - Set the output flag for a given label

1474:   Not Collective

1476:   Input Parameters:
1477: + dm     - The DMPlex object
1478: . name   - The label name
1479: - output - The flag for output

1481:   Level: developer

1483: .keywords: mesh
1484: .seealso: DMPlexGetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1485: @*/
1486: PetscErrorCode DMPlexSetLabelOutput(DM dm, const char name[], PetscBool output)
1487: {
1488:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1489:   PlexLabel      next = mesh->labels;

1495:   while (next) {
1496:     PetscBool flg;

1498:     PetscStrcmp(name, next->label->name, &flg);
1499:     if (flg) {next->output = output; return(0);}
1500:     next = next->next;
1501:   }
1502:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1503: }