Actual source code: dmlabel.c

petsc-master 2019-08-18
Report Typos and Errors
  1:  #include <petscdm.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/isimpl.h>
  4:  #include <petscsf.h>

  6: /*@C
  7:   DMLabelCreate - Create a DMLabel object, which is a multimap

  9:   Collective

 11:   Input parameters:
 12: + comm - The communicator, usually PETSC_COMM_SELF
 13: - name - The label name

 15:   Output parameter:
 16: . label - The DMLabel

 18:   Level: beginner

 20: .seealso: DMLabelDestroy()
 21: @*/
 22: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
 23: {

 28:   DMInitializePackage();

 30:   PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);

 32:   (*label)->numStrata      = 0;
 33:   (*label)->defaultValue   = -1;
 34:   (*label)->stratumValues  = NULL;
 35:   (*label)->validIS        = NULL;
 36:   (*label)->stratumSizes   = NULL;
 37:   (*label)->points         = NULL;
 38:   (*label)->ht             = NULL;
 39:   (*label)->pStart         = -1;
 40:   (*label)->pEnd           = -1;
 41:   (*label)->bt             = NULL;
 42:   PetscHMapICreate(&(*label)->hmap);
 43:   PetscObjectSetName((PetscObject) *label, name);
 44:   return(0);
 45: }

 47: /*
 48:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 50:   Not collective

 52:   Input parameter:
 53: + label - The DMLabel
 54: - v - The stratum value

 56:   Output parameter:
 57: . label - The DMLabel with stratum in sorted list format

 59:   Level: developer

 61: .seealso: DMLabelCreate()
 62: */
 63: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 64: {
 65:   IS             is;
 66:   PetscInt       off = 0, *pointArray, p;

 69:   if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
 71:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
 72:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
 73:   PetscMalloc1(label->stratumSizes[v], &pointArray);
 74:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
 75:   PetscHSetIClear(label->ht[v]);
 76:   PetscSortInt(label->stratumSizes[v], pointArray);
 77:   if (label->bt) {
 78:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 79:       const PetscInt point = pointArray[p];
 80:       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);
 81:       PetscBTSet(label->bt, point - label->pStart);
 82:     }
 83:   }
 84:   ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
 85:   PetscObjectSetName((PetscObject) is, "indices");
 86:   label->points[v]  = is;
 87:   label->validIS[v] = PETSC_TRUE;
 88:   PetscObjectStateIncrease((PetscObject) label);
 89:   return(0);
 90: }

 92: /*
 93:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

 95:   Not collective

 97:   Input parameter:
 98: . label - The DMLabel

100:   Output parameter:
101: . label - The DMLabel with all strata in sorted list format

103:   Level: developer

105: .seealso: DMLabelCreate()
106: */
107: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
108: {
109:   PetscInt       v;

113:   for (v = 0; v < label->numStrata; v++){
114:     DMLabelMakeValid_Private(label, v);
115:   }
116:   return(0);
117: }

119: /*
120:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

122:   Not collective

124:   Input parameter:
125: + label - The DMLabel
126: - v - The stratum value

128:   Output parameter:
129: . label - The DMLabel with stratum in hash format

131:   Level: developer

133: .seealso: DMLabelCreate()
134: */
135: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
136: {
137:   PetscInt       p;
138:   const PetscInt *points;

141:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
143:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
144:   if (label->points[v]) {
145:     ISGetIndices(label->points[v], &points);
146:     for (p = 0; p < label->stratumSizes[v]; ++p) {
147:       PetscHSetIAdd(label->ht[v], points[p]);
148:     }
149:     ISRestoreIndices(label->points[v],&points);
150:     ISDestroy(&(label->points[v]));
151:   }
152:   label->validIS[v] = PETSC_FALSE;
153:   return(0);
154: }

156: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
157: #define DMLABEL_LOOKUP_THRESHOLD 16
158: #endif

160: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
161: {
162:   PetscInt       v;

166:   *index = -1;
167:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
168:     for (v = 0; v < label->numStrata; ++v)
169:       if (label->stratumValues[v] == value) {*index = v; break;}
170:   } else {
171:     PetscHMapIGet(label->hmap, value, index);
172:   }
173: #if defined(PETSC_USE_DEBUG)
174:   { /* Check strata hash map consistency */
175:     PetscInt len, loc = -1;
176:     PetscHMapIGetSize(label->hmap, &len);
177:     if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
178:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
179:       PetscHMapIGet(label->hmap, value, &loc);
180:     } else {
181:       for (v = 0; v < label->numStrata; ++v)
182:         if (label->stratumValues[v] == value) {loc = v; break;}
183:     }
184:     if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
185:   }
186: #endif
187:   return(0);
188: }

190: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
191: {
192:   PetscInt       v;
193:   PetscInt      *tmpV;
194:   PetscInt      *tmpS;
195:   PetscHSetI    *tmpH, ht;
196:   IS            *tmpP, is;
197:   PetscBool     *tmpB;
198:   PetscHMapI     hmap = label->hmap;

202:   v    = label->numStrata;
203:   tmpV = label->stratumValues;
204:   tmpS = label->stratumSizes;
205:   tmpH = label->ht;
206:   tmpP = label->points;
207:   tmpB = label->validIS;
208:   { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free  */
209:     PetscInt   *oldV = tmpV;
210:     PetscInt   *oldS = tmpS;
211:     PetscHSetI *oldH = tmpH;
212:     IS         *oldP = tmpP;
213:     PetscBool  *oldB = tmpB;
214:     PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
215:     PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
216:     PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
217:     PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
218:     PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
219:     PetscArraycpy(tmpV, oldV, v);
220:     PetscArraycpy(tmpS, oldS, v);
221:     PetscArraycpy(tmpH, oldH, v);
222:     PetscArraycpy(tmpP, oldP, v);
223:     PetscArraycpy(tmpB, oldB, v);
224:     PetscFree(oldV);
225:     PetscFree(oldS);
226:     PetscFree(oldH);
227:     PetscFree(oldP);
228:     PetscFree(oldB);
229:   }
230:   label->numStrata     = v+1;
231:   label->stratumValues = tmpV;
232:   label->stratumSizes  = tmpS;
233:   label->ht            = tmpH;
234:   label->points        = tmpP;
235:   label->validIS       = tmpB;
236:   PetscHSetICreate(&ht);
237:   ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
238:   PetscHMapISet(hmap, value, v);
239:   tmpV[v] = value;
240:   tmpS[v] = 0;
241:   tmpH[v] = ht;
242:   tmpP[v] = is;
243:   tmpB[v] = PETSC_TRUE;
244:   PetscObjectStateIncrease((PetscObject) label);
245:   *index = v;
246:   return(0);
247: }

249: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
250: {
253:   DMLabelLookupStratum(label, value, index);
254:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
255:   return(0);
256: }

258: /*@
259:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

261:   Input Parameter:
262: + label - The DMLabel
263: - value - The stratum value

265:   Level: beginner

267: .seealso:  DMLabelCreate(), DMLabelDestroy()
268: @*/
269: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
270: {
271:   PetscInt       v;

276:   DMLabelLookupAddStratum(label, value, &v);
277:   return(0);
278: }

280: /*@
281:   DMLabelAddStrata - Adds new stratum values in a DMLabel

283:   Not collective

285:   Input Parameter:
286: + label - The DMLabel
287: . numStrata - The number of stratum values
288: - stratumValues - The stratum values

290:   Level: beginner

292: .seealso:  DMLabelCreate(), DMLabelDestroy()
293: @*/
294: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
295: {
296:   PetscInt       *values, v;

302:   PetscMalloc1(numStrata, &values);
303:   PetscArraycpy(values, stratumValues, numStrata);
304:   PetscSortRemoveDupsInt(&numStrata, values);
305:   if (!label->numStrata) { /* Fast preallocation */
306:     PetscInt   *tmpV;
307:     PetscInt   *tmpS;
308:     PetscHSetI *tmpH, ht;
309:     IS         *tmpP, is;
310:     PetscBool  *tmpB;
311:     PetscHMapI  hmap = label->hmap;

313:     PetscMalloc1(numStrata, &tmpV);
314:     PetscMalloc1(numStrata, &tmpS);
315:     PetscMalloc1(numStrata, &tmpH);
316:     PetscMalloc1(numStrata, &tmpP);
317:     PetscMalloc1(numStrata, &tmpB);
318:     label->numStrata     = numStrata;
319:     label->stratumValues = tmpV;
320:     label->stratumSizes  = tmpS;
321:     label->ht            = tmpH;
322:     label->points        = tmpP;
323:     label->validIS       = tmpB;
324:     for (v = 0; v < numStrata; ++v) {
325:       PetscHSetICreate(&ht);
326:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
327:       PetscHMapISet(hmap, values[v], v);
328:       tmpV[v] = values[v];
329:       tmpS[v] = 0;
330:       tmpH[v] = ht;
331:       tmpP[v] = is;
332:       tmpB[v] = PETSC_TRUE;
333:     }
334:     PetscObjectStateIncrease((PetscObject) label);
335:   } else {
336:     for (v = 0; v < numStrata; ++v) {
337:       DMLabelAddStratum(label, values[v]);
338:     }
339:   }
340:   PetscFree(values);
341:   return(0);
342: }

344: /*@
345:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

347:   Not collective

349:   Input Parameter:
350: + label - The DMLabel
351: - valueIS - Index set with stratum values

353:   Level: beginner

355: .seealso:  DMLabelCreate(), DMLabelDestroy()
356: @*/
357: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
358: {
359:   PetscInt       numStrata;
360:   const PetscInt *stratumValues;

366:   ISGetLocalSize(valueIS, &numStrata);
367:   ISGetIndices(valueIS, &stratumValues);
368:   DMLabelAddStrata(label, numStrata, stratumValues);
369:   return(0);
370: }

372: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
373: {
374:   PetscInt       v;
375:   PetscMPIInt    rank;

379:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
380:   PetscViewerASCIIPushSynchronized(viewer);
381:   if (label) {
382:     const char *name;

384:     PetscObjectGetName((PetscObject) label, &name);
385:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
386:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
387:     for (v = 0; v < label->numStrata; ++v) {
388:       const PetscInt value = label->stratumValues[v];
389:       const PetscInt *points;
390:       PetscInt       p;

392:       ISGetIndices(label->points[v], &points);
393:       for (p = 0; p < label->stratumSizes[v]; ++p) {
394:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
395:       }
396:       ISRestoreIndices(label->points[v],&points);
397:     }
398:   }
399:   PetscViewerFlush(viewer);
400:   PetscViewerASCIIPopSynchronized(viewer);
401:   return(0);
402: }

404: /*@C
405:   DMLabelView - View the label

407:   Collective on viewer

409:   Input Parameters:
410: + label - The DMLabel
411: - viewer - The PetscViewer

413:   Level: intermediate

415: .seealso: DMLabelCreate(), DMLabelDestroy()
416: @*/
417: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
418: {
419:   PetscBool      iascii;

424:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
426:   if (label) {DMLabelMakeAllValid_Private(label);}
427:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
428:   if (iascii) {
429:     DMLabelView_Ascii(label, viewer);
430:   }
431:   return(0);
432: }

434: /*@
435:   DMLabelReset - Destroys internal data structures in a DMLabel

437:   Not collective

439:   Input Parameter:
440: . label - The DMLabel

442:   Level: beginner

444: .seealso: DMLabelDestroy(), DMLabelCreate()
445: @*/
446: PetscErrorCode DMLabelReset(DMLabel label)
447: {
448:   PetscInt       v;

453:   for (v = 0; v < label->numStrata; ++v) {
454:     PetscHSetIDestroy(&label->ht[v]);
455:     ISDestroy(&label->points[v]);
456:   }
457:   label->numStrata = 0;
458:   PetscFree(label->stratumValues);
459:   PetscFree(label->stratumSizes);
460:   PetscFree(label->ht);
461:   PetscFree(label->points);
462:   PetscFree(label->validIS);
463:   PetscHMapIReset(label->hmap);
464:   label->pStart = -1;
465:   label->pEnd   = -1;
466:   PetscBTDestroy(&label->bt);
467:   return(0);
468: }

470: /*@
471:   DMLabelDestroy - Destroys a DMLabel

473:   Collective on label

475:   Input Parameter:
476: . label - The DMLabel

478:   Level: beginner

480: .seealso: DMLabelReset(), DMLabelCreate()
481: @*/
482: PetscErrorCode DMLabelDestroy(DMLabel *label)
483: {

487:   if (!*label) return(0);
489:   if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
490:   DMLabelReset(*label);
491:   PetscHMapIDestroy(&(*label)->hmap);
492:   PetscHeaderDestroy(label);
493:   return(0);
494: }

496: /*@
497:   DMLabelDuplicate - Duplicates a DMLabel

499:   Collective on label

501:   Input Parameter:
502: . label - The DMLabel

504:   Output Parameter:
505: . labelnew - location to put new vector

507:   Level: intermediate

509: .seealso: DMLabelCreate(), DMLabelDestroy()
510: @*/
511: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
512: {
513:   const char    *name;
514:   PetscInt       v;

519:   DMLabelMakeAllValid_Private(label);
520:   PetscObjectGetName((PetscObject) label, &name);
521:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

523:   (*labelnew)->numStrata    = label->numStrata;
524:   (*labelnew)->defaultValue = label->defaultValue;
525:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
526:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
527:   PetscMalloc1(label->numStrata, &(*labelnew)->ht);
528:   PetscMalloc1(label->numStrata, &(*labelnew)->points);
529:   PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
530:   for (v = 0; v < label->numStrata; ++v) {
531:     PetscHSetICreate(&(*labelnew)->ht[v]);
532:     (*labelnew)->stratumValues[v]  = label->stratumValues[v];
533:     (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
534:     PetscObjectReference((PetscObject) (label->points[v]));
535:     (*labelnew)->points[v]         = label->points[v];
536:     (*labelnew)->validIS[v]        = PETSC_TRUE;
537:   }
538:   PetscHMapIDestroy(&(*labelnew)->hmap);
539:   PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
540:   (*labelnew)->pStart = -1;
541:   (*labelnew)->pEnd   = -1;
542:   (*labelnew)->bt     = NULL;
543:   return(0);
544: }

546: /*@
547:   DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds

549:   Not collective

551:   Input Parameter:
552: . label  - The DMLabel

554:   Level: intermediate

556: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
557: @*/
558: PetscErrorCode DMLabelComputeIndex(DMLabel label)
559: {
560:   PetscInt       pStart = PETSC_MAX_INT, pEnd = -1, v;

565:   DMLabelMakeAllValid_Private(label);
566:   for (v = 0; v < label->numStrata; ++v) {
567:     const PetscInt *points;
568:     PetscInt       i;

570:     ISGetIndices(label->points[v], &points);
571:     for (i = 0; i < label->stratumSizes[v]; ++i) {
572:       const PetscInt point = points[i];

574:       pStart = PetscMin(point,   pStart);
575:       pEnd   = PetscMax(point+1, pEnd);
576:     }
577:     ISRestoreIndices(label->points[v], &points);
578:   }
579:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
580:   label->pEnd   = pEnd;
581:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
582:   return(0);
583: }

585: /*@
586:   DMLabelCreateIndex - Create an index structure for membership determination

588:   Not collective

590:   Input Parameters:
591: + label  - The DMLabel
592: . pStart - The smallest point
593: - pEnd   - The largest point + 1

595:   Level: intermediate

597: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
598: @*/
599: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
600: {
601:   PetscInt       v;

606:   DMLabelDestroyIndex(label);
607:   DMLabelMakeAllValid_Private(label);
608:   label->pStart = pStart;
609:   label->pEnd   = pEnd;
610:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
611:   PetscBTCreate(pEnd - pStart, &label->bt);
612:   for (v = 0; v < label->numStrata; ++v) {
613:     const PetscInt *points;
614:     PetscInt       i;

616:     ISGetIndices(label->points[v], &points);
617:     for (i = 0; i < label->stratumSizes[v]; ++i) {
618:       const PetscInt point = points[i];

620:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
621:       PetscBTSet(label->bt, point - pStart);
622:     }
623:     ISRestoreIndices(label->points[v], &points);
624:   }
625:   return(0);
626: }

628: /*@
629:   DMLabelDestroyIndex - Destroy the index structure

631:   Not collective

633:   Input Parameter:
634: . label - the DMLabel

636:   Level: intermediate

638: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
639: @*/
640: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
641: {

646:   label->pStart = -1;
647:   label->pEnd   = -1;
648:   PetscBTDestroy(&label->bt);
649:   return(0);
650: }

652: /*@
653:   DMLabelGetBounds - Return the smallest and largest point in the label

655:   Not collective

657:   Input Parameter:
658: . label - the DMLabel

660:   Output Parameters:
661: + pStart - The smallest point
662: - pEnd   - The largest point + 1

664:   Note: This will compute an index for the label if one does not exist.

666:   Level: intermediate

668: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
669: @*/
670: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
671: {

676:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
677:   if (pStart) {
679:     *pStart = label->pStart;
680:   }
681:   if (pEnd) {
683:     *pEnd = label->pEnd;
684:   }
685:   return(0);
686: }

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

691:   Not collective

693:   Input Parameters:
694: + label - the DMLabel
695: - value - the value

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

700:   Level: developer

702: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
703: @*/
704: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
705: {
706:   PetscInt v;

712:   DMLabelLookupStratum(label, value, &v);
713:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
714:   return(0);
715: }

717: /*@
718:   DMLabelHasPoint - Determine whether a label assigns a value to a point

720:   Not collective

722:   Input Parameters:
723: + label - the DMLabel
724: - point - the point

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

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

731:   Level: developer

733: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
734: @*/
735: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
736: {

742:   DMLabelMakeAllValid_Private(label);
743: #if defined(PETSC_USE_DEBUG)
744:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
745:   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);
746: #endif
747:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
748:   return(0);
749: }

751: /*@
752:   DMLabelStratumHasPoint - Return true if the stratum contains a point

754:   Not collective

756:   Input Parameters:
757: + label - the DMLabel
758: . value - the stratum value
759: - point - the point

761:   Output Parameter:
762: . contains - true if the stratum contains the point

764:   Level: intermediate

766: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
767: @*/
768: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
769: {
770:   PetscInt       v;

776:   *contains = PETSC_FALSE;
777:   DMLabelLookupStratum(label, value, &v);
778:   if (v < 0) return(0);

780:   if (label->validIS[v]) {
781:     PetscInt i;

783:     ISLocate(label->points[v], point, &i);
784:     if (i >= 0) *contains = PETSC_TRUE;
785:   } else {
786:     PetscBool has;

788:     PetscHSetIHas(label->ht[v], point, &has);
789:     if (has) *contains = PETSC_TRUE;
790:   }
791:   return(0);
792: }

794: /*@
795:   DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
796:   When a label is created, it is initialized to -1.

798:   Not collective

800:   Input parameter:
801: . label - a DMLabel object

803:   Output parameter:
804: . defaultValue - the default value

806:   Level: beginner

808: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
809: @*/
810: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
811: {
814:   *defaultValue = label->defaultValue;
815:   return(0);
816: }

818: /*@
819:   DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
820:   When a label is created, it is initialized to -1.

822:   Not collective

824:   Input parameter:
825: . label - a DMLabel object

827:   Output parameter:
828: . defaultValue - the default value

830:   Level: beginner

832: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
833: @*/
834: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
835: {
838:   label->defaultValue = defaultValue;
839:   return(0);
840: }

842: /*@
843:   DMLabelGetValue - Return the value a label assigns to a point, or the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue())

845:   Not collective

847:   Input Parameters:
848: + label - the DMLabel
849: - point - the point

851:   Output Parameter:
852: . value - The point value, or the default value (-1 by default)

854:   Level: intermediate

856: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
857: @*/
858: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
859: {
860:   PetscInt       v;

866:   *value = label->defaultValue;
867:   for (v = 0; v < label->numStrata; ++v) {
868:     if (label->validIS[v]) {
869:       PetscInt i;

871:       ISLocate(label->points[v], point, &i);
872:       if (i >= 0) {
873:         *value = label->stratumValues[v];
874:         break;
875:       }
876:     } else {
877:       PetscBool has;

879:       PetscHSetIHas(label->ht[v], point, &has);
880:       if (has) {
881:         *value = label->stratumValues[v];
882:         break;
883:       }
884:     }
885:   }
886:   return(0);
887: }

889: /*@
890:   DMLabelSetValue - Set the value a label assigns to a point.  If the value is the same as the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue() to something different), then this function will do nothing.

892:   Not collective

894:   Input Parameters:
895: + label - the DMLabel
896: . point - the point
897: - value - The point value

899:   Level: intermediate

901: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
902: @*/
903: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
904: {
905:   PetscInt       v;

910:   /* Find label value, add new entry if needed */
911:   if (value == label->defaultValue) return(0);
912:   DMLabelLookupAddStratum(label, value, &v);
913:   /* Set key */
914:   DMLabelMakeInvalid_Private(label, v);
915:   PetscHSetIAdd(label->ht[v], point);
916:   return(0);
917: }

919: /*@
920:   DMLabelClearValue - Clear the value a label assigns to a point

922:   Not collective

924:   Input Parameters:
925: + label - the DMLabel
926: . point - the point
927: - value - The point value

929:   Level: intermediate

931: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
932: @*/
933: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
934: {
935:   PetscInt       v;

940:   /* Find label value */
941:   DMLabelLookupStratum(label, value, &v);
942:   if (v < 0) return(0);

944:   if (label->bt) {
945:     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);
946:     PetscBTClear(label->bt, point - label->pStart);
947:   }

949:   /* Delete key */
950:   DMLabelMakeInvalid_Private(label, v);
951:   PetscHSetIDel(label->ht[v], point);
952:   return(0);
953: }

955: /*@
956:   DMLabelInsertIS - Set all points in the IS to a value

958:   Not collective

960:   Input Parameters:
961: + label - the DMLabel
962: . is    - the point IS
963: - value - The point value

965:   Level: intermediate

967: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
968: @*/
969: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
970: {
971:   PetscInt        v, n, p;
972:   const PetscInt *points;
973:   PetscErrorCode  ierr;

978:   /* Find label value, add new entry if needed */
979:   if (value == label->defaultValue) return(0);
980:   DMLabelLookupAddStratum(label, value, &v);
981:   /* Set keys */
982:   DMLabelMakeInvalid_Private(label, v);
983:   ISGetLocalSize(is, &n);
984:   ISGetIndices(is, &points);
985:   for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
986:   ISRestoreIndices(is, &points);
987:   return(0);
988: }

990: /*@
991:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

993:   Not collective

995:   Input Parameter:
996: . label - the DMLabel

998:   Output Paramater:
999: . numValues - the number of values

1001:   Level: intermediate

1003: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1004: @*/
1005: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1006: {
1010:   *numValues = label->numStrata;
1011:   return(0);
1012: }

1014: /*@
1015:   DMLabelGetValueIS - Get an IS of all values that the DMlabel takes

1017:   Not collective

1019:   Input Parameter:
1020: . label - the DMLabel

1022:   Output Paramater:
1023: . is    - the value IS

1025:   Level: intermediate

1027: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1028: @*/
1029: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1030: {

1036:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1037:   return(0);
1038: }

1040: /*@
1041:   DMLabelHasStratum - Determine whether points exist with the given value

1043:   Not collective

1045:   Input Parameters:
1046: + label - the DMLabel
1047: - value - the stratum value

1049:   Output Paramater:
1050: . exists - Flag saying whether points exist

1052:   Level: intermediate

1054: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1055: @*/
1056: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1057: {
1058:   PetscInt       v;

1064:   DMLabelLookupStratum(label, value, &v);
1065:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1066:   return(0);
1067: }

1069: /*@
1070:   DMLabelGetStratumSize - Get the size of a stratum

1072:   Not collective

1074:   Input Parameters:
1075: + label - the DMLabel
1076: - value - the stratum value

1078:   Output Paramater:
1079: . size - The number of points in the stratum

1081:   Level: intermediate

1083: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1084: @*/
1085: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1086: {
1087:   PetscInt       v;

1093:   *size = 0;
1094:   DMLabelLookupStratum(label, value, &v);
1095:   if (v < 0) return(0);
1096:   DMLabelMakeValid_Private(label, v);
1097:   *size = label->stratumSizes[v];
1098:   return(0);
1099: }

1101: /*@
1102:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1104:   Not collective

1106:   Input Parameters:
1107: + label - the DMLabel
1108: - value - the stratum value

1110:   Output Paramaters:
1111: + start - the smallest point in the stratum
1112: - end - the largest point in the stratum

1114:   Level: intermediate

1116: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1117: @*/
1118: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1119: {
1120:   PetscInt       v, min, max;

1127:   DMLabelLookupStratum(label, value, &v);
1128:   if (v < 0) return(0);
1129:   DMLabelMakeValid_Private(label, v);
1130:   if (label->stratumSizes[v] <= 0) return(0);
1131:   ISGetMinMax(label->points[v], &min, &max);
1132:   if (start) *start = min;
1133:   if (end)   *end   = max+1;
1134:   return(0);
1135: }

1137: /*@
1138:   DMLabelGetStratumIS - Get an IS with the stratum points

1140:   Not collective

1142:   Input Parameters:
1143: + label - the DMLabel
1144: - value - the stratum value

1146:   Output Paramater:
1147: . points - The stratum points

1149:   Level: intermediate

1151: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1152: @*/
1153: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1154: {
1155:   PetscInt       v;

1161:   *points = NULL;
1162:   DMLabelLookupStratum(label, value, &v);
1163:   if (v < 0) return(0);
1164:   DMLabelMakeValid_Private(label, v);
1165:   PetscObjectReference((PetscObject) label->points[v]);
1166:   *points = label->points[v];
1167:   return(0);
1168: }

1170: /*@
1171:   DMLabelSetStratumIS - Set the stratum points using an IS

1173:   Not collective

1175:   Input Parameters:
1176: + label - the DMLabel
1177: . value - the stratum value
1178: - points - The stratum points

1180:   Level: intermediate

1182: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1183: @*/
1184: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1185: {
1186:   PetscInt       v;

1192:   DMLabelLookupAddStratum(label, value, &v);
1193:   if (is == label->points[v]) return(0);
1194:   DMLabelClearStratum(label, value);
1195:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1196:   PetscObjectReference((PetscObject)is);
1197:   ISDestroy(&(label->points[v]));
1198:   label->points[v]  = is;
1199:   label->validIS[v] = PETSC_TRUE;
1200:   PetscObjectStateIncrease((PetscObject) label);
1201:   if (label->bt) {
1202:     const PetscInt *points;
1203:     PetscInt p;

1205:     ISGetIndices(is,&points);
1206:     for (p = 0; p < label->stratumSizes[v]; ++p) {
1207:       const PetscInt point = points[p];

1209:       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);
1210:       PetscBTSet(label->bt, point - label->pStart);
1211:     }
1212:   }
1213:   return(0);
1214: }

1216: /*@
1217:   DMLabelClearStratum - Remove a stratum

1219:   Not collective

1221:   Input Parameters:
1222: + label - the DMLabel
1223: - value - the stratum value

1225:   Level: intermediate

1227: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1228: @*/
1229: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1230: {
1231:   PetscInt       v;

1236:   DMLabelLookupStratum(label, value, &v);
1237:   if (v < 0) return(0);
1238:   if (label->validIS[v]) {
1239:     if (label->bt) {
1240:       PetscInt       i;
1241:       const PetscInt *points;

1243:       ISGetIndices(label->points[v], &points);
1244:       for (i = 0; i < label->stratumSizes[v]; ++i) {
1245:         const PetscInt point = points[i];

1247:         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);
1248:         PetscBTClear(label->bt, point - label->pStart);
1249:       }
1250:       ISRestoreIndices(label->points[v], &points);
1251:     }
1252:     label->stratumSizes[v] = 0;
1253:     ISDestroy(&label->points[v]);
1254:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1255:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1256:     PetscObjectStateIncrease((PetscObject) label);
1257:   } else {
1258:     PetscHSetIClear(label->ht[v]);
1259:   }
1260:   return(0);
1261: }

1263: /*@
1264:   DMLabelFilter - Remove all points outside of [start, end)

1266:   Not collective

1268:   Input Parameters:
1269: + label - the DMLabel
1270: . start - the first point
1271: - end - the last point

1273:   Level: intermediate

1275: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1276: @*/
1277: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1278: {
1279:   PetscInt       v;

1284:   DMLabelDestroyIndex(label);
1285:   DMLabelMakeAllValid_Private(label);
1286:   for (v = 0; v < label->numStrata; ++v) {
1287:     PetscInt off, q;
1288:     const PetscInt *points;
1289:     PetscInt numPointsNew = 0, *pointsNew = NULL;

1291:     ISGetIndices(label->points[v], &points);
1292:     for (q = 0; q < label->stratumSizes[v]; ++q)
1293:       if (points[q] >= start && points[q] < end)
1294:         numPointsNew++;
1295:     PetscMalloc1(numPointsNew, &pointsNew);
1296:     for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
1297:       if (points[q] >= start && points[q] < end)
1298:         pointsNew[off++] = points[q];
1299:     }
1300:     ISRestoreIndices(label->points[v],&points);

1302:     label->stratumSizes[v] = numPointsNew;
1303:     ISDestroy(&label->points[v]);
1304:     ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
1305:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1306:   }
1307:   DMLabelCreateIndex(label, start, end);
1308:   return(0);
1309: }

1311: /*@
1312:   DMLabelPermute - Create a new label with permuted points

1314:   Not collective

1316:   Input Parameters:
1317: + label - the DMLabel
1318: - permutation - the point permutation

1320:   Output Parameter:
1321: . labelnew - the new label containing the permuted points

1323:   Level: intermediate

1325: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1326: @*/
1327: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1328: {
1329:   const PetscInt *perm;
1330:   PetscInt        numValues, numPoints, v, q;
1331:   PetscErrorCode  ierr;

1336:   DMLabelMakeAllValid_Private(label);
1337:   DMLabelDuplicate(label, labelNew);
1338:   DMLabelGetNumValues(*labelNew, &numValues);
1339:   ISGetLocalSize(permutation, &numPoints);
1340:   ISGetIndices(permutation, &perm);
1341:   for (v = 0; v < numValues; ++v) {
1342:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1343:     const PetscInt *points;
1344:     PetscInt *pointsNew;

1346:     ISGetIndices((*labelNew)->points[v],&points);
1347:     PetscMalloc1(size,&pointsNew);
1348:     for (q = 0; q < size; ++q) {
1349:       const PetscInt point = points[q];

1351:       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);
1352:       pointsNew[q] = perm[point];
1353:     }
1354:     ISRestoreIndices((*labelNew)->points[v],&points);
1355:     PetscSortInt(size, pointsNew);
1356:     ISDestroy(&((*labelNew)->points[v]));
1357:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1358:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1359:       PetscFree(pointsNew);
1360:     } else {
1361:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1362:     }
1363:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1364:   }
1365:   ISRestoreIndices(permutation, &perm);
1366:   if (label->bt) {
1367:     PetscBTDestroy(&label->bt);
1368:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1369:   }
1370:   return(0);
1371: }

1373: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1374: {
1375:   MPI_Comm       comm;
1376:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1377:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1378:   PetscSection   rootSection;
1379:   PetscSF        labelSF;

1383:   if (label) {DMLabelMakeAllValid_Private(label);}
1384:   PetscObjectGetComm((PetscObject)sf, &comm);
1385:   /* Build a section of stratum values per point, generate the according SF
1386:      and distribute point-wise stratum values to leaves. */
1387:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1388:   PetscSectionCreate(comm, &rootSection);
1389:   PetscSectionSetChart(rootSection, 0, nroots);
1390:   if (label) {
1391:     for (s = 0; s < label->numStrata; ++s) {
1392:       const PetscInt *points;

1394:       ISGetIndices(label->points[s], &points);
1395:       for (l = 0; l < label->stratumSizes[s]; l++) {
1396:         PetscSectionGetDof(rootSection, points[l], &dof);
1397:         PetscSectionSetDof(rootSection, points[l], dof+1);
1398:       }
1399:       ISRestoreIndices(label->points[s], &points);
1400:     }
1401:   }
1402:   PetscSectionSetUp(rootSection);
1403:   /* Create a point-wise array of stratum values */
1404:   PetscSectionGetStorageSize(rootSection, &size);
1405:   PetscMalloc1(size, &rootStrata);
1406:   PetscCalloc1(nroots, &rootIdx);
1407:   if (label) {
1408:     for (s = 0; s < label->numStrata; ++s) {
1409:       const PetscInt *points;

1411:       ISGetIndices(label->points[s], &points);
1412:       for (l = 0; l < label->stratumSizes[s]; l++) {
1413:         const PetscInt p = points[l];
1414:         PetscSectionGetOffset(rootSection, p, &offset);
1415:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1416:       }
1417:       ISRestoreIndices(label->points[s], &points);
1418:     }
1419:   }
1420:   /* Build SF that maps label points to remote processes */
1421:   PetscSectionCreate(comm, leafSection);
1422:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1423:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1424:   PetscFree(remoteOffsets);
1425:   /* Send the strata for each point over the derived SF */
1426:   PetscSectionGetStorageSize(*leafSection, &size);
1427:   PetscMalloc1(size, leafStrata);
1428:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1429:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1430:   /* Clean up */
1431:   PetscFree(rootStrata);
1432:   PetscFree(rootIdx);
1433:   PetscSectionDestroy(&rootSection);
1434:   PetscSFDestroy(&labelSF);
1435:   return(0);
1436: }

1438: /*@
1439:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1441:   Collective on sf

1443:   Input Parameters:
1444: + label - the DMLabel
1445: - sf    - the map from old to new distribution

1447:   Output Parameter:
1448: . labelnew - the new redistributed label

1450:   Level: intermediate

1452: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1453: @*/
1454: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1455: {
1456:   MPI_Comm       comm;
1457:   PetscSection   leafSection;
1458:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1459:   PetscInt      *leafStrata, *strataIdx;
1460:   PetscInt     **points;
1461:   const char    *lname = NULL;
1462:   char          *name;
1463:   PetscInt       nameSize;
1464:   PetscHSetI     stratumHash;
1465:   size_t         len = 0;
1466:   PetscMPIInt    rank;

1471:   if (label) {
1473:     DMLabelMakeAllValid_Private(label);
1474:   }
1475:   PetscObjectGetComm((PetscObject)sf, &comm);
1476:   MPI_Comm_rank(comm, &rank);
1477:   /* Bcast name */
1478:   if (!rank) {
1479:     PetscObjectGetName((PetscObject) label, &lname);
1480:     PetscStrlen(lname, &len);
1481:   }
1482:   nameSize = len;
1483:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1484:   PetscMalloc1(nameSize+1, &name);
1485:   if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1486:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1487:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1488:   PetscFree(name);
1489:   /* Bcast defaultValue */
1490:   if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1491:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1492:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1493:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1494:   /* Determine received stratum values and initialise new label*/
1495:   PetscHSetICreate(&stratumHash);
1496:   PetscSectionGetStorageSize(leafSection, &size);
1497:   for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1498:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1499:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1500:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1501:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1502:   /* Turn leafStrata into indices rather than stratum values */
1503:   offset = 0;
1504:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1505:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1506:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1507:     PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1508:   }
1509:   for (p = 0; p < size; ++p) {
1510:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1511:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1512:     }
1513:   }
1514:   /* Rebuild the point strata on the receiver */
1515:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1516:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1517:   for (p=pStart; p<pEnd; p++) {
1518:     PetscSectionGetDof(leafSection, p, &dof);
1519:     PetscSectionGetOffset(leafSection, p, &offset);
1520:     for (s=0; s<dof; s++) {
1521:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1522:     }
1523:   }
1524:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1525:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1526:   PetscMalloc1((*labelNew)->numStrata,&points);
1527:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1528:     PetscHSetICreate(&(*labelNew)->ht[s]);
1529:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1530:   }
1531:   /* Insert points into new strata */
1532:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1533:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1534:   for (p=pStart; p<pEnd; p++) {
1535:     PetscSectionGetDof(leafSection, p, &dof);
1536:     PetscSectionGetOffset(leafSection, p, &offset);
1537:     for (s=0; s<dof; s++) {
1538:       stratum = leafStrata[offset+s];
1539:       points[stratum][strataIdx[stratum]++] = p;
1540:     }
1541:   }
1542:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1543:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1544:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1545:   }
1546:   PetscFree(points);
1547:   PetscHSetIDestroy(&stratumHash);
1548:   PetscFree(leafStrata);
1549:   PetscFree(strataIdx);
1550:   PetscSectionDestroy(&leafSection);
1551:   return(0);
1552: }

1554: /*@
1555:   DMLabelGather - Gather all label values from leafs into roots

1557:   Collective on sf

1559:   Input Parameters:
1560: + label - the DMLabel
1561: - sf - the Star Forest point communication map

1563:   Output Parameters:
1564: . labelNew - the new DMLabel with localised leaf values

1566:   Level: developer

1568:   Note: This is the inverse operation to DMLabelDistribute.

1570: .seealso: DMLabelDistribute()
1571: @*/
1572: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1573: {
1574:   MPI_Comm       comm;
1575:   PetscSection   rootSection;
1576:   PetscSF        sfLabel;
1577:   PetscSFNode   *rootPoints, *leafPoints;
1578:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1579:   const PetscInt *rootDegree, *ilocal;
1580:   PetscInt       *rootStrata;
1581:   const char    *lname;
1582:   char          *name;
1583:   PetscInt       nameSize;
1584:   size_t         len = 0;
1585:   PetscMPIInt    rank, size;

1591:   PetscObjectGetComm((PetscObject)sf, &comm);
1592:   MPI_Comm_rank(comm, &rank);
1593:   MPI_Comm_size(comm, &size);
1594:   /* Bcast name */
1595:   if (!rank) {
1596:     PetscObjectGetName((PetscObject) label, &lname);
1597:     PetscStrlen(lname, &len);
1598:   }
1599:   nameSize = len;
1600:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1601:   PetscMalloc1(nameSize+1, &name);
1602:   if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1603:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1604:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1605:   PetscFree(name);
1606:   /* Gather rank/index pairs of leaves into local roots to build
1607:      an inverse, multi-rooted SF. Note that this ignores local leaf
1608:      indexing due to the use of the multiSF in PetscSFGather. */
1609:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1610:   PetscMalloc1(nroots, &leafPoints);
1611:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1612:   for (p = 0; p < nleaves; p++) {
1613:     PetscInt ilp = ilocal ? ilocal[p] : p;

1615:     leafPoints[ilp].index = ilp;
1616:     leafPoints[ilp].rank  = rank;
1617:   }
1618:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1619:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1620:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1621:   PetscMalloc1(nmultiroots, &rootPoints);
1622:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1623:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1624:   PetscSFCreate(comm,& sfLabel);
1625:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1626:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1627:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1628:   /* Rebuild the point strata on the receiver */
1629:   for (p = 0, idx = 0; p < nroots; p++) {
1630:     for (d = 0; d < rootDegree[p]; d++) {
1631:       PetscSectionGetDof(rootSection, idx+d, &dof);
1632:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1633:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1634:     }
1635:     idx += rootDegree[p];
1636:   }
1637:   PetscFree(leafPoints);
1638:   PetscFree(rootStrata);
1639:   PetscSectionDestroy(&rootSection);
1640:   PetscSFDestroy(&sfLabel);
1641:   return(0);
1642: }

1644: /*@
1645:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1647:   Not collective

1649:   Input Parameter:
1650: . label - the DMLabel

1652:   Output Parameters:
1653: + section - the section giving offsets for each stratum
1654: - is - An IS containing all the label points

1656:   Level: developer

1658: .seealso: DMLabelDistribute()
1659: @*/
1660: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1661: {
1662:   IS              vIS;
1663:   const PetscInt *values;
1664:   PetscInt       *points;
1665:   PetscInt        nV, vS = 0, vE = 0, v, N;
1666:   PetscErrorCode  ierr;

1670:   DMLabelGetNumValues(label, &nV);
1671:   DMLabelGetValueIS(label, &vIS);
1672:   ISGetIndices(vIS, &values);
1673:   if (nV) {vS = values[0]; vE = values[0]+1;}
1674:   for (v = 1; v < nV; ++v) {
1675:     vS = PetscMin(vS, values[v]);
1676:     vE = PetscMax(vE, values[v]+1);
1677:   }
1678:   PetscSectionCreate(PETSC_COMM_SELF, section);
1679:   PetscSectionSetChart(*section, vS, vE);
1680:   for (v = 0; v < nV; ++v) {
1681:     PetscInt n;

1683:     DMLabelGetStratumSize(label, values[v], &n);
1684:     PetscSectionSetDof(*section, values[v], n);
1685:   }
1686:   PetscSectionSetUp(*section);
1687:   PetscSectionGetStorageSize(*section, &N);
1688:   PetscMalloc1(N, &points);
1689:   for (v = 0; v < nV; ++v) {
1690:     IS              is;
1691:     const PetscInt *spoints;
1692:     PetscInt        dof, off, p;

1694:     PetscSectionGetDof(*section, values[v], &dof);
1695:     PetscSectionGetOffset(*section, values[v], &off);
1696:     DMLabelGetStratumIS(label, values[v], &is);
1697:     ISGetIndices(is, &spoints);
1698:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1699:     ISRestoreIndices(is, &spoints);
1700:     ISDestroy(&is);
1701:   }
1702:   ISRestoreIndices(vIS, &values);
1703:   ISDestroy(&vIS);
1704:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1705:   return(0);
1706: }

1708: /*@
1709:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1710:   the local section and an SF describing the section point overlap.

1712:   Collective on sf

1714:   Input Parameters:
1715:   + s - The PetscSection for the local field layout
1716:   . sf - The SF describing parallel layout of the section points
1717:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1718:   . label - The label specifying the points
1719:   - labelValue - The label stratum specifying the points

1721:   Output Parameter:
1722:   . gsection - The PetscSection for the global field layout

1724:   Note: This gives negative sizes and offsets to points not owned by this process

1726:   Level: developer

1728: .seealso: PetscSectionCreate()
1729: @*/
1730: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1731: {
1732:   PetscInt      *neg = NULL, *tmpOff = NULL;
1733:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1740:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1741:   PetscSectionGetChart(s, &pStart, &pEnd);
1742:   PetscSectionSetChart(*gsection, pStart, pEnd);
1743:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1744:   if (nroots >= 0) {
1745:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1746:     PetscCalloc1(nroots, &neg);
1747:     if (nroots > pEnd-pStart) {
1748:       PetscCalloc1(nroots, &tmpOff);
1749:     } else {
1750:       tmpOff = &(*gsection)->atlasDof[-pStart];
1751:     }
1752:   }
1753:   /* Mark ghost points with negative dof */
1754:   for (p = pStart; p < pEnd; ++p) {
1755:     PetscInt value;

1757:     DMLabelGetValue(label, p, &value);
1758:     if (value != labelValue) continue;
1759:     PetscSectionGetDof(s, p, &dof);
1760:     PetscSectionSetDof(*gsection, p, dof);
1761:     PetscSectionGetConstraintDof(s, p, &cdof);
1762:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1763:     if (neg) neg[p] = -(dof+1);
1764:   }
1765:   PetscSectionSetUpBC(*gsection);
1766:   if (nroots >= 0) {
1767:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1768:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1769:     if (nroots > pEnd-pStart) {
1770:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1771:     }
1772:   }
1773:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1774:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1775:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1776:     (*gsection)->atlasOff[p] = off;
1777:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1778:   }
1779:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1780:   globalOff -= off;
1781:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1782:     (*gsection)->atlasOff[p] += globalOff;
1783:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1784:   }
1785:   /* Put in negative offsets for ghost points */
1786:   if (nroots >= 0) {
1787:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1788:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1789:     if (nroots > pEnd-pStart) {
1790:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1791:     }
1792:   }
1793:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1794:   PetscFree(neg);
1795:   return(0);
1796: }

1798: typedef struct _n_PetscSectionSym_Label
1799: {
1800:   DMLabel           label;
1801:   PetscCopyMode     *modes;
1802:   PetscInt          *sizes;
1803:   const PetscInt    ***perms;
1804:   const PetscScalar ***rots;
1805:   PetscInt          (*minMaxOrients)[2];
1806:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1807: } PetscSectionSym_Label;

1809: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1810: {
1811:   PetscInt              i, j;
1812:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1813:   PetscErrorCode        ierr;

1816:   for (i = 0; i <= sl->numStrata; i++) {
1817:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1818:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1819:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1820:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1821:       }
1822:       if (sl->perms[i]) {
1823:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1825:         PetscFree(perms);
1826:       }
1827:       if (sl->rots[i]) {
1828:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1830:         PetscFree(rots);
1831:       }
1832:     }
1833:   }
1834:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1835:   DMLabelDestroy(&sl->label);
1836:   sl->numStrata = 0;
1837:   return(0);
1838: }

1840: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1841: {

1845:   PetscSectionSymLabelReset(sym);
1846:   PetscFree(sym->data);
1847:   return(0);
1848: }

1850: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1851: {
1852:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1853:   PetscBool             isAscii;
1854:   DMLabel               label = sl->label;
1855:   const char           *name;
1856:   PetscErrorCode        ierr;

1859:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1860:   if (isAscii) {
1861:     PetscInt          i, j, k;
1862:     PetscViewerFormat format;

1864:     PetscViewerGetFormat(viewer,&format);
1865:     if (label) {
1866:       PetscViewerGetFormat(viewer,&format);
1867:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1868:         PetscViewerASCIIPushTab(viewer);
1869:         DMLabelView(label, viewer);
1870:         PetscViewerASCIIPopTab(viewer);
1871:       } else {
1872:         PetscObjectGetName((PetscObject) sl->label, &name);
1873:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
1874:       }
1875:     } else {
1876:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1877:     }
1878:     PetscViewerASCIIPushTab(viewer);
1879:     for (i = 0; i <= sl->numStrata; i++) {
1880:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1882:       if (!(sl->perms[i] || sl->rots[i])) {
1883:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1884:       } else {
1885:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1886:         PetscViewerASCIIPushTab(viewer);
1887:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1888:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1889:           PetscViewerASCIIPushTab(viewer);
1890:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1891:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1892:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1893:             } else {
1894:               PetscInt tab;

1896:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1897:               PetscViewerASCIIPushTab(viewer);
1898:               PetscViewerASCIIGetTab(viewer,&tab);
1899:               if (sl->perms[i] && sl->perms[i][j]) {
1900:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1901:                 PetscViewerASCIISetTab(viewer,0);
1902:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1903:                 PetscViewerASCIIPrintf(viewer,"\n");
1904:                 PetscViewerASCIISetTab(viewer,tab);
1905:               }
1906:               if (sl->rots[i] && sl->rots[i][j]) {
1907:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
1908:                 PetscViewerASCIISetTab(viewer,0);
1909: #if defined(PETSC_USE_COMPLEX)
1910:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f+i*%+f",PetscRealPart(sl->rots[i][j][k]),PetscImaginaryPart(sl->rots[i][j][k]));}
1911: #else
1912:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1913: #endif
1914:                 PetscViewerASCIIPrintf(viewer,"\n");
1915:                 PetscViewerASCIISetTab(viewer,tab);
1916:               }
1917:               PetscViewerASCIIPopTab(viewer);
1918:             }
1919:           }
1920:           PetscViewerASCIIPopTab(viewer);
1921:         }
1922:         PetscViewerASCIIPopTab(viewer);
1923:       }
1924:     }
1925:     PetscViewerASCIIPopTab(viewer);
1926:   }
1927:   return(0);
1928: }

1930: /*@
1931:   PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries

1933:   Logically collective on sym

1935:   Input parameters:
1936: + sym - the section symmetries
1937: - label - the DMLabel describing the types of points

1939:   Level: developer:

1941: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1942: @*/
1943: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1944: {
1945:   PetscSectionSym_Label *sl;
1946:   PetscErrorCode        ierr;

1950:   sl = (PetscSectionSym_Label *) sym->data;
1951:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1952:   if (label) {
1953:     sl->label = label;
1954:     PetscObjectReference((PetscObject) label);
1955:     DMLabelGetNumValues(label,&sl->numStrata);
1956:     PetscMalloc5(sl->numStrata+1,&sl->modes,sl->numStrata+1,&sl->sizes,sl->numStrata+1,&sl->perms,sl->numStrata+1,&sl->rots,sl->numStrata+1,&sl->minMaxOrients);
1957:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1958:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1959:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1960:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1961:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1962:   }
1963:   return(0);
1964: }

1966: /*@C
1967:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

1969:   Logically collective on sym

1971:   InputParameters:
1972: + sym       - the section symmetries
1973: . stratum   - the stratum value in the label that we are assigning symmetries for
1974: . size      - the number of dofs for points in the stratum of the label
1975: . minOrient - the smallest orientation for a point in this stratum
1976: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1977: . mode      - how sym should copy the perms and rots arrays
1978: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
1979: - rots      - NULL if there are no rotations, or (maxOrient - minOrient) sets of rotations, one for each orientation.  A NULL set of orientations is the identity

1981:   Level: developer

1983: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1984: @*/
1985: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1986: {
1987:   PetscSectionSym_Label *sl;
1988:   const char            *name;
1989:   PetscInt               i, j, k;
1990:   PetscErrorCode         ierr;

1994:   sl = (PetscSectionSym_Label *) sym->data;
1995:   if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1996:   for (i = 0; i <= sl->numStrata; i++) {
1997:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

1999:     if (stratum == value) break;
2000:   }
2001:   PetscObjectGetName((PetscObject) sl->label, &name);
2002:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2003:   sl->sizes[i] = size;
2004:   sl->modes[i] = mode;
2005:   sl->minMaxOrients[i][0] = minOrient;
2006:   sl->minMaxOrients[i][1] = maxOrient;
2007:   if (mode == PETSC_COPY_VALUES) {
2008:     if (perms) {
2009:       PetscInt    **ownPerms;

2011:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
2012:       for (j = 0; j < maxOrient-minOrient; j++) {
2013:         if (perms[j]) {
2014:           PetscMalloc1(size,&ownPerms[j]);
2015:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2016:         }
2017:       }
2018:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2019:     }
2020:     if (rots) {
2021:       PetscScalar **ownRots;

2023:       PetscCalloc1(maxOrient - minOrient,&ownRots);
2024:       for (j = 0; j < maxOrient-minOrient; j++) {
2025:         if (rots[j]) {
2026:           PetscMalloc1(size,&ownRots[j]);
2027:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2028:         }
2029:       }
2030:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2031:     }
2032:   } else {
2033:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2034:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2035:   }
2036:   return(0);
2037: }

2039: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2040: {
2041:   PetscInt              i, j, numStrata;
2042:   PetscSectionSym_Label *sl;
2043:   DMLabel               label;
2044:   PetscErrorCode        ierr;

2047:   sl = (PetscSectionSym_Label *) sym->data;
2048:   numStrata = sl->numStrata;
2049:   label     = sl->label;
2050:   for (i = 0; i < numPoints; i++) {
2051:     PetscInt point = points[2*i];
2052:     PetscInt ornt  = points[2*i+1];

2054:     for (j = 0; j < numStrata; j++) {
2055:       if (label->validIS[j]) {
2056:         PetscInt k;

2058:         ISLocate(label->points[j],point,&k);
2059:         if (k >= 0) break;
2060:       } else {
2061:         PetscBool has;

2063:         PetscHSetIHas(label->ht[j], point, &has);
2064:         if (has) break;
2065:       }
2066:     }
2067:     if ((sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) && (ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1])) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D orientation %D not in range [%D, %D) for stratum %D",point,ornt,sl->minMaxOrients[j][0],sl->minMaxOrients[j][1],j < numStrata ? label->stratumValues[j] : label->defaultValue);
2068:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2069:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2070:   }
2071:   return(0);
2072: }

2074: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2075: {
2076:   PetscSectionSym_Label *sl;
2077:   PetscErrorCode        ierr;

2080:   PetscNewLog(sym,&sl);
2081:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2082:   sym->ops->view      = PetscSectionSymView_Label;
2083:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
2084:   sym->data           = (void *) sl;
2085:   return(0);
2086: }

2088: /*@
2089:   PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label

2091:   Collective

2093:   Input Parameters:
2094: + comm - the MPI communicator for the new symmetry
2095: - label - the label defining the strata

2097:   Output Parameters:
2098: . sym - the section symmetries

2100:   Level: developer

2102: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2103: @*/
2104: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2105: {
2106:   PetscErrorCode        ierr;

2109:   DMInitializePackage();
2110:   PetscSectionSymCreate(comm,sym);
2111:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2112:   PetscSectionSymLabelSetLabel(*sym,label);
2113:   return(0);
2114: }