Actual source code: plexlabel.c

petsc-master 2014-12-19
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     = PETSC_TRUE;
 18:   (*label)->stratumOffsets = NULL;
 19:   (*label)->stratumSizes   = NULL;
 20:   (*label)->points         = NULL;
 21:   (*label)->ht             = NULL;
 22:   (*label)->pStart         = -1;
 23:   (*label)->pEnd           = -1;
 24:   (*label)->bt             = NULL;
 25:   return(0);
 26: }

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

 35:   if (label->arrayValid) return 0;
 37:   PetscMalloc2(label->numStrata,&label->stratumSizes,label->numStrata+1,&label->stratumOffsets);
 38:   for (v = 0; v < label->numStrata; ++v) {
 39:     PetscInt size = 0;
 40:     PetscHashISize(label->ht[v], size);
 41:     label->stratumSizes[v]   = size;
 42:     label->stratumOffsets[v] = off;
 43:     off += size;
 44:   }
 45:   label->stratumOffsets[v] = off;
 46:   PetscMalloc1(off, &label->points);
 47:   off = 0;
 48:   for (v = 0; v < label->numStrata; ++v) {
 49:     PetscHashIGetKeys(label->ht[v], &off, label->points);
 50:     if (off != label->stratumOffsets[v+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of contributed points %d from value %d should be %d", off-label->stratumOffsets[v], label->stratumValues[v], label->stratumOffsets[v+1]-label->stratumOffsets[v]);
 51:     PetscHashIDestroy(label->ht[v]);
 52:     PetscSortInt(label->stratumSizes[v], &label->points[label->stratumOffsets[v]]);
 53:     if (label->bt) {
 54:       PetscInt p;

 56:       for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) {
 57:         const PetscInt point = label->points[p];

 59:         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);
 60:         PetscBTSet(label->bt, point - label->pStart);
 61:       }
 62:     }
 63:   }
 64:   PetscFree(label->ht);
 65:   label->ht = NULL;
 66:   label->arrayValid = PETSC_TRUE;
 67:   return(0);
 68: }

 72: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label)
 73: {
 74:   PetscInt       v;

 77:   if (!label->arrayValid) return 0;
 79:   PetscMalloc1(label->numStrata, &label->ht);
 80:   for (v = 0; v < label->numStrata; ++v) {
 81:     PETSC_UNUSED PetscHashIIter ret, iter;
 82:     PetscInt                    p;

 84:     PetscHashICreate(label->ht[v]);
 85:     for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[p], ret, iter);
 86:   }
 87:   PetscFree2(label->stratumSizes,label->stratumOffsets);
 88:   PetscFree(label->points);
 89:   label->arrayValid = PETSC_FALSE;
 90:   return(0);
 91: }

 95: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
 96: {
 99:   *name = label->name;
100:   return(0);
101: }

105: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
106: {
107:   PetscInt       v;
108:   PetscMPIInt    rank;

112:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
113:   if (label) {
114:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
115:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%d, %d)\n", label->pStart, label->pEnd);}
116:     for (v = 0; v < label->numStrata; ++v) {
117:       const PetscInt value = label->stratumValues[v];
118:       PetscInt       p;

120:       for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) {
121:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D (%D)\n", rank, label->points[p], value);
122:       }
123:     }
124:   }
125:   PetscViewerFlush(viewer);
126:   return(0);
127: }

131: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
132: {
133:   PetscBool      iascii;

138:   if (label) {DMLabelMakeValid_Private(label);}
139:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
140:   if (iascii) {
141:     DMLabelView_Ascii(label, viewer);
142:   } else SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not supported by this mesh object", ((PetscObject) viewer)->type_name);
143:   return(0);
144: }

148: PetscErrorCode DMLabelDestroy(DMLabel *label)
149: {

153:   if (!(*label)) return(0);
154:   if (--(*label)->refct > 0) return(0);
155:   PetscFree((*label)->name);
156:   PetscFree((*label)->stratumValues);
157:   PetscFree2((*label)->stratumSizes,(*label)->stratumOffsets);
158:   PetscFree((*label)->points);
159:   if ((*label)->ht) {
160:     PetscInt v;
161:     for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
162:     PetscFree((*label)->ht);
163:   }
164:   PetscBTDestroy(&(*label)->bt);
165:   PetscFree(*label);
166:   return(0);
167: }

171: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
172: {
173:   PetscInt       v, q;

177:   DMLabelMakeValid_Private(label);
178:   PetscNew(labelnew);
179:   PetscStrallocpy(label->name, &(*labelnew)->name);

181:   (*labelnew)->refct      = 1;
182:   (*labelnew)->numStrata  = label->numStrata;
183:   (*labelnew)->arrayValid = PETSC_TRUE;
184:   if (label->numStrata) {
185:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
186:     PetscMalloc2(label->numStrata,&(*labelnew)->stratumSizes,label->numStrata+1,&(*labelnew)->stratumOffsets);
187:     PetscMalloc1(label->stratumOffsets[label->numStrata], &(*labelnew)->points);
188:     /* Could eliminate unused space here */
189:     for (v = 0; v < label->numStrata; ++v) {
190:       (*labelnew)->stratumValues[v]  = label->stratumValues[v];
191:       (*labelnew)->stratumOffsets[v] = label->stratumOffsets[v];
192:       (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
193:       for (q = label->stratumOffsets[v]; q < label->stratumOffsets[v]+label->stratumSizes[v]; ++q) {
194:         (*labelnew)->points[q] = label->points[q];
195:       }
196:     }
197:     (*labelnew)->stratumOffsets[label->numStrata] = label->stratumOffsets[label->numStrata];
198:   }
199:   (*labelnew)->ht     = NULL;
200:   (*labelnew)->pStart = -1;
201:   (*labelnew)->pEnd   = -1;
202:   (*labelnew)->bt     = NULL;
203:   return(0);
204: }

208: /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
209: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
210: {
211:   PetscInt       v;

215:   DMLabelMakeValid_Private(label);
216:   if (label->bt) {PetscBTDestroy(&label->bt);}
217:   label->pStart = pStart;
218:   label->pEnd   = pEnd;
219:   PetscBTCreate(pEnd - pStart, &label->bt);
220:   PetscBTMemzero(pEnd - pStart, label->bt);
221:   for (v = 0; v < label->numStrata; ++v) {
222:     PetscInt i;

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

227:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, pStart, pEnd);
228:       PetscBTSet(label->bt, point - pStart);
229:     }
230:   }
231:   return(0);
232: }

236: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
237: {

241:   label->pStart = -1;
242:   label->pEnd   = -1;
243:   if (label->bt) {PetscBTDestroy(&label->bt);}
244:   return(0);
245: }

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

252:   Input Parameters:
253: + label - the DMLabel
254: - value - the value

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

259:   Level: developer

261: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
262: @*/
263: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
264: {
265:   PetscInt v;

269:   for (v = 0; v < label->numStrata; ++v) {
270:     if (value == label->stratumValues[v]) break;
271:   }
272:   *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
273:   return(0);
274: }

278: /*@
279:   DMLabelHasPoint - Determine whether a label assigns a value to a point

281:   Input Parameters:
282: + label - the DMLabel
283: - point - the point

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

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

290:   Level: developer

292: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
293: @*/
294: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
295: {

300:   DMLabelMakeValid_Private(label);
301: #if defined(PETSC_USE_DEBUG)
302:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
303:   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);
304: #endif
305:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
306:   return(0);
307: }

311: /*@
312:   DMLabelStratumHasPoint - Return true if the stratum contains a point

314:   Input Parameters:
315: + label - the DMLabel
316: . value - the stratum value
317: - point - the point

319:   Output Parameter:
320: . contains - true if the stratum contains the point

322:   Level: intermediate

324: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
325: @*/
326: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
327: {
328:   PetscInt       v;

333:   *contains = PETSC_FALSE;
334:   for (v = 0; v < label->numStrata; ++v) {
335:     if (label->stratumValues[v] == value) {
336:       if (label->arrayValid) {
337:         PetscInt i;

339:         PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &i);
340:         if (i >= 0) {
341:           *contains = PETSC_TRUE;
342:           break;
343:         }
344:       } else {
345:         PetscBool has;

347:         PetscHashIHasKey(label->ht[v], point, has);
348:         if (has) {
349:           *contains = PETSC_TRUE;
350:           break;
351:         }
352:       }
353:     }
354:   }
355:   return(0);
356: }

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

363:   Input Parameters:
364: + label - the DMLabel
365: - point - the point

367:   Output Parameter:
368: . value - The point value, or -1

370:   Level: intermediate

372: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
373: @*/
374: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
375: {
376:   PetscInt       v;

381:   *value = -1;
382:   for (v = 0; v < label->numStrata; ++v) {
383:     if (label->arrayValid) {
384:       PetscInt i;

386:       PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &i);
387:       if (i >= 0) {
388:         *value = label->stratumValues[v];
389:         break;
390:       }
391:     } else {
392:       PetscBool has;

394:       PetscHashIHasKey(label->ht[v], point, has);
395:       if (has) {
396:         *value = label->stratumValues[v];
397:         break;
398:       }
399:     }
400:   }
401:   return(0);
402: }

406: /*@
407:   DMLabelSetValue - Set the value a label assigns to a point

409:   Input Parameters:
410: + label - the DMLabel
411: . point - the point
412: - value - The point value

414:   Level: intermediate

416: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue()
417: @*/
418: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
419: {
420:   PETSC_UNUSED PetscHashIIter iter, ret;
421:   PetscInt                    v;
422:   PetscErrorCode              ierr;

425:   DMLabelMakeInvalid_Private(label);
426:   /* Find, or add, label value */
427:   for (v = 0; v < label->numStrata; ++v) {
428:     if (label->stratumValues[v] == value) break;
429:   }
430:   /* Create new table */
431:   if (v >= label->numStrata) {
432:     PetscInt   *tmpV;
433:     PetscHashI *tmpH;

435:     PetscMalloc1(label->numStrata+1, &tmpV);
436:     PetscMalloc1(label->numStrata+1, &tmpH);
437:     for (v = 0; v < label->numStrata; ++v) {
438:       tmpV[v] = label->stratumValues[v];
439:       tmpH[v] = label->ht[v];
440:     }
441:     tmpV[v] = value;
442:     PetscHashICreate(tmpH[v]);
443:     ++label->numStrata;
444:     PetscFree(label->stratumValues);
445:     PetscFree(label->ht);
446:     label->stratumValues = tmpV;
447:     label->ht            = tmpH;
448:   }
449:   /* Set key */
450:   PetscHashIPut(label->ht[v], point, ret, iter);
451:   return(0);
452: }

456: /*@
457:   DMLabelClearValue - Clear the value a label assigns to a point

459:   Input Parameters:
460: + label - the DMLabel
461: . point - the point
462: - value - The point value

464:   Level: intermediate

466: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
467: @*/
468: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
469: {
470:   PetscInt       v, p;

474:   /* Find label value */
475:   for (v = 0; v < label->numStrata; ++v) {
476:     if (label->stratumValues[v] == value) break;
477:   }
478:   if (v >= label->numStrata) return(0);
479:   if (label->arrayValid) {
480:     /* Check whether point exists */
481:     PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &p);
482:     if (p >= 0) {
483:       PetscMemmove(&label->points[p+label->stratumOffsets[v]], &label->points[p+label->stratumOffsets[v]+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
484:       --label->stratumSizes[v];
485:       if (label->bt) {
486:         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);
487:         PetscBTClear(label->bt, point - label->pStart);
488:       }
489:     }
490:   } else {
491:     PetscHashIDelKey(label->ht[v], point);
492:   }
493:   return(0);
494: }

498: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
499: {

504:   DMLabelMakeValid_Private(label);
505:   *numValues = label->numStrata;
506:   return(0);
507: }

511: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
512: {

517:   DMLabelMakeValid_Private(label);
518:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
519:   return(0);
520: }

524: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
525: {
526:   PetscInt       v;

531:   DMLabelMakeValid_Private(label);
532:   *size = 0;
533:   for (v = 0; v < label->numStrata; ++v) {
534:     if (label->stratumValues[v] == value) {
535:       *size = label->stratumSizes[v];
536:       break;
537:     }
538:   }
539:   return(0);
540: }

544: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
545: {
546:   PetscInt       v;

552:   DMLabelMakeValid_Private(label);
553:   for (v = 0; v < label->numStrata; ++v) {
554:     if (label->stratumValues[v] != value) continue;
555:     if (start) *start = label->points[label->stratumOffsets[v]];
556:     if (end)   *end   = label->points[label->stratumOffsets[v]+label->stratumSizes[v]-1]+1;
557:     break;
558:   }
559:   return(0);
560: }

564: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
565: {
566:   PetscInt       v;

571:   DMLabelMakeValid_Private(label);
572:   *points = NULL;
573:   for (v = 0; v < label->numStrata; ++v) {
574:     if (label->stratumValues[v] == value) {
575:       if (label->arrayValid) {
576:         ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], PETSC_COPY_VALUES, points);
577:         PetscObjectSetName((PetscObject) *points, "indices");
578:       } else {
579:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
580:       }
581:       break;
582:     }
583:   }
584:   return(0);
585: }

589: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
590: {
591:   PetscInt       v;

595:   for (v = 0; v < label->numStrata; ++v) {
596:     if (label->stratumValues[v] == value) break;
597:   }
598:   if (v >= label->numStrata) return(0);
599:   if (label->bt) {
600:     PetscInt i;

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

605:       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);
606:       PetscBTClear(label->bt, point - label->pStart);
607:     }
608:   }
609:   if (label->arrayValid) {
610:     label->stratumSizes[v] = 0;
611:   } else {
612:     PetscHashIClear(label->ht[v]);
613:   }
614:   return(0);
615: }

619: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
620: {
621:   PetscInt       v;

625:   DMLabelMakeValid_Private(label);
626:   label->pStart = start;
627:   label->pEnd   = end;
628:   if (label->bt) {PetscBTDestroy(&label->bt);}
629:   /* Could squish offsets, but would only make sense if I reallocate the storage */
630:   for (v = 0; v < label->numStrata; ++v) {
631:     const PetscInt offset = label->stratumOffsets[v];
632:     const PetscInt size   = label->stratumSizes[v];
633:     PetscInt       off    = offset, q;

635:     for (q = offset; q < offset+size; ++q) {
636:       const PetscInt point = label->points[q];

638:       if ((point < start) || (point >= end)) continue;
639:       label->points[off++] = point;
640:     }
641:     label->stratumSizes[v] = off-offset;
642:   }
643:   DMLabelCreateIndex(label, start, end);
644:   return(0);
645: }

649: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
650: {
651:   const PetscInt *perm;
652:   PetscInt        numValues, numPoints, v, q;
653:   PetscErrorCode  ierr;

656:   DMLabelMakeValid_Private(label);
657:   DMLabelDuplicate(label, labelNew);
658:   DMLabelGetNumValues(*labelNew, &numValues);
659:   ISGetLocalSize(permutation, &numPoints);
660:   ISGetIndices(permutation, &perm);
661:   for (v = 0; v < numValues; ++v) {
662:     const PetscInt offset = (*labelNew)->stratumOffsets[v];
663:     const PetscInt size   = (*labelNew)->stratumSizes[v];

665:     for (q = offset; q < offset+size; ++q) {
666:       const PetscInt point = (*labelNew)->points[q];

668:       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);
669:       (*labelNew)->points[q] = perm[point];
670:     }
671:     PetscSortInt(size, &(*labelNew)->points[offset]);
672:   }
673:   ISRestoreIndices(permutation, &perm);
674:   if (label->bt) {
675:     PetscBTDestroy(&label->bt);
676:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
677:   }
678:   return(0);
679: }

683: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
684: {
685:   MPI_Comm       comm;
686:   PetscSection   rootSection, leafSection;
687:   PetscSF        labelSF;
688:   PetscInt       p, pStart, pEnd, l, lStart, lEnd, s, nroots, nleaves, size, dof, offset, stratum;
689:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx, *leafStrata, *strataIdx;
690:   char          *name;
691:   PetscInt       nameSize;
692:   size_t         len = 0;
693:   PetscMPIInt    rank, numProcs;

697:   if (label) {DMLabelMakeValid_Private(label);}
698:   PetscObjectGetComm((PetscObject)sf, &comm);
699:   MPI_Comm_rank(comm, &rank);
700:   MPI_Comm_size(comm, &numProcs);
701:   /* Bcast name */
702:   if (!rank) {PetscStrlen(label->name, &len);}
703:   nameSize = len;
704:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
705:   PetscMalloc1(nameSize+1, &name);
706:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
707:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
708:   DMLabelCreate(name, labelNew);
709:   PetscFree(name);
710:   /* Bcast numStrata */
711:   if (!rank) (*labelNew)->numStrata = label->numStrata;
712:   MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
713:   /* Bcast stratumValues */
714:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
715:   if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
716:   MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);

718:   /* Build a section detailing strata-per-point, distribute and build SF
719:      from that and then send our points. */
720:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
721:   PetscSectionCreate(comm, &rootSection);
722:   PetscSectionSetChart(rootSection, 0, nroots);
723:   if (label) {
724:     for (s = 0; s < label->numStrata; ++s) {
725:       lStart = label->stratumOffsets[s];
726:       lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
727:       for (l=lStart; l<lEnd; l++) {
728:         PetscSectionGetDof(rootSection, label->points[l], &dof);
729:         PetscSectionSetDof(rootSection, label->points[l], dof+1);
730:       }
731:     }
732:   }
733:   PetscSectionSetUp(rootSection);

735:   /* Create a point-wise array of point strata */
736:   PetscSectionGetStorageSize(rootSection, &size);
737:   PetscMalloc1(size, &rootStrata);
738:   PetscCalloc1(nroots, &rootIdx);
739:   if (label) {
740:     for (s = 0; s < label->numStrata; ++s) {
741:       lStart = label->stratumOffsets[s];
742:       lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
743:       for (l=lStart; l<lEnd; l++) {
744:         p = label->points[l];
745:         PetscSectionGetOffset(rootSection, p, &offset);
746:         rootStrata[offset+rootIdx[p]++] = s;
747:       }
748:     }
749:   }

751:   /* Build SF that maps label points to remote processes */
752:   PetscSectionCreate(comm, &leafSection);
753:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, leafSection);
754:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, leafSection, &labelSF);

756:   /* Send the strata for each point over the derived SF */
757:   PetscSectionGetStorageSize(leafSection, &size);
758:   PetscMalloc1(size, &leafStrata);
759:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, leafStrata);
760:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, leafStrata);

762:   /* Rebuild the point strata on the receiver */
763:   PetscCalloc2((*labelNew)->numStrata,&(*labelNew)->stratumSizes,(*labelNew)->numStrata+1,&(*labelNew)->stratumOffsets);
764:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
765:   for (p=pStart; p<pEnd; p++) {
766:     PetscSectionGetDof(leafSection, p, &dof);
767:     PetscSectionGetOffset(leafSection, p, &offset);
768:     for (s=0; s<dof; s++) {
769:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
770:     }
771:   }
772:   /* Calculate stratumOffsets */
773:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
774:     (*labelNew)->stratumOffsets[s+1] = (*labelNew)->stratumSizes[s] + (*labelNew)->stratumOffsets[s];
775:   }
776:   PetscMalloc1((*labelNew)->stratumOffsets[(*labelNew)->numStrata],&(*labelNew)->points);

778:   /* Insert points into new strata */
779:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
780:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
781:   for (p=pStart; p<pEnd; p++) {
782:     PetscSectionGetDof(leafSection, p, &dof);
783:     PetscSectionGetOffset(leafSection, p, &offset);
784:     for (s=0; s<dof; s++) {
785:       stratum = leafStrata[offset+s];
786:       (*labelNew)->points[(*labelNew)->stratumOffsets[stratum] + strataIdx[stratum]++] = p;
787:     }
788:   }
789:   PetscFree2(rootStrata, leafStrata);
790:   PetscFree2(rootIdx, strataIdx);
791:   PetscSectionDestroy(&rootSection);
792:   PetscSectionDestroy(&leafSection);
793:   PetscSFDestroy(&labelSF);
794:   return(0);
795: }


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

803:   Not Collective

805:   Input Parameters:
806: + dm   - The DMPlex object
807: - name - The label name

809:   Level: intermediate

811: .keywords: mesh
812: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
813: @*/
814: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
815: {
816:   DM_Plex        *mesh = (DM_Plex*) dm->data;
817:   PlexLabel      next  = mesh->labels;
818:   PetscBool      flg   = PETSC_FALSE;

824:   while (next) {
825:     PetscStrcmp(name, next->label->name, &flg);
826:     if (flg) break;
827:     next = next->next;
828:   }
829:   if (!flg) {
830:     PlexLabel tmpLabel;

832:     PetscCalloc1(1, &tmpLabel);
833:     DMLabelCreate(name, &tmpLabel->label);
834:     tmpLabel->output = PETSC_TRUE;
835:     tmpLabel->next   = mesh->labels;
836:     mesh->labels     = tmpLabel;
837:   }
838:   return(0);
839: }

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

846:   Not Collective

848:   Input Parameters:
849: + dm   - The DMPlex object
850: . name - The label name
851: - point - The mesh point

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

856:   Level: beginner

858: .keywords: mesh
859: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
860: @*/
861: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
862: {
863:   DMLabel        label;

869:   DMPlexGetLabel(dm, name, &label);
870:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
871:   DMLabelGetValue(label, point, value);
872:   return(0);
873: }

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

880:   Not Collective

882:   Input Parameters:
883: + dm   - The DMPlex object
884: . name - The label name
885: . point - The mesh point
886: - value - The label value for this point

888:   Output Parameter:

890:   Level: beginner

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

903:   DMPlexGetLabel(dm, name, &label);
904:   if (!label) {
905:     DMPlexCreateLabel(dm, name);
906:     DMPlexGetLabel(dm, name, &label);
907:   }
908:   DMLabelSetValue(label, point, value);
909:   return(0);
910: }

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

917:   Not Collective

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

925:   Output Parameter:

927:   Level: beginner

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

940:   DMPlexGetLabel(dm, name, &label);
941:   if (!label) return(0);
942:   DMLabelClearValue(label, point, value);
943:   return(0);
944: }

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

951:   Not Collective

953:   Input Parameters:
954: + dm   - The DMPlex object
955: - name - The label name

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

960:   Level: beginner

962: .keywords: mesh
963: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
964: @*/
965: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
966: {
967:   DMLabel        label;

974:   DMPlexGetLabel(dm, name, &label);
975:   *size = 0;
976:   if (!label) return(0);
977:   DMLabelGetNumValues(label, size);
978:   return(0);
979: }

983: /*@C
984:   DMPlexGetLabelIdIS - Get the integer ids in a label

986:   Not Collective

988:   Input Parameters:
989: + mesh - The DMPlex object
990: - name - The label name

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

995:   Level: beginner

997: .keywords: mesh
998: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
999: @*/
1000: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
1001: {
1002:   DMLabel        label;

1009:   DMPlexGetLabel(dm, name, &label);
1010:   *ids = NULL;
1011:   if (!label) return(0);
1012:   DMLabelGetValueIS(label, ids);
1013:   return(0);
1014: }

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

1021:   Not Collective

1023:   Input Parameters:
1024: + dm - The DMPlex object
1025: . name - The label name
1026: - value - The stratum value

1028:   Output Parameter:
1029: . size - The stratum size

1031:   Level: beginner

1033: .keywords: mesh
1034: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
1035: @*/
1036: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1037: {
1038:   DMLabel        label;

1045:   DMPlexGetLabel(dm, name, &label);
1046:   *size = 0;
1047:   if (!label) return(0);
1048:   DMLabelGetStratumSize(label, value, size);
1049:   return(0);
1050: }

1054: /*@C
1055:   DMPlexGetStratumIS - Get the 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: . points - The stratum points, or NULL if the label does not exist or does not have that value

1067:   Level: beginner

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

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

1090: /*@C
1091:   DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label

1093:   Not Collective

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

1100:   Output Parameter:

1102:   Level: beginner

1104: .keywords: mesh
1105: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1106: @*/
1107: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1108: {
1109:   DMLabel        label;

1115:   DMPlexGetLabel(dm, name, &label);
1116:   if (!label) return(0);
1117:   DMLabelClearStratum(label, value);
1118:   return(0);
1119: }

1123: /*@
1124:   DMPlexGetNumLabels - Return the number of labels defined by the mesh

1126:   Not Collective

1128:   Input Parameter:
1129: . dm   - The DMPlex object

1131:   Output Parameter:
1132: . numLabels - the number of Labels

1134:   Level: intermediate

1136: .keywords: mesh
1137: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1138: @*/
1139: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1140: {
1141:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1142:   PlexLabel next = mesh->labels;
1143:   PetscInt  n    = 0;

1148:   while (next) {++n; next = next->next;}
1149:   *numLabels = n;
1150:   return(0);
1151: }

1155: /*@C
1156:   DMPlexGetLabelName - Return the name of nth label

1158:   Not Collective

1160:   Input Parameters:
1161: + dm - The DMPlex object
1162: - n  - the label number

1164:   Output Parameter:
1165: . name - the label name

1167:   Level: intermediate

1169: .keywords: mesh
1170: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1171: @*/
1172: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1173: {
1174:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1175:   PlexLabel next = mesh->labels;
1176:   PetscInt  l    = 0;

1181:   while (next) {
1182:     if (l == n) {
1183:       *name = next->label->name;
1184:       return(0);
1185:     }
1186:     ++l;
1187:     next = next->next;
1188:   }
1189:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1190: }

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

1197:   Not Collective

1199:   Input Parameters:
1200: + dm   - The DMPlex object
1201: - name - The label name

1203:   Output Parameter:
1204: . hasLabel - PETSC_TRUE if the label is present

1206:   Level: intermediate

1208: .keywords: mesh
1209: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1210: @*/
1211: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1212: {
1213:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1214:   PlexLabel      next = mesh->labels;

1221:   *hasLabel = PETSC_FALSE;
1222:   while (next) {
1223:     PetscStrcmp(name, next->label->name, hasLabel);
1224:     if (*hasLabel) break;
1225:     next = next->next;
1226:   }
1227:   return(0);
1228: }

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

1235:   Not Collective

1237:   Input Parameters:
1238: + dm   - The DMPlex object
1239: - name - The label name

1241:   Output Parameter:
1242: . label - The DMLabel, or NULL if the label is absent

1244:   Level: intermediate

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

1260:   *label = NULL;
1261:   while (next) {
1262:     PetscStrcmp(name, next->label->name, &hasLabel);
1263:     if (hasLabel) {
1264:       *label = next->label;
1265:       break;
1266:     }
1267:     next = next->next;
1268:   }
1269:   return(0);
1270: }

1274: /*@C
1275:   DMPlexGetLabelByNum - Return the nth label

1277:   Not Collective

1279:   Input Parameters:
1280: + dm - The DMPlex object
1281: - n  - the label number

1283:   Output Parameter:
1284: . label - the label

1286:   Level: intermediate

1288: .keywords: mesh
1289: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1290: @*/
1291: PetscErrorCode DMPlexGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
1292: {
1293:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1294:   PlexLabel next = mesh->labels;
1295:   PetscInt  l    = 0;

1300:   while (next) {
1301:     if (l == n) {
1302:       *label = next->label;
1303:       return(0);
1304:     }
1305:     ++l;
1306:     next = next->next;
1307:   }
1308:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1309: }

1313: /*@C
1314:   DMPlexAddLabel - Add the label to this mesh

1316:   Not Collective

1318:   Input Parameters:
1319: + dm   - The DMPlex object
1320: - label - The DMLabel

1322:   Level: developer

1324: .keywords: mesh
1325: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1326: @*/
1327: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1328: {
1329:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1330:   PlexLabel      tmpLabel;
1331:   PetscBool      hasLabel;

1336:   DMPlexHasLabel(dm, label->name, &hasLabel);
1337:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1338:   PetscCalloc1(1, &tmpLabel);
1339:   tmpLabel->label  = label;
1340:   tmpLabel->output = PETSC_TRUE;
1341:   tmpLabel->next   = mesh->labels;
1342:   mesh->labels     = tmpLabel;
1343:   return(0);
1344: }

1348: /*@C
1349:   DMPlexRemoveLabel - Remove the label from this mesh

1351:   Not Collective

1353:   Input Parameters:
1354: + dm   - The DMPlex object
1355: - name - The label name

1357:   Output Parameter:
1358: . label - The DMLabel, or NULL if the label is absent

1360:   Level: developer

1362: .keywords: mesh
1363: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1364: @*/
1365: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1366: {
1367:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1368:   PlexLabel      next = mesh->labels;
1369:   PlexLabel      last = NULL;
1370:   PetscBool      hasLabel;

1375:   DMPlexHasLabel(dm, name, &hasLabel);
1376:   *label = NULL;
1377:   if (!hasLabel) return(0);
1378:   while (next) {
1379:     PetscStrcmp(name, next->label->name, &hasLabel);
1380:     if (hasLabel) {
1381:       if (last) last->next   = next->next;
1382:       else      mesh->labels = next->next;
1383:       next->next = NULL;
1384:       *label     = next->label;
1385:       PetscFree(next);
1386:       break;
1387:     }
1388:     last = next;
1389:     next = next->next;
1390:   }
1391:   return(0);
1392: }

1396: /*@C
1397:   DMPlexGetLabelOutput - Get the output flag for a given label

1399:   Not Collective

1401:   Input Parameters:
1402: + dm   - The DMPlex object
1403: - name - The label name

1405:   Output Parameter:
1406: . output - The flag for output

1408:   Level: developer

1410: .keywords: mesh
1411: .seealso: DMPlexSetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1412: @*/
1413: PetscErrorCode DMPlexGetLabelOutput(DM dm, const char name[], PetscBool *output)
1414: {
1415:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1416:   PlexLabel      next = mesh->labels;

1423:   while (next) {
1424:     PetscBool flg;

1426:     PetscStrcmp(name, next->label->name, &flg);
1427:     if (flg) {*output = next->output; return(0);}
1428:     next = next->next;
1429:   }
1430:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1431: }

1435: /*@C
1436:   DMPlexSetLabelOutput - Set the output flag for a given label

1438:   Not Collective

1440:   Input Parameters:
1441: + dm     - The DMPlex object
1442: . name   - The label name
1443: - output - The flag for output

1445:   Level: developer

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

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

1462:     PetscStrcmp(name, next->label->name, &flg);
1463:     if (flg) {next->output = 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: }