Actual source code: dmlabel.c

petsc-master 2019-12-13
Report Typos and Errors
  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:   ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
 90:   PetscObjectSetName((PetscObject) is, "indices");
 91:   label->points[v]  = is;
 92:   label->validIS[v] = PETSC_TRUE;
 93:   PetscObjectStateIncrease((PetscObject) label);
 94:   return(0);
 95: }

 97: /*
 98:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

100:   Not collective

102:   Input parameter:
103: . label - The DMLabel

105:   Output parameter:
106: . label - The DMLabel with all strata in sorted list format

108:   Level: developer

110: .seealso: DMLabelCreate()
111: */
112: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
113: {
114:   PetscInt       v;

118:   for (v = 0; v < label->numStrata; v++){
119:     DMLabelMakeValid_Private(label, v);
120:   }
121:   return(0);
122: }

124: /*
125:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

127:   Not collective

129:   Input parameter:
130: + label - The DMLabel
131: - v - The stratum value

133:   Output parameter:
134: . label - The DMLabel with stratum in hash format

136:   Level: developer

138: .seealso: DMLabelCreate()
139: */
140: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
141: {
142:   PetscInt       p;
143:   const PetscInt *points;

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

161: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
162: #define DMLABEL_LOOKUP_THRESHOLD 16
163: #endif

165: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
166: {
167:   PetscInt       v;

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

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

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

254: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
255: {
258:   DMLabelLookupStratum(label, value, index);
259:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
260:   return(0);
261: }

263: /*@
264:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

266:   Input Parameter:
267: + label - The DMLabel
268: - value - The stratum value

270:   Level: beginner

272: .seealso:  DMLabelCreate(), DMLabelDestroy()
273: @*/
274: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
275: {
276:   PetscInt       v;

281:   DMLabelLookupAddStratum(label, value, &v);
282:   return(0);
283: }

285: /*@
286:   DMLabelAddStrata - Adds new stratum values in a DMLabel

288:   Not collective

290:   Input Parameter:
291: + label - The DMLabel
292: . numStrata - The number of stratum values
293: - stratumValues - The stratum values

295:   Level: beginner

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

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

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

349: /*@
350:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

352:   Not collective

354:   Input Parameter:
355: + label - The DMLabel
356: - valueIS - Index set with stratum values

358:   Level: beginner

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

371:   ISGetLocalSize(valueIS, &numStrata);
372:   ISGetIndices(valueIS, &stratumValues);
373:   DMLabelAddStrata(label, numStrata, stratumValues);
374:   return(0);
375: }

377: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
378: {
379:   PetscInt       v;
380:   PetscMPIInt    rank;

384:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
385:   PetscViewerASCIIPushSynchronized(viewer);
386:   if (label) {
387:     const char *name;

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

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

409: /*@C
410:   DMLabelView - View the label

412:   Collective on viewer

414:   Input Parameters:
415: + label - The DMLabel
416: - viewer - The PetscViewer

418:   Level: intermediate

420: .seealso: DMLabelCreate(), DMLabelDestroy()
421: @*/
422: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
423: {
424:   PetscBool      iascii;

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

439: /*@
440:   DMLabelReset - Destroys internal data structures in a DMLabel

442:   Not collective

444:   Input Parameter:
445: . label - The DMLabel

447:   Level: beginner

449: .seealso: DMLabelDestroy(), DMLabelCreate()
450: @*/
451: PetscErrorCode DMLabelReset(DMLabel label)
452: {
453:   PetscInt       v;

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

475: /*@
476:   DMLabelDestroy - Destroys a DMLabel

478:   Collective on label

480:   Input Parameter:
481: . label - The DMLabel

483:   Level: beginner

485: .seealso: DMLabelReset(), DMLabelCreate()
486: @*/
487: PetscErrorCode DMLabelDestroy(DMLabel *label)
488: {

492:   if (!*label) return(0);
494:   if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
495:   DMLabelReset(*label);
496:   PetscHMapIDestroy(&(*label)->hmap);
497:   PetscHeaderDestroy(label);
498:   return(0);
499: }

501: /*@
502:   DMLabelDuplicate - Duplicates a DMLabel

504:   Collective on label

506:   Input Parameter:
507: . label - The DMLabel

509:   Output Parameter:
510: . labelnew - location to put new vector

512:   Level: intermediate

514: .seealso: DMLabelCreate(), DMLabelDestroy()
515: @*/
516: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
517: {
518:   const char    *name;
519:   PetscInt       v;

524:   DMLabelMakeAllValid_Private(label);
525:   PetscObjectGetName((PetscObject) label, &name);
526:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

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

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

554:   Not collective

556:   Input Parameter:
557: . label  - The DMLabel

559:   Level: intermediate

561: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
562: @*/
563: PetscErrorCode DMLabelComputeIndex(DMLabel label)
564: {
565:   PetscInt       pStart = PETSC_MAX_INT, pEnd = -1, v;

570:   DMLabelMakeAllValid_Private(label);
571:   for (v = 0; v < label->numStrata; ++v) {
572:     const PetscInt *points;
573:     PetscInt       i;

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

579:       pStart = PetscMin(point,   pStart);
580:       pEnd   = PetscMax(point+1, pEnd);
581:     }
582:     ISRestoreIndices(label->points[v], &points);
583:   }
584:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
585:   label->pEnd   = pEnd;
586:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
587:   return(0);
588: }

590: /*@
591:   DMLabelCreateIndex - Create an index structure for membership determination

593:   Not collective

595:   Input Parameters:
596: + label  - The DMLabel
597: . pStart - The smallest point
598: - pEnd   - The largest point + 1

600:   Level: intermediate

602: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
603: @*/
604: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
605: {
606:   PetscInt       v;

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

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

625:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
626:       PetscBTSet(label->bt, point - pStart);
627:     }
628:     ISRestoreIndices(label->points[v], &points);
629:   }
630:   return(0);
631: }

633: /*@
634:   DMLabelDestroyIndex - Destroy the index structure

636:   Not collective

638:   Input Parameter:
639: . label - the DMLabel

641:   Level: intermediate

643: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
644: @*/
645: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
646: {

651:   label->pStart = -1;
652:   label->pEnd   = -1;
653:   PetscBTDestroy(&label->bt);
654:   return(0);
655: }

657: /*@
658:   DMLabelGetBounds - Return the smallest and largest point in the label

660:   Not collective

662:   Input Parameter:
663: . label - the DMLabel

665:   Output Parameters:
666: + pStart - The smallest point
667: - pEnd   - The largest point + 1

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

671:   Level: intermediate

673: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
674: @*/
675: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
676: {

681:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
682:   if (pStart) {
684:     *pStart = label->pStart;
685:   }
686:   if (pEnd) {
688:     *pEnd = label->pEnd;
689:   }
690:   return(0);
691: }

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

696:   Not collective

698:   Input Parameters:
699: + label - the DMLabel
700: - value - the value

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

705:   Level: developer

707: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
708: @*/
709: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
710: {
711:   PetscInt v;

717:   DMLabelLookupStratum(label, value, &v);
718:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
719:   return(0);
720: }

722: /*@
723:   DMLabelHasPoint - Determine whether a label assigns a value to a point

725:   Not collective

727:   Input Parameters:
728: + label - the DMLabel
729: - point - the point

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

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

736:   Level: developer

738: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
739: @*/
740: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
741: {

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

756: /*@
757:   DMLabelStratumHasPoint - Return true if the stratum contains a point

759:   Not collective

761:   Input Parameters:
762: + label - the DMLabel
763: . value - the stratum value
764: - point - the point

766:   Output Parameter:
767: . contains - true if the stratum contains the point

769:   Level: intermediate

771: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
772: @*/
773: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
774: {
775:   PetscInt       v;

781:   *contains = PETSC_FALSE;
782:   DMLabelLookupStratum(label, value, &v);
783:   if (v < 0) return(0);

785:   if (label->validIS[v]) {
786:     PetscInt i;

788:     ISLocate(label->points[v], point, &i);
789:     if (i >= 0) *contains = PETSC_TRUE;
790:   } else {
791:     PetscBool has;

793:     PetscHSetIHas(label->ht[v], point, &has);
794:     if (has) *contains = PETSC_TRUE;
795:   }
796:   return(0);
797: }

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

803:   Not collective

805:   Input parameter:
806: . label - a DMLabel object

808:   Output parameter:
809: . defaultValue - the default value

811:   Level: beginner

813: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
814: @*/
815: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
816: {
819:   *defaultValue = label->defaultValue;
820:   return(0);
821: }

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

827:   Not collective

829:   Input parameter:
830: . label - a DMLabel object

832:   Output parameter:
833: . defaultValue - the default value

835:   Level: beginner

837: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
838: @*/
839: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
840: {
843:   label->defaultValue = defaultValue;
844:   return(0);
845: }

847: /*@
848:   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())

850:   Not collective

852:   Input Parameters:
853: + label - the DMLabel
854: - point - the point

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

859:   Level: intermediate

861: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
862: @*/
863: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
864: {
865:   PetscInt       v;

871:   *value = label->defaultValue;
872:   for (v = 0; v < label->numStrata; ++v) {
873:     if (label->validIS[v]) {
874:       PetscInt i;

876:       ISLocate(label->points[v], point, &i);
877:       if (i >= 0) {
878:         *value = label->stratumValues[v];
879:         break;
880:       }
881:     } else {
882:       PetscBool has;

884:       PetscHSetIHas(label->ht[v], point, &has);
885:       if (has) {
886:         *value = label->stratumValues[v];
887:         break;
888:       }
889:     }
890:   }
891:   return(0);
892: }

894: /*@
895:   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.

897:   Not collective

899:   Input Parameters:
900: + label - the DMLabel
901: . point - the point
902: - value - The point value

904:   Level: intermediate

906: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
907: @*/
908: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
909: {
910:   PetscInt       v;

915:   /* Find label value, add new entry if needed */
916:   if (value == label->defaultValue) return(0);
917:   DMLabelLookupAddStratum(label, value, &v);
918:   /* Set key */
919:   DMLabelMakeInvalid_Private(label, v);
920:   PetscHSetIAdd(label->ht[v], point);
921:   return(0);
922: }

924: /*@
925:   DMLabelClearValue - Clear the value a label assigns to a point

927:   Not collective

929:   Input Parameters:
930: + label - the DMLabel
931: . point - the point
932: - value - The point value

934:   Level: intermediate

936: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
937: @*/
938: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
939: {
940:   PetscInt       v;

945:   /* Find label value */
946:   DMLabelLookupStratum(label, value, &v);
947:   if (v < 0) return(0);

949:   if (label->bt) {
950:     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);
951:     PetscBTClear(label->bt, point - label->pStart);
952:   }

954:   /* Delete key */
955:   DMLabelMakeInvalid_Private(label, v);
956:   PetscHSetIDel(label->ht[v], point);
957:   return(0);
958: }

960: /*@
961:   DMLabelInsertIS - Set all points in the IS to a value

963:   Not collective

965:   Input Parameters:
966: + label - the DMLabel
967: . is    - the point IS
968: - value - The point value

970:   Level: intermediate

972: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
973: @*/
974: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
975: {
976:   PetscInt        v, n, p;
977:   const PetscInt *points;
978:   PetscErrorCode  ierr;

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

995: /*@
996:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

998:   Not collective

1000:   Input Parameter:
1001: . label - the DMLabel

1003:   Output Paramater:
1004: . numValues - the number of values

1006:   Level: intermediate

1008: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1009: @*/
1010: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1011: {
1015:   *numValues = label->numStrata;
1016:   return(0);
1017: }

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

1022:   Not collective

1024:   Input Parameter:
1025: . label - the DMLabel

1027:   Output Paramater:
1028: . is    - the value IS

1030:   Level: intermediate

1032: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1033: @*/
1034: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1035: {

1041:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1042:   return(0);
1043: }

1045: /*@
1046:   DMLabelHasStratum - Determine whether points exist with the given value

1048:   Not collective

1050:   Input Parameters:
1051: + label - the DMLabel
1052: - value - the stratum value

1054:   Output Paramater:
1055: . exists - Flag saying whether points exist

1057:   Level: intermediate

1059: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1060: @*/
1061: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1062: {
1063:   PetscInt       v;

1069:   DMLabelLookupStratum(label, value, &v);
1070:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1071:   return(0);
1072: }

1074: /*@
1075:   DMLabelGetStratumSize - Get the size of a stratum

1077:   Not collective

1079:   Input Parameters:
1080: + label - the DMLabel
1081: - value - the stratum value

1083:   Output Paramater:
1084: . size - The number of points in the stratum

1086:   Level: intermediate

1088: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1089: @*/
1090: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1091: {
1092:   PetscInt       v;

1098:   *size = 0;
1099:   DMLabelLookupStratum(label, value, &v);
1100:   if (v < 0) return(0);
1101:   DMLabelMakeValid_Private(label, v);
1102:   *size = label->stratumSizes[v];
1103:   return(0);
1104: }

1106: /*@
1107:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1109:   Not collective

1111:   Input Parameters:
1112: + label - the DMLabel
1113: - value - the stratum value

1115:   Output Paramaters:
1116: + start - the smallest point in the stratum
1117: - end - the largest point in the stratum

1119:   Level: intermediate

1121: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1122: @*/
1123: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1124: {
1125:   PetscInt       v, min, max;

1132:   DMLabelLookupStratum(label, value, &v);
1133:   if (v < 0) return(0);
1134:   DMLabelMakeValid_Private(label, v);
1135:   if (label->stratumSizes[v] <= 0) return(0);
1136:   ISGetMinMax(label->points[v], &min, &max);
1137:   if (start) *start = min;
1138:   if (end)   *end   = max+1;
1139:   return(0);
1140: }

1142: /*@
1143:   DMLabelGetStratumIS - Get an IS with the stratum points

1145:   Not collective

1147:   Input Parameters:
1148: + label - the DMLabel
1149: - value - the stratum value

1151:   Output Paramater:
1152: . points - The stratum points

1154:   Level: intermediate

1156:   Notes:
1157:   The output IS should be destroyed when no longer needed.
1158:   Returns NULL if the stratum is empty.

1160: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1161: @*/
1162: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1163: {
1164:   PetscInt       v;

1170:   *points = NULL;
1171:   DMLabelLookupStratum(label, value, &v);
1172:   if (v < 0) return(0);
1173:   DMLabelMakeValid_Private(label, v);
1174:   PetscObjectReference((PetscObject) label->points[v]);
1175:   *points = label->points[v];
1176:   return(0);
1177: }

1179: /*@
1180:   DMLabelSetStratumIS - Set the stratum points using an IS

1182:   Not collective

1184:   Input Parameters:
1185: + label - the DMLabel
1186: . value - the stratum value
1187: - points - The stratum points

1189:   Level: intermediate

1191: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1192: @*/
1193: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1194: {
1195:   PetscInt       v;

1201:   DMLabelLookupAddStratum(label, value, &v);
1202:   if (is == label->points[v]) return(0);
1203:   DMLabelClearStratum(label, value);
1204:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1205:   PetscObjectReference((PetscObject)is);
1206:   ISDestroy(&(label->points[v]));
1207:   label->points[v]  = is;
1208:   label->validIS[v] = PETSC_TRUE;
1209:   PetscObjectStateIncrease((PetscObject) label);
1210:   if (label->bt) {
1211:     const PetscInt *points;
1212:     PetscInt p;

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

1218:       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);
1219:       PetscBTSet(label->bt, point - label->pStart);
1220:     }
1221:   }
1222:   return(0);
1223: }

1225: /*@
1226:   DMLabelClearStratum - Remove a stratum

1228:   Not collective

1230:   Input Parameters:
1231: + label - the DMLabel
1232: - value - the stratum value

1234:   Level: intermediate

1236: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1237: @*/
1238: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1239: {
1240:   PetscInt       v;

1245:   DMLabelLookupStratum(label, value, &v);
1246:   if (v < 0) return(0);
1247:   if (label->validIS[v]) {
1248:     if (label->bt) {
1249:       PetscInt       i;
1250:       const PetscInt *points;

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

1256:         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);
1257:         PetscBTClear(label->bt, point - label->pStart);
1258:       }
1259:       ISRestoreIndices(label->points[v], &points);
1260:     }
1261:     label->stratumSizes[v] = 0;
1262:     ISDestroy(&label->points[v]);
1263:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1264:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1265:     PetscObjectStateIncrease((PetscObject) label);
1266:   } else {
1267:     PetscHSetIClear(label->ht[v]);
1268:   }
1269:   return(0);
1270: }

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

1275:   Not collective

1277:   Input Parameters:
1278: + label - the DMLabel
1279: . start - the first point kept
1280: - end - one more than the last point kept

1282:   Level: intermediate

1284: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1285: @*/
1286: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1287: {
1288:   PetscInt       v;

1293:   DMLabelDestroyIndex(label);
1294:   DMLabelMakeAllValid_Private(label);
1295:   for (v = 0; v < label->numStrata; ++v) {
1296:     ISGeneralFilter(label->points[v], start, end);
1297:   }
1298:   DMLabelCreateIndex(label, start, end);
1299:   return(0);
1300: }

1302: /*@
1303:   DMLabelPermute - Create a new label with permuted points

1305:   Not collective

1307:   Input Parameters:
1308: + label - the DMLabel
1309: - permutation - the point permutation

1311:   Output Parameter:
1312: . labelnew - the new label containing the permuted points

1314:   Level: intermediate

1316: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1317: @*/
1318: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1319: {
1320:   const PetscInt *perm;
1321:   PetscInt        numValues, numPoints, v, q;
1322:   PetscErrorCode  ierr;

1327:   DMLabelMakeAllValid_Private(label);
1328:   DMLabelDuplicate(label, labelNew);
1329:   DMLabelGetNumValues(*labelNew, &numValues);
1330:   ISGetLocalSize(permutation, &numPoints);
1331:   ISGetIndices(permutation, &perm);
1332:   for (v = 0; v < numValues; ++v) {
1333:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1334:     const PetscInt *points;
1335:     PetscInt *pointsNew;

1337:     ISGetIndices((*labelNew)->points[v],&points);
1338:     PetscMalloc1(size,&pointsNew);
1339:     for (q = 0; q < size; ++q) {
1340:       const PetscInt point = points[q];

1342:       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);
1343:       pointsNew[q] = perm[point];
1344:     }
1345:     ISRestoreIndices((*labelNew)->points[v],&points);
1346:     PetscSortInt(size, pointsNew);
1347:     ISDestroy(&((*labelNew)->points[v]));
1348:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1349:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1350:       PetscFree(pointsNew);
1351:     } else {
1352:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1353:     }
1354:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1355:   }
1356:   ISRestoreIndices(permutation, &perm);
1357:   if (label->bt) {
1358:     PetscBTDestroy(&label->bt);
1359:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1360:   }
1361:   return(0);
1362: }

1364: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1365: {
1366:   MPI_Comm       comm;
1367:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1368:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1369:   PetscSection   rootSection;
1370:   PetscSF        labelSF;

1374:   if (label) {DMLabelMakeAllValid_Private(label);}
1375:   PetscObjectGetComm((PetscObject)sf, &comm);
1376:   /* Build a section of stratum values per point, generate the according SF
1377:      and distribute point-wise stratum values to leaves. */
1378:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1379:   PetscSectionCreate(comm, &rootSection);
1380:   PetscSectionSetChart(rootSection, 0, nroots);
1381:   if (label) {
1382:     for (s = 0; s < label->numStrata; ++s) {
1383:       const PetscInt *points;

1385:       ISGetIndices(label->points[s], &points);
1386:       for (l = 0; l < label->stratumSizes[s]; l++) {
1387:         PetscSectionGetDof(rootSection, points[l], &dof);
1388:         PetscSectionSetDof(rootSection, points[l], dof+1);
1389:       }
1390:       ISRestoreIndices(label->points[s], &points);
1391:     }
1392:   }
1393:   PetscSectionSetUp(rootSection);
1394:   /* Create a point-wise array of stratum values */
1395:   PetscSectionGetStorageSize(rootSection, &size);
1396:   PetscMalloc1(size, &rootStrata);
1397:   PetscCalloc1(nroots, &rootIdx);
1398:   if (label) {
1399:     for (s = 0; s < label->numStrata; ++s) {
1400:       const PetscInt *points;

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

1429: /*@
1430:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1432:   Collective on sf

1434:   Input Parameters:
1435: + label - the DMLabel
1436: - sf    - the map from old to new distribution

1438:   Output Parameter:
1439: . labelnew - the new redistributed label

1441:   Level: intermediate

1443: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1444: @*/
1445: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1446: {
1447:   MPI_Comm       comm;
1448:   PetscSection   leafSection;
1449:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1450:   PetscInt      *leafStrata, *strataIdx;
1451:   PetscInt     **points;
1452:   const char    *lname = NULL;
1453:   char          *name;
1454:   PetscInt       nameSize;
1455:   PetscHSetI     stratumHash;
1456:   size_t         len = 0;
1457:   PetscMPIInt    rank;

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

1545: /*@
1546:   DMLabelGather - Gather all label values from leafs into roots

1548:   Collective on sf

1550:   Input Parameters:
1551: + label - the DMLabel
1552: - sf - the Star Forest point communication map

1554:   Output Parameters:
1555: . labelNew - the new DMLabel with localised leaf values

1557:   Level: developer

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

1561: .seealso: DMLabelDistribute()
1562: @*/
1563: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1564: {
1565:   MPI_Comm       comm;
1566:   PetscSection   rootSection;
1567:   PetscSF        sfLabel;
1568:   PetscSFNode   *rootPoints, *leafPoints;
1569:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1570:   const PetscInt *rootDegree, *ilocal;
1571:   PetscInt       *rootStrata;
1572:   const char    *lname;
1573:   char          *name;
1574:   PetscInt       nameSize;
1575:   size_t         len = 0;
1576:   PetscMPIInt    rank, size;

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

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

1635: /*@
1636:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1638:   Not collective

1640:   Input Parameter:
1641: . label - the DMLabel

1643:   Output Parameters:
1644: + section - the section giving offsets for each stratum
1645: - is - An IS containing all the label points

1647:   Level: developer

1649: .seealso: DMLabelDistribute()
1650: @*/
1651: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1652: {
1653:   IS              vIS;
1654:   const PetscInt *values;
1655:   PetscInt       *points;
1656:   PetscInt        nV, vS = 0, vE = 0, v, N;
1657:   PetscErrorCode  ierr;

1661:   DMLabelGetNumValues(label, &nV);
1662:   DMLabelGetValueIS(label, &vIS);
1663:   ISGetIndices(vIS, &values);
1664:   if (nV) {vS = values[0]; vE = values[0]+1;}
1665:   for (v = 1; v < nV; ++v) {
1666:     vS = PetscMin(vS, values[v]);
1667:     vE = PetscMax(vE, values[v]+1);
1668:   }
1669:   PetscSectionCreate(PETSC_COMM_SELF, section);
1670:   PetscSectionSetChart(*section, vS, vE);
1671:   for (v = 0; v < nV; ++v) {
1672:     PetscInt n;

1674:     DMLabelGetStratumSize(label, values[v], &n);
1675:     PetscSectionSetDof(*section, values[v], n);
1676:   }
1677:   PetscSectionSetUp(*section);
1678:   PetscSectionGetStorageSize(*section, &N);
1679:   PetscMalloc1(N, &points);
1680:   for (v = 0; v < nV; ++v) {
1681:     IS              is;
1682:     const PetscInt *spoints;
1683:     PetscInt        dof, off, p;

1685:     PetscSectionGetDof(*section, values[v], &dof);
1686:     PetscSectionGetOffset(*section, values[v], &off);
1687:     DMLabelGetStratumIS(label, values[v], &is);
1688:     ISGetIndices(is, &spoints);
1689:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1690:     ISRestoreIndices(is, &spoints);
1691:     ISDestroy(&is);
1692:   }
1693:   ISRestoreIndices(vIS, &values);
1694:   ISDestroy(&vIS);
1695:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1696:   return(0);
1697: }

1699: /*@
1700:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1701:   the local section and an SF describing the section point overlap.

1703:   Collective on sf

1705:   Input Parameters:
1706:   + s - The PetscSection for the local field layout
1707:   . sf - The SF describing parallel layout of the section points
1708:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1709:   . label - The label specifying the points
1710:   - labelValue - The label stratum specifying the points

1712:   Output Parameter:
1713:   . gsection - The PetscSection for the global field layout

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

1717:   Level: developer

1719: .seealso: PetscSectionCreate()
1720: @*/
1721: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1722: {
1723:   PetscInt      *neg = NULL, *tmpOff = NULL;
1724:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1731:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1732:   PetscSectionGetChart(s, &pStart, &pEnd);
1733:   PetscSectionSetChart(*gsection, pStart, pEnd);
1734:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1735:   if (nroots >= 0) {
1736:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1737:     PetscCalloc1(nroots, &neg);
1738:     if (nroots > pEnd-pStart) {
1739:       PetscCalloc1(nroots, &tmpOff);
1740:     } else {
1741:       tmpOff = &(*gsection)->atlasDof[-pStart];
1742:     }
1743:   }
1744:   /* Mark ghost points with negative dof */
1745:   for (p = pStart; p < pEnd; ++p) {
1746:     PetscInt value;

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

1789: typedef struct _n_PetscSectionSym_Label
1790: {
1791:   DMLabel           label;
1792:   PetscCopyMode     *modes;
1793:   PetscInt          *sizes;
1794:   const PetscInt    ***perms;
1795:   const PetscScalar ***rots;
1796:   PetscInt          (*minMaxOrients)[2];
1797:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1798: } PetscSectionSym_Label;

1800: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1801: {
1802:   PetscInt              i, j;
1803:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1804:   PetscErrorCode        ierr;

1807:   for (i = 0; i <= sl->numStrata; i++) {
1808:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1809:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1810:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1811:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1812:       }
1813:       if (sl->perms[i]) {
1814:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1816:         PetscFree(perms);
1817:       }
1818:       if (sl->rots[i]) {
1819:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1821:         PetscFree(rots);
1822:       }
1823:     }
1824:   }
1825:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1826:   DMLabelDestroy(&sl->label);
1827:   sl->numStrata = 0;
1828:   return(0);
1829: }

1831: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1832: {

1836:   PetscSectionSymLabelReset(sym);
1837:   PetscFree(sym->data);
1838:   return(0);
1839: }

1841: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1842: {
1843:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1844:   PetscBool             isAscii;
1845:   DMLabel               label = sl->label;
1846:   const char           *name;
1847:   PetscErrorCode        ierr;

1850:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1851:   if (isAscii) {
1852:     PetscInt          i, j, k;
1853:     PetscViewerFormat format;

1855:     PetscViewerGetFormat(viewer,&format);
1856:     if (label) {
1857:       PetscViewerGetFormat(viewer,&format);
1858:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1859:         PetscViewerASCIIPushTab(viewer);
1860:         DMLabelView(label, viewer);
1861:         PetscViewerASCIIPopTab(viewer);
1862:       } else {
1863:         PetscObjectGetName((PetscObject) sl->label, &name);
1864:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
1865:       }
1866:     } else {
1867:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1868:     }
1869:     PetscViewerASCIIPushTab(viewer);
1870:     for (i = 0; i <= sl->numStrata; i++) {
1871:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1873:       if (!(sl->perms[i] || sl->rots[i])) {
1874:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1875:       } else {
1876:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1877:         PetscViewerASCIIPushTab(viewer);
1878:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1879:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1880:           PetscViewerASCIIPushTab(viewer);
1881:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1882:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1883:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1884:             } else {
1885:               PetscInt tab;

1887:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1888:               PetscViewerASCIIPushTab(viewer);
1889:               PetscViewerASCIIGetTab(viewer,&tab);
1890:               if (sl->perms[i] && sl->perms[i][j]) {
1891:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1892:                 PetscViewerASCIISetTab(viewer,0);
1893:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1894:                 PetscViewerASCIIPrintf(viewer,"\n");
1895:                 PetscViewerASCIISetTab(viewer,tab);
1896:               }
1897:               if (sl->rots[i] && sl->rots[i][j]) {
1898:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
1899:                 PetscViewerASCIISetTab(viewer,0);
1900: #if defined(PETSC_USE_COMPLEX)
1901:                 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]));}
1902: #else
1903:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1904: #endif
1905:                 PetscViewerASCIIPrintf(viewer,"\n");
1906:                 PetscViewerASCIISetTab(viewer,tab);
1907:               }
1908:               PetscViewerASCIIPopTab(viewer);
1909:             }
1910:           }
1911:           PetscViewerASCIIPopTab(viewer);
1912:         }
1913:         PetscViewerASCIIPopTab(viewer);
1914:       }
1915:     }
1916:     PetscViewerASCIIPopTab(viewer);
1917:   }
1918:   return(0);
1919: }

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

1924:   Logically collective on sym

1926:   Input parameters:
1927: + sym - the section symmetries
1928: - label - the DMLabel describing the types of points

1930:   Level: developer:

1932: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1933: @*/
1934: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1935: {
1936:   PetscSectionSym_Label *sl;
1937:   PetscErrorCode        ierr;

1941:   sl = (PetscSectionSym_Label *) sym->data;
1942:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1943:   if (label) {
1944:     sl->label = label;
1945:     PetscObjectReference((PetscObject) label);
1946:     DMLabelGetNumValues(label,&sl->numStrata);
1947:     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);
1948:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1949:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1950:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1951:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1952:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1953:   }
1954:   return(0);
1955: }

1957: /*@C
1958:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

1960:   Logically collective on sym

1962:   InputParameters:
1963: + sym       - the section symmetries
1964: . stratum   - the stratum value in the label that we are assigning symmetries for
1965: . size      - the number of dofs for points in the stratum of the label
1966: . minOrient - the smallest orientation for a point in this stratum
1967: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1968: . mode      - how sym should copy the perms and rots arrays
1969: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
1970: - 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

1972:   Level: developer

1974: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1975: @*/
1976: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1977: {
1978:   PetscSectionSym_Label *sl;
1979:   const char            *name;
1980:   PetscInt               i, j, k;
1981:   PetscErrorCode         ierr;

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

1990:     if (stratum == value) break;
1991:   }
1992:   PetscObjectGetName((PetscObject) sl->label, &name);
1993:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
1994:   sl->sizes[i] = size;
1995:   sl->modes[i] = mode;
1996:   sl->minMaxOrients[i][0] = minOrient;
1997:   sl->minMaxOrients[i][1] = maxOrient;
1998:   if (mode == PETSC_COPY_VALUES) {
1999:     if (perms) {
2000:       PetscInt    **ownPerms;

2002:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
2003:       for (j = 0; j < maxOrient-minOrient; j++) {
2004:         if (perms[j]) {
2005:           PetscMalloc1(size,&ownPerms[j]);
2006:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2007:         }
2008:       }
2009:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2010:     }
2011:     if (rots) {
2012:       PetscScalar **ownRots;

2014:       PetscCalloc1(maxOrient - minOrient,&ownRots);
2015:       for (j = 0; j < maxOrient-minOrient; j++) {
2016:         if (rots[j]) {
2017:           PetscMalloc1(size,&ownRots[j]);
2018:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2019:         }
2020:       }
2021:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2022:     }
2023:   } else {
2024:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2025:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2026:   }
2027:   return(0);
2028: }

2030: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2031: {
2032:   PetscInt              i, j, numStrata;
2033:   PetscSectionSym_Label *sl;
2034:   DMLabel               label;
2035:   PetscErrorCode        ierr;

2038:   sl = (PetscSectionSym_Label *) sym->data;
2039:   numStrata = sl->numStrata;
2040:   label     = sl->label;
2041:   for (i = 0; i < numPoints; i++) {
2042:     PetscInt point = points[2*i];
2043:     PetscInt ornt  = points[2*i+1];

2045:     for (j = 0; j < numStrata; j++) {
2046:       if (label->validIS[j]) {
2047:         PetscInt k;

2049:         ISLocate(label->points[j],point,&k);
2050:         if (k >= 0) break;
2051:       } else {
2052:         PetscBool has;

2054:         PetscHSetIHas(label->ht[j], point, &has);
2055:         if (has) break;
2056:       }
2057:     }
2058:     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);
2059:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2060:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2061:   }
2062:   return(0);
2063: }

2065: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2066: {
2067:   PetscSectionSym_Label *sl;
2068:   PetscErrorCode        ierr;

2071:   PetscNewLog(sym,&sl);
2072:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2073:   sym->ops->view      = PetscSectionSymView_Label;
2074:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
2075:   sym->data           = (void *) sl;
2076:   return(0);
2077: }

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

2082:   Collective

2084:   Input Parameters:
2085: + comm - the MPI communicator for the new symmetry
2086: - label - the label defining the strata

2088:   Output Parameters:
2089: . sym - the section symmetries

2091:   Level: developer

2093: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2094: @*/
2095: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2096: {
2097:   PetscErrorCode        ierr;

2100:   DMInitializePackage();
2101:   PetscSectionSymCreate(comm,sym);
2102:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2103:   PetscSectionSymLabelSetLabel(*sym,label);
2104:   return(0);
2105: }