Actual source code: dmlabel.c

  1: #include <petscdm.h>
  2: #include <petsc/private/dmlabelimpl.h>
  3: #include <petsc/private/sectionimpl.h>
  4: #include <petscsf.h>
  5: #include <petscsection.h>

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

 10:   Collective

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

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

 19:   Level: beginner

 21:   Notes:
 22:   The label name is actually usual PetscObject name.
 23:   One can get/set it with PetscObjectGetName()/PetscObjectSetName().

 25: .seealso: DMLabelDestroy()
 26: @*/
 27: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
 28: {

 33:   DMInitializePackage();

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

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

 52: /*
 53:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 55:   Not collective

 57:   Input parameter:
 58: + label - The DMLabel
 59: - v - The stratum value

 61:   Output parameter:
 62: . label - The DMLabel with stratum in sorted list format

 64:   Level: developer

 66: .seealso: DMLabelCreate()
 67: */
 68: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 69: {
 70:   IS             is;
 71:   PetscInt       off = 0, *pointArray, p;

 74:   if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
 76:   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);
 77:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
 78:   PetscMalloc1(label->stratumSizes[v], &pointArray);
 79:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
 80:   PetscHSetIClear(label->ht[v]);
 81:   PetscSortInt(label->stratumSizes[v], pointArray);
 82:   if (label->bt) {
 83:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 84:       const PetscInt point = pointArray[p];
 85:       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);
 86:       PetscBTSet(label->bt, point - label->pStart);
 87:     }
 88:   }
 89:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
 90:     ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
 91:     PetscFree(pointArray);
 92:   } else {
 93:     ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
 94:   }
 95:   PetscObjectSetName((PetscObject) is, "indices");
 96:   label->points[v]  = is;
 97:   label->validIS[v] = PETSC_TRUE;
 98:   PetscObjectStateIncrease((PetscObject) label);
 99:   return(0);
100: }

102: /*
103:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

105:   Not collective

107:   Input parameter:
108: . label - The DMLabel

110:   Output parameter:
111: . label - The DMLabel with all strata in sorted list format

113:   Level: developer

115: .seealso: DMLabelCreate()
116: */
117: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
118: {
119:   PetscInt       v;

123:   for (v = 0; v < label->numStrata; v++) {
124:     DMLabelMakeValid_Private(label, v);
125:   }
126:   return(0);
127: }

129: /*
130:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

132:   Not collective

134:   Input parameter:
135: + label - The DMLabel
136: - v - The stratum value

138:   Output parameter:
139: . label - The DMLabel with stratum in hash format

141:   Level: developer

143: .seealso: DMLabelCreate()
144: */
145: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
146: {
147:   PetscInt       p;
148:   const PetscInt *points;

151:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
153:   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);
154:   if (label->points[v]) {
155:     ISGetIndices(label->points[v], &points);
156:     for (p = 0; p < label->stratumSizes[v]; ++p) {
157:       PetscHSetIAdd(label->ht[v], points[p]);
158:     }
159:     ISRestoreIndices(label->points[v],&points);
160:     ISDestroy(&(label->points[v]));
161:   }
162:   label->validIS[v] = PETSC_FALSE;
163:   return(0);
164: }

166: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
167: #define DMLABEL_LOOKUP_THRESHOLD 16
168: #endif

170: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
171: {
172:   PetscInt       v;

176:   *index = -1;
177:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
178:     for (v = 0; v < label->numStrata; ++v)
179:       if (label->stratumValues[v] == value) {*index = v; break;}
180:   } else {
181:     PetscHMapIGet(label->hmap, value, index);
182:   }
183:   if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
184:     PetscInt len, loc = -1;
185:     PetscHMapIGetSize(label->hmap, &len);
186:     if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
187:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
188:       PetscHMapIGet(label->hmap, value, &loc);
189:     } else {
190:       for (v = 0; v < label->numStrata; ++v)
191:         if (label->stratumValues[v] == value) {loc = v; break;}
192:     }
193:     if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
194:   }
195:   return(0);
196: }

198: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
199: {
200:   PetscInt       v;
201:   PetscInt      *tmpV;
202:   PetscInt      *tmpS;
203:   PetscHSetI    *tmpH, ht;
204:   IS            *tmpP, is;
205:   PetscBool     *tmpB;
206:   PetscHMapI     hmap = label->hmap;

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

257: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
258: {
261:   DMLabelLookupStratum(label, value, index);
262:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
263:   return(0);
264: }

266: /*@
267:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

269:   Input Parameters:
270: + label - The DMLabel
271: - value - The stratum value

273:   Level: beginner

275: .seealso:  DMLabelCreate(), DMLabelDestroy()
276: @*/
277: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
278: {
279:   PetscInt       v;

284:   DMLabelLookupAddStratum(label, value, &v);
285:   return(0);
286: }

288: /*@
289:   DMLabelAddStrata - Adds new stratum values in a DMLabel

291:   Not collective

293:   Input Parameters:
294: + label - The DMLabel
295: . numStrata - The number of stratum values
296: - stratumValues - The stratum values

298:   Level: beginner

300: .seealso:  DMLabelCreate(), DMLabelDestroy()
301: @*/
302: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
303: {
304:   PetscInt       *values, v;

310:   PetscMalloc1(numStrata, &values);
311:   PetscArraycpy(values, stratumValues, numStrata);
312:   PetscSortRemoveDupsInt(&numStrata, values);
313:   if (!label->numStrata) { /* Fast preallocation */
314:     PetscInt   *tmpV;
315:     PetscInt   *tmpS;
316:     PetscHSetI *tmpH, ht;
317:     IS         *tmpP, is;
318:     PetscBool  *tmpB;
319:     PetscHMapI  hmap = label->hmap;

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

352: /*@
353:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

355:   Not collective

357:   Input Parameters:
358: + label - The DMLabel
359: - valueIS - Index set with stratum values

361:   Level: beginner

363: .seealso:  DMLabelCreate(), DMLabelDestroy()
364: @*/
365: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
366: {
367:   PetscInt       numStrata;
368:   const PetscInt *stratumValues;

374:   ISGetLocalSize(valueIS, &numStrata);
375:   ISGetIndices(valueIS, &stratumValues);
376:   DMLabelAddStrata(label, numStrata, stratumValues);
377:   return(0);
378: }

380: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
381: {
382:   PetscInt       v;
383:   PetscMPIInt    rank;

387:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
388:   PetscViewerASCIIPushSynchronized(viewer);
389:   if (label) {
390:     const char *name;

392:     PetscObjectGetName((PetscObject) label, &name);
393:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
394:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
395:     for (v = 0; v < label->numStrata; ++v) {
396:       const PetscInt value = label->stratumValues[v];
397:       const PetscInt *points;
398:       PetscInt       p;

400:       ISGetIndices(label->points[v], &points);
401:       for (p = 0; p < label->stratumSizes[v]; ++p) {
402:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
403:       }
404:       ISRestoreIndices(label->points[v],&points);
405:     }
406:   }
407:   PetscViewerFlush(viewer);
408:   PetscViewerASCIIPopSynchronized(viewer);
409:   return(0);
410: }

412: /*@C
413:   DMLabelView - View the label

415:   Collective on viewer

417:   Input Parameters:
418: + label - The DMLabel
419: - viewer - The PetscViewer

421:   Level: intermediate

423: .seealso: DMLabelCreate(), DMLabelDestroy()
424: @*/
425: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
426: {
427:   PetscBool      iascii;

432:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
434:   if (label) {DMLabelMakeAllValid_Private(label);}
435:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
436:   if (iascii) {
437:     DMLabelView_Ascii(label, viewer);
438:   }
439:   return(0);
440: }

442: /*@
443:   DMLabelReset - Destroys internal data structures in a DMLabel

445:   Not collective

447:   Input Parameter:
448: . label - The DMLabel

450:   Level: beginner

452: .seealso: DMLabelDestroy(), DMLabelCreate()
453: @*/
454: PetscErrorCode DMLabelReset(DMLabel label)
455: {
456:   PetscInt       v;

461:   for (v = 0; v < label->numStrata; ++v) {
462:     PetscHSetIDestroy(&label->ht[v]);
463:     ISDestroy(&label->points[v]);
464:   }
465:   label->numStrata = 0;
466:   PetscFree(label->stratumValues);
467:   PetscFree(label->stratumSizes);
468:   PetscFree(label->ht);
469:   PetscFree(label->points);
470:   PetscFree(label->validIS);
471:   label->stratumValues = NULL;
472:   label->stratumSizes  = NULL;
473:   label->ht            = NULL;
474:   label->points        = NULL;
475:   label->validIS       = NULL;
476:   PetscHMapIReset(label->hmap);
477:   label->pStart = -1;
478:   label->pEnd   = -1;
479:   PetscBTDestroy(&label->bt);
480:   return(0);
481: }

483: /*@
484:   DMLabelDestroy - Destroys a DMLabel

486:   Collective on label

488:   Input Parameter:
489: . label - The DMLabel

491:   Level: beginner

493: .seealso: DMLabelReset(), DMLabelCreate()
494: @*/
495: PetscErrorCode DMLabelDestroy(DMLabel *label)
496: {

500:   if (!*label) return(0);
502:   if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
503:   DMLabelReset(*label);
504:   PetscHMapIDestroy(&(*label)->hmap);
505:   PetscHeaderDestroy(label);
506:   return(0);
507: }

509: /*@
510:   DMLabelDuplicate - Duplicates a DMLabel

512:   Collective on label

514:   Input Parameter:
515: . label - The DMLabel

517:   Output Parameter:
518: . labelnew - location to put new vector

520:   Level: intermediate

522: .seealso: DMLabelCreate(), DMLabelDestroy()
523: @*/
524: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
525: {
526:   const char    *name;
527:   PetscInt       v;

532:   DMLabelMakeAllValid_Private(label);
533:   PetscObjectGetName((PetscObject) label, &name);
534:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

536:   (*labelnew)->numStrata    = label->numStrata;
537:   (*labelnew)->defaultValue = label->defaultValue;
538:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
539:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
540:   PetscMalloc1(label->numStrata, &(*labelnew)->ht);
541:   PetscMalloc1(label->numStrata, &(*labelnew)->points);
542:   PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
543:   for (v = 0; v < label->numStrata; ++v) {
544:     PetscHSetICreate(&(*labelnew)->ht[v]);
545:     (*labelnew)->stratumValues[v]  = label->stratumValues[v];
546:     (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
547:     PetscObjectReference((PetscObject) (label->points[v]));
548:     (*labelnew)->points[v]         = label->points[v];
549:     (*labelnew)->validIS[v]        = PETSC_TRUE;
550:   }
551:   PetscHMapIDestroy(&(*labelnew)->hmap);
552:   PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
553:   (*labelnew)->pStart = -1;
554:   (*labelnew)->pEnd   = -1;
555:   (*labelnew)->bt     = NULL;
556:   return(0);
557: }

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

562:   Not collective

564:   Input Parameter:
565: . label  - The DMLabel

567:   Level: intermediate

569: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
570: @*/
571: PetscErrorCode DMLabelComputeIndex(DMLabel label)
572: {
573:   PetscInt       pStart = PETSC_MAX_INT, pEnd = -1, v;

578:   DMLabelMakeAllValid_Private(label);
579:   for (v = 0; v < label->numStrata; ++v) {
580:     const PetscInt *points;
581:     PetscInt       i;

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

587:       pStart = PetscMin(point,   pStart);
588:       pEnd   = PetscMax(point+1, pEnd);
589:     }
590:     ISRestoreIndices(label->points[v], &points);
591:   }
592:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
593:   label->pEnd   = pEnd;
594:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
595:   return(0);
596: }

598: /*@
599:   DMLabelCreateIndex - Create an index structure for membership determination

601:   Not collective

603:   Input Parameters:
604: + label  - The DMLabel
605: . pStart - The smallest point
606: - pEnd   - The largest point + 1

608:   Level: intermediate

610: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
611: @*/
612: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
613: {
614:   PetscInt       v;

619:   DMLabelDestroyIndex(label);
620:   DMLabelMakeAllValid_Private(label);
621:   label->pStart = pStart;
622:   label->pEnd   = pEnd;
623:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
624:   PetscBTCreate(pEnd - pStart, &label->bt);
625:   for (v = 0; v < label->numStrata; ++v) {
626:     const PetscInt *points;
627:     PetscInt       i;

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

633:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
634:       PetscBTSet(label->bt, point - pStart);
635:     }
636:     ISRestoreIndices(label->points[v], &points);
637:   }
638:   return(0);
639: }

641: /*@
642:   DMLabelDestroyIndex - Destroy the index structure

644:   Not collective

646:   Input Parameter:
647: . label - the DMLabel

649:   Level: intermediate

651: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
652: @*/
653: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
654: {

659:   label->pStart = -1;
660:   label->pEnd   = -1;
661:   PetscBTDestroy(&label->bt);
662:   return(0);
663: }

665: /*@
666:   DMLabelGetBounds - Return the smallest and largest point in the label

668:   Not collective

670:   Input Parameter:
671: . label - the DMLabel

673:   Output Parameters:
674: + pStart - The smallest point
675: - pEnd   - The largest point + 1

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

679:   Level: intermediate

681: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
682: @*/
683: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
684: {

689:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
690:   if (pStart) {
692:     *pStart = label->pStart;
693:   }
694:   if (pEnd) {
696:     *pEnd = label->pEnd;
697:   }
698:   return(0);
699: }

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

704:   Not collective

706:   Input Parameters:
707: + label - the DMLabel
708: - value - the value

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

713:   Level: developer

715: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
716: @*/
717: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
718: {
719:   PetscInt v;

725:   DMLabelLookupStratum(label, value, &v);
726:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
727:   return(0);
728: }

730: /*@
731:   DMLabelHasPoint - Determine whether a label assigns a value to a point

733:   Not collective

735:   Input Parameters:
736: + label - the DMLabel
737: - point - the point

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

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

744:   Level: developer

746: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
747: @*/
748: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
749: {

755:   DMLabelMakeAllValid_Private(label);
756:   if (PetscDefined(USE_DEBUG)) {
757:     if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
758:     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);
759:   }
760:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
761:   return(0);
762: }

764: /*@
765:   DMLabelStratumHasPoint - Return true if the stratum contains a point

767:   Not collective

769:   Input Parameters:
770: + label - the DMLabel
771: . value - the stratum value
772: - point - the point

774:   Output Parameter:
775: . contains - true if the stratum contains the point

777:   Level: intermediate

779: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
780: @*/
781: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
782: {
783:   PetscInt       v;

789:   *contains = PETSC_FALSE;
790:   DMLabelLookupStratum(label, value, &v);
791:   if (v < 0) return(0);

793:   if (label->validIS[v]) {
794:     PetscInt i;

796:     ISLocate(label->points[v], point, &i);
797:     if (i >= 0) *contains = PETSC_TRUE;
798:   } else {
799:     PetscBool has;

801:     PetscHSetIHas(label->ht[v], point, &has);
802:     if (has) *contains = PETSC_TRUE;
803:   }
804:   return(0);
805: }

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

811:   Not collective

813:   Input parameter:
814: . label - a DMLabel object

816:   Output parameter:
817: . defaultValue - the default value

819:   Level: beginner

821: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
822: @*/
823: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
824: {
827:   *defaultValue = label->defaultValue;
828:   return(0);
829: }

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

835:   Not collective

837:   Input parameter:
838: . label - a DMLabel object

840:   Output parameter:
841: . defaultValue - the default value

843:   Level: beginner

845: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
846: @*/
847: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
848: {
851:   label->defaultValue = defaultValue;
852:   return(0);
853: }

855: /*@
856:   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())

858:   Not collective

860:   Input Parameters:
861: + label - the DMLabel
862: - point - the point

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

867:   Level: intermediate

869: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
870: @*/
871: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
872: {
873:   PetscInt       v;

879:   *value = label->defaultValue;
880:   for (v = 0; v < label->numStrata; ++v) {
881:     if (label->validIS[v]) {
882:       PetscInt i;

884:       ISLocate(label->points[v], point, &i);
885:       if (i >= 0) {
886:         *value = label->stratumValues[v];
887:         break;
888:       }
889:     } else {
890:       PetscBool has;

892:       PetscHSetIHas(label->ht[v], point, &has);
893:       if (has) {
894:         *value = label->stratumValues[v];
895:         break;
896:       }
897:     }
898:   }
899:   return(0);
900: }

902: /*@
903:   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.

905:   Not collective

907:   Input Parameters:
908: + label - the DMLabel
909: . point - the point
910: - value - The point value

912:   Level: intermediate

914: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
915: @*/
916: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
917: {
918:   PetscInt       v;

923:   /* Find label value, add new entry if needed */
924:   if (value == label->defaultValue) return(0);
925:   DMLabelLookupAddStratum(label, value, &v);
926:   /* Set key */
927:   DMLabelMakeInvalid_Private(label, v);
928:   PetscHSetIAdd(label->ht[v], point);
929:   return(0);
930: }

932: /*@
933:   DMLabelClearValue - Clear the value a label assigns to a point

935:   Not collective

937:   Input Parameters:
938: + label - the DMLabel
939: . point - the point
940: - value - The point value

942:   Level: intermediate

944: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
945: @*/
946: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
947: {
948:   PetscInt       v;

953:   /* Find label value */
954:   DMLabelLookupStratum(label, value, &v);
955:   if (v < 0) return(0);

957:   if (label->bt) {
958:     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);
959:     PetscBTClear(label->bt, point - label->pStart);
960:   }

962:   /* Delete key */
963:   DMLabelMakeInvalid_Private(label, v);
964:   PetscHSetIDel(label->ht[v], point);
965:   return(0);
966: }

968: /*@
969:   DMLabelInsertIS - Set all points in the IS to a value

971:   Not collective

973:   Input Parameters:
974: + label - the DMLabel
975: . is    - the point IS
976: - value - The point value

978:   Level: intermediate

980: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
981: @*/
982: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
983: {
984:   PetscInt        v, n, p;
985:   const PetscInt *points;
986:   PetscErrorCode  ierr;

991:   /* Find label value, add new entry if needed */
992:   if (value == label->defaultValue) return(0);
993:   DMLabelLookupAddStratum(label, value, &v);
994:   /* Set keys */
995:   DMLabelMakeInvalid_Private(label, v);
996:   ISGetLocalSize(is, &n);
997:   ISGetIndices(is, &points);
998:   for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
999:   ISRestoreIndices(is, &points);
1000:   return(0);
1001: }

1003: /*@
1004:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

1006:   Not collective

1008:   Input Parameter:
1009: . label - the DMLabel

1011:   Output Parameter:
1012: . numValues - the number of values

1014:   Level: intermediate

1016: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1017: @*/
1018: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1019: {
1023:   *numValues = label->numStrata;
1024:   return(0);
1025: }

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

1030:   Not collective

1032:   Input Parameter:
1033: . label - the DMLabel

1035:   Output Parameter:
1036: . is    - the value IS

1038:   Level: intermediate

1040:   Notes:
1041:   The output IS should be destroyed when no longer needed.

1043: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1044: @*/
1045: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1046: {

1052:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1053:   return(0);
1054: }

1056: /*@
1057:   DMLabelGetValueIndex - Get the index of a given value in the list of values for the DMlabel, or -1 if it is not present

1059:   Not collective

1061:   Input Parameters:
1062: + label - the DMLabel
1063: - value - the value

1065:   Output Parameter:
1066: . index - the index of value in the list of values

1068:   Level: intermediate

1070: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1071: @*/
1072: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1073: {
1074:   PetscInt v;

1079:   /* Do not assume they are sorted */
1080:   for (v = 0; v < label->numStrata; ++v) if (label->stratumValues[v] == value) break;
1081:   if (v >= label->numStrata) *index = -1;
1082:   else                       *index = v;
1083:   return(0);
1084: }

1086: /*@
1087:   DMLabelHasStratum - Determine whether points exist with the given value

1089:   Not collective

1091:   Input Parameters:
1092: + label - the DMLabel
1093: - value - the stratum value

1095:   Output Parameter:
1096: . exists - Flag saying whether points exist

1098:   Level: intermediate

1100: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1101: @*/
1102: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1103: {
1104:   PetscInt       v;

1110:   DMLabelLookupStratum(label, value, &v);
1111:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1112:   return(0);
1113: }

1115: /*@
1116:   DMLabelGetStratumSize - Get the size of a stratum

1118:   Not collective

1120:   Input Parameters:
1121: + label - the DMLabel
1122: - value - the stratum value

1124:   Output Parameter:
1125: . size - The number of points in the stratum

1127:   Level: intermediate

1129: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1130: @*/
1131: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1132: {
1133:   PetscInt       v;

1139:   *size = 0;
1140:   DMLabelLookupStratum(label, value, &v);
1141:   if (v < 0) return(0);
1142:   DMLabelMakeValid_Private(label, v);
1143:   *size = label->stratumSizes[v];
1144:   return(0);
1145: }

1147: /*@
1148:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1150:   Not collective

1152:   Input Parameters:
1153: + label - the DMLabel
1154: - value - the stratum value

1156:   Output Parameters:
1157: + start - the smallest point in the stratum
1158: - end - the largest point in the stratum

1160:   Level: intermediate

1162: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1163: @*/
1164: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1165: {
1166:   PetscInt       v, min, max;

1173:   DMLabelLookupStratum(label, value, &v);
1174:   if (v < 0) return(0);
1175:   DMLabelMakeValid_Private(label, v);
1176:   if (label->stratumSizes[v] <= 0) return(0);
1177:   ISGetMinMax(label->points[v], &min, &max);
1178:   if (start) *start = min;
1179:   if (end)   *end   = max+1;
1180:   return(0);
1181: }

1183: /*@
1184:   DMLabelGetStratumIS - Get an IS with the stratum points

1186:   Not collective

1188:   Input Parameters:
1189: + label - the DMLabel
1190: - value - the stratum value

1192:   Output Parameter:
1193: . points - The stratum points

1195:   Level: intermediate

1197:   Notes:
1198:   The output IS should be destroyed when no longer needed.
1199:   Returns NULL if the stratum is empty.

1201: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1202: @*/
1203: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1204: {
1205:   PetscInt       v;

1211:   *points = NULL;
1212:   DMLabelLookupStratum(label, value, &v);
1213:   if (v < 0) return(0);
1214:   DMLabelMakeValid_Private(label, v);
1215:   PetscObjectReference((PetscObject) label->points[v]);
1216:   *points = label->points[v];
1217:   return(0);
1218: }

1220: /*@
1221:   DMLabelSetStratumIS - Set the stratum points using an IS

1223:   Not collective

1225:   Input Parameters:
1226: + label - the DMLabel
1227: . value - the stratum value
1228: - points - The stratum points

1230:   Level: intermediate

1232: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1233: @*/
1234: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1235: {
1236:   PetscInt       v;

1242:   DMLabelLookupAddStratum(label, value, &v);
1243:   if (is == label->points[v]) return(0);
1244:   DMLabelClearStratum(label, value);
1245:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1246:   PetscObjectReference((PetscObject)is);
1247:   ISDestroy(&(label->points[v]));
1248:   label->points[v]  = is;
1249:   label->validIS[v] = PETSC_TRUE;
1250:   PetscObjectStateIncrease((PetscObject) label);
1251:   if (label->bt) {
1252:     const PetscInt *points;
1253:     PetscInt p;

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

1259:       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);
1260:       PetscBTSet(label->bt, point - label->pStart);
1261:     }
1262:   }
1263:   return(0);
1264: }

1266: /*@
1267:   DMLabelClearStratum - Remove a stratum

1269:   Not collective

1271:   Input Parameters:
1272: + label - the DMLabel
1273: - value - the stratum value

1275:   Level: intermediate

1277: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1278: @*/
1279: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1280: {
1281:   PetscInt       v;

1286:   DMLabelLookupStratum(label, value, &v);
1287:   if (v < 0) return(0);
1288:   if (label->validIS[v]) {
1289:     if (label->bt) {
1290:       PetscInt       i;
1291:       const PetscInt *points;

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

1297:         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);
1298:         PetscBTClear(label->bt, point - label->pStart);
1299:       }
1300:       ISRestoreIndices(label->points[v], &points);
1301:     }
1302:     label->stratumSizes[v] = 0;
1303:     ISDestroy(&label->points[v]);
1304:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1305:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1306:     PetscObjectStateIncrease((PetscObject) label);
1307:   } else {
1308:     PetscHSetIClear(label->ht[v]);
1309:   }
1310:   return(0);
1311: }

1313: /*@
1314:   DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value

1316:   Not collective

1318:   Input Parameters:
1319: + label  - The DMLabel
1320: . value  - The label value for all points
1321: . pStart - The first point
1322: - pEnd   - A point beyond all marked points

1324:   Note: The marks points are [pStart, pEnd), and only the bounds are stored.

1326:   Level: intermediate

1328: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1329: @*/
1330: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1331: {
1332:   IS             pIS;

1336:   ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1337:   DMLabelSetStratumIS(label, value, pIS);
1338:   ISDestroy(&pIS);
1339:   return(0);
1340: }

1342: /*@
1343:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

1345:   Not collective

1347:   Input Parameters:
1348: + label  - The DMLabel
1349: . value  - The label value
1350: - p      - A point with this value

1352:   Output Parameter:
1353: . index  - The index of this point in the stratum, or -1 if the point is not in the stratum or the stratum does not exist

1355:   Level: intermediate

1357: .seealso: DMLabelGetValueIndex(), DMLabelGetStratumIS(), DMLabelCreate()
1358: @*/
1359: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1360: {
1361:   const PetscInt *indices;
1362:   PetscInt        v;
1363:   PetscErrorCode  ierr;

1368:   *index = -1;
1369:   DMLabelLookupStratum(label, value, &v);
1370:   if (v < 0) return(0);
1371:   DMLabelMakeValid_Private(label, v);
1372:   ISGetIndices(label->points[v], &indices);
1373:   PetscFindInt(p, label->stratumSizes[v], indices, index);
1374:   ISRestoreIndices(label->points[v], &indices);
1375:   return(0);
1376: }

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

1381:   Not collective

1383:   Input Parameters:
1384: + label - the DMLabel
1385: . start - the first point kept
1386: - end - one more than the last point kept

1388:   Level: intermediate

1390: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1391: @*/
1392: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1393: {
1394:   PetscInt       v;

1399:   DMLabelDestroyIndex(label);
1400:   DMLabelMakeAllValid_Private(label);
1401:   for (v = 0; v < label->numStrata; ++v) {
1402:     ISGeneralFilter(label->points[v], start, end);
1403:   }
1404:   DMLabelCreateIndex(label, start, end);
1405:   return(0);
1406: }

1408: /*@
1409:   DMLabelPermute - Create a new label with permuted points

1411:   Not collective

1413:   Input Parameters:
1414: + label - the DMLabel
1415: - permutation - the point permutation

1417:   Output Parameter:
1418: . labelnew - the new label containing the permuted points

1420:   Level: intermediate

1422: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1423: @*/
1424: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1425: {
1426:   const PetscInt *perm;
1427:   PetscInt        numValues, numPoints, v, q;
1428:   PetscErrorCode  ierr;

1433:   DMLabelMakeAllValid_Private(label);
1434:   DMLabelDuplicate(label, labelNew);
1435:   DMLabelGetNumValues(*labelNew, &numValues);
1436:   ISGetLocalSize(permutation, &numPoints);
1437:   ISGetIndices(permutation, &perm);
1438:   for (v = 0; v < numValues; ++v) {
1439:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1440:     const PetscInt *points;
1441:     PetscInt *pointsNew;

1443:     ISGetIndices((*labelNew)->points[v],&points);
1444:     PetscMalloc1(size,&pointsNew);
1445:     for (q = 0; q < size; ++q) {
1446:       const PetscInt point = points[q];

1448:       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);
1449:       pointsNew[q] = perm[point];
1450:     }
1451:     ISRestoreIndices((*labelNew)->points[v],&points);
1452:     PetscSortInt(size, pointsNew);
1453:     ISDestroy(&((*labelNew)->points[v]));
1454:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1455:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1456:       PetscFree(pointsNew);
1457:     } else {
1458:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1459:     }
1460:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1461:   }
1462:   ISRestoreIndices(permutation, &perm);
1463:   if (label->bt) {
1464:     PetscBTDestroy(&label->bt);
1465:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1466:   }
1467:   return(0);
1468: }

1470: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1471: {
1472:   MPI_Comm       comm;
1473:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1474:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1475:   PetscSection   rootSection;
1476:   PetscSF        labelSF;

1480:   if (label) {DMLabelMakeAllValid_Private(label);}
1481:   PetscObjectGetComm((PetscObject)sf, &comm);
1482:   /* Build a section of stratum values per point, generate the according SF
1483:      and distribute point-wise stratum values to leaves. */
1484:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1485:   PetscSectionCreate(comm, &rootSection);
1486:   PetscSectionSetChart(rootSection, 0, nroots);
1487:   if (label) {
1488:     for (s = 0; s < label->numStrata; ++s) {
1489:       const PetscInt *points;

1491:       ISGetIndices(label->points[s], &points);
1492:       for (l = 0; l < label->stratumSizes[s]; l++) {
1493:         PetscSectionGetDof(rootSection, points[l], &dof);
1494:         PetscSectionSetDof(rootSection, points[l], dof+1);
1495:       }
1496:       ISRestoreIndices(label->points[s], &points);
1497:     }
1498:   }
1499:   PetscSectionSetUp(rootSection);
1500:   /* Create a point-wise array of stratum values */
1501:   PetscSectionGetStorageSize(rootSection, &size);
1502:   PetscMalloc1(size, &rootStrata);
1503:   PetscCalloc1(nroots, &rootIdx);
1504:   if (label) {
1505:     for (s = 0; s < label->numStrata; ++s) {
1506:       const PetscInt *points;

1508:       ISGetIndices(label->points[s], &points);
1509:       for (l = 0; l < label->stratumSizes[s]; l++) {
1510:         const PetscInt p = points[l];
1511:         PetscSectionGetOffset(rootSection, p, &offset);
1512:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1513:       }
1514:       ISRestoreIndices(label->points[s], &points);
1515:     }
1516:   }
1517:   /* Build SF that maps label points to remote processes */
1518:   PetscSectionCreate(comm, leafSection);
1519:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1520:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1521:   PetscFree(remoteOffsets);
1522:   /* Send the strata for each point over the derived SF */
1523:   PetscSectionGetStorageSize(*leafSection, &size);
1524:   PetscMalloc1(size, leafStrata);
1525:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1526:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1527:   /* Clean up */
1528:   PetscFree(rootStrata);
1529:   PetscFree(rootIdx);
1530:   PetscSectionDestroy(&rootSection);
1531:   PetscSFDestroy(&labelSF);
1532:   return(0);
1533: }

1535: /*@
1536:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1538:   Collective on sf

1540:   Input Parameters:
1541: + label - the DMLabel
1542: - sf    - the map from old to new distribution

1544:   Output Parameter:
1545: . labelnew - the new redistributed label

1547:   Level: intermediate

1549: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1550: @*/
1551: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1552: {
1553:   MPI_Comm       comm;
1554:   PetscSection   leafSection;
1555:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1556:   PetscInt      *leafStrata, *strataIdx;
1557:   PetscInt     **points;
1558:   const char    *lname = NULL;
1559:   char          *name;
1560:   PetscInt       nameSize;
1561:   PetscHSetI     stratumHash;
1562:   size_t         len = 0;
1563:   PetscMPIInt    rank;

1568:   if (label) {
1570:     DMLabelMakeAllValid_Private(label);
1571:   }
1572:   PetscObjectGetComm((PetscObject)sf, &comm);
1573:   MPI_Comm_rank(comm, &rank);
1574:   /* Bcast name */
1575:   if (rank == 0) {
1576:     PetscObjectGetName((PetscObject) label, &lname);
1577:     PetscStrlen(lname, &len);
1578:   }
1579:   nameSize = len;
1580:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1581:   PetscMalloc1(nameSize+1, &name);
1582:   if (rank == 0) {PetscArraycpy(name, lname, nameSize+1);}
1583:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1584:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1585:   PetscFree(name);
1586:   /* Bcast defaultValue */
1587:   if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1588:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1589:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1590:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1591:   /* Determine received stratum values and initialise new label*/
1592:   PetscHSetICreate(&stratumHash);
1593:   PetscSectionGetStorageSize(leafSection, &size);
1594:   for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1595:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1596:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1597:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1598:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1599:   /* Turn leafStrata into indices rather than stratum values */
1600:   offset = 0;
1601:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1602:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1603:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1604:     PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1605:   }
1606:   for (p = 0; p < size; ++p) {
1607:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1608:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1609:     }
1610:   }
1611:   /* Rebuild the point strata on the receiver */
1612:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1613:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1614:   for (p=pStart; p<pEnd; p++) {
1615:     PetscSectionGetDof(leafSection, p, &dof);
1616:     PetscSectionGetOffset(leafSection, p, &offset);
1617:     for (s=0; s<dof; s++) {
1618:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1619:     }
1620:   }
1621:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1622:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1623:   PetscMalloc1((*labelNew)->numStrata,&points);
1624:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1625:     PetscHSetICreate(&(*labelNew)->ht[s]);
1626:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1627:   }
1628:   /* Insert points into new strata */
1629:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1630:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1631:   for (p=pStart; p<pEnd; p++) {
1632:     PetscSectionGetDof(leafSection, p, &dof);
1633:     PetscSectionGetOffset(leafSection, p, &offset);
1634:     for (s=0; s<dof; s++) {
1635:       stratum = leafStrata[offset+s];
1636:       points[stratum][strataIdx[stratum]++] = p;
1637:     }
1638:   }
1639:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1640:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1641:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1642:   }
1643:   PetscFree(points);
1644:   PetscHSetIDestroy(&stratumHash);
1645:   PetscFree(leafStrata);
1646:   PetscFree(strataIdx);
1647:   PetscSectionDestroy(&leafSection);
1648:   return(0);
1649: }

1651: /*@
1652:   DMLabelGather - Gather all label values from leafs into roots

1654:   Collective on sf

1656:   Input Parameters:
1657: + label - the DMLabel
1658: - sf - the Star Forest point communication map

1660:   Output Parameters:
1661: . labelNew - the new DMLabel with localised leaf values

1663:   Level: developer

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

1667: .seealso: DMLabelDistribute()
1668: @*/
1669: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1670: {
1671:   MPI_Comm       comm;
1672:   PetscSection   rootSection;
1673:   PetscSF        sfLabel;
1674:   PetscSFNode   *rootPoints, *leafPoints;
1675:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1676:   const PetscInt *rootDegree, *ilocal;
1677:   PetscInt       *rootStrata;
1678:   const char    *lname;
1679:   char          *name;
1680:   PetscInt       nameSize;
1681:   size_t         len = 0;
1682:   PetscMPIInt    rank, size;

1688:   PetscObjectGetComm((PetscObject)sf, &comm);
1689:   MPI_Comm_rank(comm, &rank);
1690:   MPI_Comm_size(comm, &size);
1691:   /* Bcast name */
1692:   if (rank == 0) {
1693:     PetscObjectGetName((PetscObject) label, &lname);
1694:     PetscStrlen(lname, &len);
1695:   }
1696:   nameSize = len;
1697:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1698:   PetscMalloc1(nameSize+1, &name);
1699:   if (rank == 0) {PetscArraycpy(name, lname, nameSize+1);}
1700:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1701:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1702:   PetscFree(name);
1703:   /* Gather rank/index pairs of leaves into local roots to build
1704:      an inverse, multi-rooted SF. Note that this ignores local leaf
1705:      indexing due to the use of the multiSF in PetscSFGather. */
1706:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1707:   PetscMalloc1(nroots, &leafPoints);
1708:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1709:   for (p = 0; p < nleaves; p++) {
1710:     PetscInt ilp = ilocal ? ilocal[p] : p;

1712:     leafPoints[ilp].index = ilp;
1713:     leafPoints[ilp].rank  = rank;
1714:   }
1715:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1716:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1717:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1718:   PetscMalloc1(nmultiroots, &rootPoints);
1719:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1720:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1721:   PetscSFCreate(comm,& sfLabel);
1722:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1723:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1724:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1725:   /* Rebuild the point strata on the receiver */
1726:   for (p = 0, idx = 0; p < nroots; p++) {
1727:     for (d = 0; d < rootDegree[p]; d++) {
1728:       PetscSectionGetDof(rootSection, idx+d, &dof);
1729:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1730:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1731:     }
1732:     idx += rootDegree[p];
1733:   }
1734:   PetscFree(leafPoints);
1735:   PetscFree(rootStrata);
1736:   PetscSectionDestroy(&rootSection);
1737:   PetscSFDestroy(&sfLabel);
1738:   return(0);
1739: }

1741: /*@
1742:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1744:   Not collective

1746:   Input Parameter:
1747: . label - the DMLabel

1749:   Output Parameters:
1750: + section - the section giving offsets for each stratum
1751: - is - An IS containing all the label points

1753:   Level: developer

1755: .seealso: DMLabelDistribute()
1756: @*/
1757: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1758: {
1759:   IS              vIS;
1760:   const PetscInt *values;
1761:   PetscInt       *points;
1762:   PetscInt        nV, vS = 0, vE = 0, v, N;
1763:   PetscErrorCode  ierr;

1767:   DMLabelGetNumValues(label, &nV);
1768:   DMLabelGetValueIS(label, &vIS);
1769:   ISGetIndices(vIS, &values);
1770:   if (nV) {vS = values[0]; vE = values[0]+1;}
1771:   for (v = 1; v < nV; ++v) {
1772:     vS = PetscMin(vS, values[v]);
1773:     vE = PetscMax(vE, values[v]+1);
1774:   }
1775:   PetscSectionCreate(PETSC_COMM_SELF, section);
1776:   PetscSectionSetChart(*section, vS, vE);
1777:   for (v = 0; v < nV; ++v) {
1778:     PetscInt n;

1780:     DMLabelGetStratumSize(label, values[v], &n);
1781:     PetscSectionSetDof(*section, values[v], n);
1782:   }
1783:   PetscSectionSetUp(*section);
1784:   PetscSectionGetStorageSize(*section, &N);
1785:   PetscMalloc1(N, &points);
1786:   for (v = 0; v < nV; ++v) {
1787:     IS              is;
1788:     const PetscInt *spoints;
1789:     PetscInt        dof, off, p;

1791:     PetscSectionGetDof(*section, values[v], &dof);
1792:     PetscSectionGetOffset(*section, values[v], &off);
1793:     DMLabelGetStratumIS(label, values[v], &is);
1794:     ISGetIndices(is, &spoints);
1795:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1796:     ISRestoreIndices(is, &spoints);
1797:     ISDestroy(&is);
1798:   }
1799:   ISRestoreIndices(vIS, &values);
1800:   ISDestroy(&vIS);
1801:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1802:   return(0);
1803: }

1805: /*@
1806:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1807:   the local section and an SF describing the section point overlap.

1809:   Collective on sf

1811:   Input Parameters:
1812:   + s - The PetscSection for the local field layout
1813:   . sf - The SF describing parallel layout of the section points
1814:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1815:   . label - The label specifying the points
1816:   - labelValue - The label stratum specifying the points

1818:   Output Parameter:
1819:   . gsection - The PetscSection for the global field layout

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

1823:   Level: developer

1825: .seealso: PetscSectionCreate()
1826: @*/
1827: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1828: {
1829:   PetscInt      *neg = NULL, *tmpOff = NULL;
1830:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1837:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1838:   PetscSectionGetChart(s, &pStart, &pEnd);
1839:   PetscSectionSetChart(*gsection, pStart, pEnd);
1840:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1841:   if (nroots >= 0) {
1842:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1843:     PetscCalloc1(nroots, &neg);
1844:     if (nroots > pEnd-pStart) {
1845:       PetscCalloc1(nroots, &tmpOff);
1846:     } else {
1847:       tmpOff = &(*gsection)->atlasDof[-pStart];
1848:     }
1849:   }
1850:   /* Mark ghost points with negative dof */
1851:   for (p = pStart; p < pEnd; ++p) {
1852:     PetscInt value;

1854:     DMLabelGetValue(label, p, &value);
1855:     if (value != labelValue) continue;
1856:     PetscSectionGetDof(s, p, &dof);
1857:     PetscSectionSetDof(*gsection, p, dof);
1858:     PetscSectionGetConstraintDof(s, p, &cdof);
1859:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1860:     if (neg) neg[p] = -(dof+1);
1861:   }
1862:   PetscSectionSetUpBC(*gsection);
1863:   if (nroots >= 0) {
1864:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1865:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1866:     if (nroots > pEnd-pStart) {
1867:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1868:     }
1869:   }
1870:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1871:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1872:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1873:     (*gsection)->atlasOff[p] = off;
1874:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1875:   }
1876:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1877:   globalOff -= off;
1878:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1879:     (*gsection)->atlasOff[p] += globalOff;
1880:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1881:   }
1882:   /* Put in negative offsets for ghost points */
1883:   if (nroots >= 0) {
1884:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1885:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1886:     if (nroots > pEnd-pStart) {
1887:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1888:     }
1889:   }
1890:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1891:   PetscFree(neg);
1892:   return(0);
1893: }

1895: typedef struct _n_PetscSectionSym_Label
1896: {
1897:   DMLabel           label;
1898:   PetscCopyMode     *modes;
1899:   PetscInt          *sizes;
1900:   const PetscInt    ***perms;
1901:   const PetscScalar ***rots;
1902:   PetscInt          (*minMaxOrients)[2];
1903:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1904: } PetscSectionSym_Label;

1906: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1907: {
1908:   PetscInt              i, j;
1909:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1910:   PetscErrorCode        ierr;

1913:   for (i = 0; i <= sl->numStrata; i++) {
1914:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1915:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1916:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1917:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1918:       }
1919:       if (sl->perms[i]) {
1920:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1922:         PetscFree(perms);
1923:       }
1924:       if (sl->rots[i]) {
1925:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1927:         PetscFree(rots);
1928:       }
1929:     }
1930:   }
1931:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1932:   DMLabelDestroy(&sl->label);
1933:   sl->numStrata = 0;
1934:   return(0);
1935: }

1937: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1938: {

1942:   PetscSectionSymLabelReset(sym);
1943:   PetscFree(sym->data);
1944:   return(0);
1945: }

1947: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1948: {
1949:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1950:   PetscBool             isAscii;
1951:   DMLabel               label = sl->label;
1952:   const char           *name;
1953:   PetscErrorCode        ierr;

1956:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1957:   if (isAscii) {
1958:     PetscInt          i, j, k;
1959:     PetscViewerFormat format;

1961:     PetscViewerGetFormat(viewer,&format);
1962:     if (label) {
1963:       PetscViewerGetFormat(viewer,&format);
1964:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1965:         PetscViewerASCIIPushTab(viewer);
1966:         DMLabelView(label, viewer);
1967:         PetscViewerASCIIPopTab(viewer);
1968:       } else {
1969:         PetscObjectGetName((PetscObject) sl->label, &name);
1970:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
1971:       }
1972:     } else {
1973:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1974:     }
1975:     PetscViewerASCIIPushTab(viewer);
1976:     for (i = 0; i <= sl->numStrata; i++) {
1977:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1979:       if (!(sl->perms[i] || sl->rots[i])) {
1980:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1981:       } else {
1982:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1983:         PetscViewerASCIIPushTab(viewer);
1984:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1985:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1986:           PetscViewerASCIIPushTab(viewer);
1987:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1988:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1989:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1990:             } else {
1991:               PetscInt tab;

1993:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1994:               PetscViewerASCIIPushTab(viewer);
1995:               PetscViewerASCIIGetTab(viewer,&tab);
1996:               if (sl->perms[i] && sl->perms[i][j]) {
1997:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1998:                 PetscViewerASCIISetTab(viewer,0);
1999:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
2000:                 PetscViewerASCIIPrintf(viewer,"\n");
2001:                 PetscViewerASCIISetTab(viewer,tab);
2002:               }
2003:               if (sl->rots[i] && sl->rots[i][j]) {
2004:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
2005:                 PetscViewerASCIISetTab(viewer,0);
2006: #if defined(PETSC_USE_COMPLEX)
2007:                 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]));}
2008: #else
2009:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
2010: #endif
2011:                 PetscViewerASCIIPrintf(viewer,"\n");
2012:                 PetscViewerASCIISetTab(viewer,tab);
2013:               }
2014:               PetscViewerASCIIPopTab(viewer);
2015:             }
2016:           }
2017:           PetscViewerASCIIPopTab(viewer);
2018:         }
2019:         PetscViewerASCIIPopTab(viewer);
2020:       }
2021:     }
2022:     PetscViewerASCIIPopTab(viewer);
2023:   }
2024:   return(0);
2025: }

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

2030:   Logically collective on sym

2032:   Input parameters:
2033: + sym - the section symmetries
2034: - label - the DMLabel describing the types of points

2036:   Level: developer:

2038: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
2039: @*/
2040: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2041: {
2042:   PetscSectionSym_Label *sl;
2043:   PetscErrorCode        ierr;

2047:   sl = (PetscSectionSym_Label *) sym->data;
2048:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
2049:   if (label) {
2050:     sl->label = label;
2051:     PetscObjectReference((PetscObject) label);
2052:     DMLabelGetNumValues(label,&sl->numStrata);
2053:     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);
2054:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
2055:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
2056:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
2057:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
2058:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
2059:   }
2060:   return(0);
2061: }

2063: /*@C
2064:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

2066:   Logically collective on sym

2068:   InputParameters:
2069: + sym       - the section symmetries
2070: . stratum   - the stratum value in the label that we are assigning symmetries for
2071: . size      - the number of dofs for points in the stratum of the label
2072: . minOrient - the smallest orientation for a point in this stratum
2073: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2074: . mode      - how sym should copy the perms and rots arrays
2075: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2076: - 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

2078:   Level: developer

2080: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2081: @*/
2082: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2083: {
2084:   PetscSectionSym_Label *sl;
2085:   const char            *name;
2086:   PetscInt               i, j, k;
2087:   PetscErrorCode         ierr;

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

2096:     if (stratum == value) break;
2097:   }
2098:   PetscObjectGetName((PetscObject) sl->label, &name);
2099:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2100:   sl->sizes[i] = size;
2101:   sl->modes[i] = mode;
2102:   sl->minMaxOrients[i][0] = minOrient;
2103:   sl->minMaxOrients[i][1] = maxOrient;
2104:   if (mode == PETSC_COPY_VALUES) {
2105:     if (perms) {
2106:       PetscInt    **ownPerms;

2108:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
2109:       for (j = 0; j < maxOrient-minOrient; j++) {
2110:         if (perms[j]) {
2111:           PetscMalloc1(size,&ownPerms[j]);
2112:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2113:         }
2114:       }
2115:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2116:     }
2117:     if (rots) {
2118:       PetscScalar **ownRots;

2120:       PetscCalloc1(maxOrient - minOrient,&ownRots);
2121:       for (j = 0; j < maxOrient-minOrient; j++) {
2122:         if (rots[j]) {
2123:           PetscMalloc1(size,&ownRots[j]);
2124:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2125:         }
2126:       }
2127:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2128:     }
2129:   } else {
2130:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2131:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2132:   }
2133:   return(0);
2134: }

2136: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2137: {
2138:   PetscInt              i, j, numStrata;
2139:   PetscSectionSym_Label *sl;
2140:   DMLabel               label;
2141:   PetscErrorCode        ierr;

2144:   sl = (PetscSectionSym_Label *) sym->data;
2145:   numStrata = sl->numStrata;
2146:   label     = sl->label;
2147:   for (i = 0; i < numPoints; i++) {
2148:     PetscInt point = points[2*i];
2149:     PetscInt ornt  = points[2*i+1];

2151:     for (j = 0; j < numStrata; j++) {
2152:       if (label->validIS[j]) {
2153:         PetscInt k;

2155:         ISLocate(label->points[j],point,&k);
2156:         if (k >= 0) break;
2157:       } else {
2158:         PetscBool has;

2160:         PetscHSetIHas(label->ht[j], point, &has);
2161:         if (has) break;
2162:       }
2163:     }
2164:     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);
2165:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2166:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2167:   }
2168:   return(0);
2169: }

2171: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2172: {
2173:   PetscSectionSym_Label *sl;
2174:   PetscErrorCode        ierr;

2177:   PetscNewLog(sym,&sl);
2178:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2179:   sym->ops->view      = PetscSectionSymView_Label;
2180:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
2181:   sym->data           = (void *) sl;
2182:   return(0);
2183: }

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

2188:   Collective

2190:   Input Parameters:
2191: + comm - the MPI communicator for the new symmetry
2192: - label - the label defining the strata

2194:   Output Parameters:
2195: . sym - the section symmetries

2197:   Level: developer

2199: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2200: @*/
2201: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2202: {
2203:   PetscErrorCode        ierr;

2206:   DMInitializePackage();
2207:   PetscSectionSymCreate(comm,sym);
2208:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2209:   PetscSectionSymLabelSetLabel(*sym,label);
2210:   return(0);
2211: }