Actual source code: dmlabel.c

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

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

  9:   Input parameters:
 10: + comm - The communicator, usually PETSC_COMM_SELF
 11: - name - The label name

 13:   Output parameter:
 14: . label - The DMLabel

 16:   Level: beginner

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

 26:   DMInitializePackage();

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

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

 45: /*
 46:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 48:   Input parameter:
 49: + label - The DMLabel
 50: - v - The stratum value

 52:   Output parameter:
 53: . label - The DMLabel with stratum in sorted list format

 55:   Level: developer

 57: .seealso: DMLabelCreate()
 58: */
 59: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 60: {
 61:   PetscInt       off = 0, *pointArray, p;

 64:   if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
 66:   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);
 67:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
 68:   PetscMalloc1(label->stratumSizes[v], &pointArray);
 69:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
 70:   PetscHSetIClear(label->ht[v]);
 71:   PetscSortInt(label->stratumSizes[v], pointArray);
 72:   if (label->bt) {
 73:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 74:       const PetscInt point = pointArray[p];
 75:       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);
 76:       PetscBTSet(label->bt, point - label->pStart);
 77:     }
 78:   }
 79:   ISCreateGeneral(PETSC_COMM_SELF,label->stratumSizes[v],pointArray,PETSC_OWN_POINTER,&(label->points[v]));
 80:   PetscObjectSetName((PetscObject) (label->points[v]), "indices");
 81:   label->validIS[v] = PETSC_TRUE;
 82:   PetscObjectStateIncrease((PetscObject) label);
 83:   return(0);
 84: }

 86: /*
 87:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

 89:   Input parameter:
 90: . label - The DMLabel

 92:   Output parameter:
 93: . label - The DMLabel with all strata in sorted list format

 95:   Level: developer

 97: .seealso: DMLabelCreate()
 98: */
 99: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
100: {
101:   PetscInt       v;

105:   for (v = 0; v < label->numStrata; v++){
106:     DMLabelMakeValid_Private(label, v);
107:   }
108:   return(0);
109: }

111: /*
112:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

114:   Input parameter:
115: + label - The DMLabel
116: - v - The stratum value

118:   Output parameter:
119: . label - The DMLabel with stratum in hash format

121:   Level: developer

123: .seealso: DMLabelCreate()
124: */
125: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
126: {
127:   PetscInt       p;
128:   const PetscInt *points;

131:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
133:   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);
134:   if (label->points[v]) {
135:     ISGetIndices(label->points[v], &points);
136:     for (p = 0; p < label->stratumSizes[v]; ++p) {
137:       PetscHSetIAdd(label->ht[v], points[p]);
138:     }
139:     ISRestoreIndices(label->points[v],&points);
140:     ISDestroy(&(label->points[v]));
141:   }
142:   label->validIS[v] = PETSC_FALSE;
143:   return(0);
144: }

146: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
147: #define DMLABEL_LOOKUP_THRESHOLD 16
148: #endif

150: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
151: {
152:   PetscInt       v;

156:   *index = -1;
157:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
158:     for (v = 0; v < label->numStrata; ++v)
159:       if (label->stratumValues[v] == value) {*index = v; break;}
160:   } else {
161:     PetscHMapIGet(label->hmap, value, index);
162:   }
163: #if defined(PETSC_USE_DEBUG)
164:   { /* Check strata hash map consistency */
165:     PetscInt len, loc = -1;
166:     PetscHMapIGetSize(label->hmap, &len);
167:     if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
168:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
169:       PetscHMapIGet(label->hmap, value, &loc);
170:     } else {
171:       for (v = 0; v < label->numStrata; ++v)
172:         if (label->stratumValues[v] == value) {loc = v; break;}
173:     }
174:     if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
175:   }
176: #endif
177:   return(0);
178: }

180: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
181: {
182:   PetscInt       v;
183:   PetscInt      *tmpV;
184:   PetscInt      *tmpS;
185:   PetscHSetI    *tmpH, ht;
186:   IS            *tmpP, is;
187:   PetscBool     *tmpB;
188:   PetscHMapI     hmap = label->hmap;

192:   v    = label->numStrata;
193:   tmpV = label->stratumValues;
194:   tmpS = label->stratumSizes;
195:   tmpH = label->ht;
196:   tmpP = label->points;
197:   tmpB = label->validIS;
198:   { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free  */
199:     PetscInt   *oldV = tmpV;
200:     PetscInt   *oldS = tmpS;
201:     PetscHSetI *oldH = tmpH;
202:     IS         *oldP = tmpP;
203:     PetscBool  *oldB = tmpB;
204:     PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
205:     PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
206:     PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
207:     PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
208:     PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
209:     PetscMemcpy(tmpV, oldV, v*sizeof(*tmpV));
210:     PetscMemcpy(tmpS, oldS, v*sizeof(*tmpS));
211:     PetscMemcpy(tmpH, oldH, v*sizeof(*tmpH));
212:     PetscMemcpy(tmpP, oldP, v*sizeof(*tmpP));
213:     PetscMemcpy(tmpB, oldB, v*sizeof(*tmpB));
214:     PetscFree(oldV);
215:     PetscFree(oldS);
216:     PetscFree(oldH);
217:     PetscFree(oldP);
218:     PetscFree(oldB);
219:   }
220:   label->numStrata     = v+1;
221:   label->stratumValues = tmpV;
222:   label->stratumSizes  = tmpS;
223:   label->ht            = tmpH;
224:   label->points        = tmpP;
225:   label->validIS       = tmpB;
226:   PetscHSetICreate(&ht);
227:   ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
228:   PetscHMapISet(hmap, value, v);
229:   tmpV[v] = value;
230:   tmpS[v] = 0;
231:   tmpH[v] = ht;
232:   tmpP[v] = is;
233:   tmpB[v] = PETSC_TRUE;
234:   *index = v;
235:   return(0);
236: }

238: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
239: {
242:   DMLabelLookupStratum(label, value, index);
243:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
244:   return(0);
245: }

247: /*@
248:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

250:   Input Parameter:
251: + label - The DMLabel
252: - value - The stratum value

254:   Level: beginner

256: .seealso:  DMLabelCreate(), DMLabelDestroy()
257: @*/
258: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
259: {
260:   PetscInt       v;

265:   DMLabelLookupAddStratum(label, value, &v);
266:   return(0);
267: }

269: /*@
270:   DMLabelAddStrata - Adds new stratum values in a DMLabel

272:   Input Parameter:
273: + label - The DMLabel
274: . numStrata - The number of stratum values
275: - stratumValues - The stratum values

277:   Level: beginner

279: .seealso:  DMLabelCreate(), DMLabelDestroy()
280: @*/
281: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
282: {
283:   PetscInt       *values, v;

289:   PetscMalloc1(numStrata, &values);
290:   PetscMemcpy(values, stratumValues, numStrata*sizeof(PetscInt));
291:   PetscSortRemoveDupsInt(&numStrata, values);
292:   if (!label->numStrata) { /* Fast preallocation */
293:     PetscInt   *tmpV;
294:     PetscInt   *tmpS;
295:     PetscHSetI *tmpH, ht;
296:     IS         *tmpP, is;
297:     PetscBool  *tmpB;
298:     PetscHMapI  hmap = label->hmap;

300:     PetscMalloc1(numStrata, &tmpV);
301:     PetscMalloc1(numStrata, &tmpS);
302:     PetscMalloc1(numStrata, &tmpH);
303:     PetscMalloc1(numStrata, &tmpP);
304:     PetscMalloc1(numStrata, &tmpB);
305:     label->numStrata     = numStrata;
306:     label->stratumValues = tmpV;
307:     label->stratumSizes  = tmpS;
308:     label->ht            = tmpH;
309:     label->points        = tmpP;
310:     label->validIS       = tmpB;
311:     for (v = 0; v < numStrata; ++v) {
312:       PetscHSetICreate(&ht);
313:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
314:       PetscHMapISet(hmap, values[v], v);
315:       tmpV[v] = values[v];
316:       tmpS[v] = 0;
317:       tmpH[v] = ht;
318:       tmpP[v] = is;
319:       tmpB[v] = PETSC_TRUE;
320:     }
321:   } else {
322:     for (v = 0; v < numStrata; ++v) {
323:       DMLabelAddStratum(label, values[v]);
324:     }
325:   }
326:   PetscFree(values);
327:   return(0);
328: }

330: /*@
331:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

333:   Input Parameter:
334: + label - The DMLabel
335: - valueIS - Index set with stratum values

337:   Level: beginner

339: .seealso:  DMLabelCreate(), DMLabelDestroy()
340: @*/
341: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
342: {
343:   PetscInt       numStrata;
344:   const PetscInt *stratumValues;

350:   ISGetLocalSize(valueIS, &numStrata);
351:   ISGetIndices(valueIS, &stratumValues);
352:   DMLabelAddStrata(label, numStrata, stratumValues);
353:   return(0);
354: }

356: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
357: {
358:   PetscInt       v;
359:   PetscMPIInt    rank;

363:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
364:   PetscViewerASCIIPushSynchronized(viewer);
365:   if (label) {
366:     const char *name;

368:     PetscObjectGetName((PetscObject) label, &name);
369:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
370:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
371:     for (v = 0; v < label->numStrata; ++v) {
372:       const PetscInt value = label->stratumValues[v];
373:       const PetscInt *points;
374:       PetscInt       p;

376:       ISGetIndices(label->points[v], &points);
377:       for (p = 0; p < label->stratumSizes[v]; ++p) {
378:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
379:       }
380:       ISRestoreIndices(label->points[v],&points);
381:     }
382:   }
383:   PetscViewerFlush(viewer);
384:   PetscViewerASCIIPopSynchronized(viewer);
385:   return(0);
386: }

388: /*@C
389:   DMLabelView - View the label

391:   Input Parameters:
392: + label - The DMLabel
393: - viewer - The PetscViewer

395:   Level: intermediate

397: .seealso: DMLabelCreate(), DMLabelDestroy()
398: @*/
399: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
400: {
401:   PetscBool      iascii;

406:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
408:   if (label) {DMLabelMakeAllValid_Private(label);}
409:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
410:   if (iascii) {
411:     DMLabelView_Ascii(label, viewer);
412:   }
413:   return(0);
414: }

416: /*@
417:   DMLabelReset - Destroys internal data structures in a DMLabel

419:   Input Parameter:
420: . label - The DMLabel

422:   Level: beginner

424: .seealso: DMLabelDestroy(), DMLabelCreate()
425: @*/
426: PetscErrorCode DMLabelReset(DMLabel label)
427: {
428:   PetscInt       v;

433:   for (v = 0; v < label->numStrata; ++v) {
434:     PetscHSetIDestroy(&label->ht[v]);
435:     ISDestroy(&label->points[v]);
436:   }
437:   label->numStrata = 0;
438:   PetscFree(label->stratumValues);
439:   PetscFree(label->stratumSizes);
440:   PetscFree(label->ht);
441:   PetscFree(label->points);
442:   PetscFree(label->validIS);
443:   PetscHMapIReset(label->hmap);
444:   label->pStart = -1;
445:   label->pEnd   = -1;
446:   PetscBTDestroy(&label->bt);
447:   return(0);
448: }

450: /*@
451:   DMLabelDestroy - Destroys a DMLabel

453:   Input Parameter:
454: . label - The DMLabel

456:   Level: beginner

458: .seealso: DMLabelReset(), DMLabelCreate()
459: @*/
460: PetscErrorCode DMLabelDestroy(DMLabel *label)
461: {

465:   if (!*label) return(0);
467:   if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
468:   DMLabelReset(*label);
469:   PetscHMapIDestroy(&(*label)->hmap);
470:   PetscHeaderDestroy(label);
471:   return(0);
472: }

474: /*@
475:   DMLabelDuplicate - Duplicates a DMLabel

477:   Input Parameter:
478: . label - The DMLabel

480:   Output Parameter:
481: . labelnew - location to put new vector

483:   Level: intermediate

485: .seealso: DMLabelCreate(), DMLabelDestroy()
486: @*/
487: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
488: {
489:   const char    *name;
490:   PetscInt       v;

495:   DMLabelMakeAllValid_Private(label);
496:   PetscObjectGetName((PetscObject) label, &name);
497:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

499:   (*labelnew)->numStrata    = label->numStrata;
500:   (*labelnew)->defaultValue = label->defaultValue;
501:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
502:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
503:   PetscMalloc1(label->numStrata, &(*labelnew)->ht);
504:   PetscMalloc1(label->numStrata, &(*labelnew)->points);
505:   PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
506:   for (v = 0; v < label->numStrata; ++v) {
507:     PetscHSetICreate(&(*labelnew)->ht[v]);
508:     (*labelnew)->stratumValues[v]  = label->stratumValues[v];
509:     (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
510:     PetscObjectReference((PetscObject) (label->points[v]));
511:     (*labelnew)->points[v]         = label->points[v];
512:     (*labelnew)->validIS[v]        = PETSC_TRUE;
513:   }
514:   PetscHMapIDestroy(&(*labelnew)->hmap);
515:   PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
516:   (*labelnew)->pStart = -1;
517:   (*labelnew)->pEnd   = -1;
518:   (*labelnew)->bt     = NULL;
519:   return(0);
520: }

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

525:   Input Parameter:
526: . label  - The DMLabel

528:   Level: intermediate

530: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
531: @*/
532: PetscErrorCode DMLabelComputeIndex(DMLabel label)
533: {
534:   PetscInt       pStart = PETSC_MAX_INT, pEnd = -1, v;

539:   DMLabelMakeAllValid_Private(label);
540:   for (v = 0; v < label->numStrata; ++v) {
541:     const PetscInt *points;
542:     PetscInt       i;

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

548:       pStart = PetscMin(point,   pStart);
549:       pEnd   = PetscMax(point+1, pEnd);
550:     }
551:     ISRestoreIndices(label->points[v], &points);
552:   }
553:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
554:   label->pEnd   = pEnd;
555:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
556:   return(0);
557: }

559: /*@
560:   DMLabelCreateIndex - Create an index structure for membership determination

562:   Input Parameters:
563: + label  - The DMLabel
564: . pStart - The smallest point
565: - pEnd   - The largest point + 1

567:   Level: intermediate

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

578:   DMLabelDestroyIndex(label);
579:   DMLabelMakeAllValid_Private(label);
580:   label->pStart = pStart;
581:   label->pEnd   = pEnd;
582:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
583:   PetscBTCreate(pEnd - pStart, &label->bt);
584:   for (v = 0; v < label->numStrata; ++v) {
585:     const PetscInt *points;
586:     PetscInt       i;

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

592:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
593:       PetscBTSet(label->bt, point - pStart);
594:     }
595:     ISRestoreIndices(label->points[v], &points);
596:   }
597:   return(0);
598: }

600: /*@
601:   DMLabelDestroyIndex - Destroy the index structure

603:   Input Parameter:
604: . label - the DMLabel

606:   Level: intermediate

608: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
609: @*/
610: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
611: {

616:   label->pStart = -1;
617:   label->pEnd   = -1;
618:   PetscBTDestroy(&label->bt);
619:   return(0);
620: }

622: /*@
623:   DMLabelGetBounds - Return the smallest and largest point in the label

625:   Input Parameter:
626: . label - the DMLabel

628:   Output Parameters:
629: + pStart - The smallest point
630: - pEnd   - The largest point + 1

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

634:   Level: intermediate

636: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
637: @*/
638: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
639: {

644:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
645:   if (pStart) {
647:     *pStart = label->pStart;
648:   }
649:   if (pEnd) {
651:     *pEnd = label->pEnd;
652:   }
653:   return(0);
654: }

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

659:   Input Parameters:
660: + label - the DMLabel
661: - value - the value

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

666:   Level: developer

668: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
669: @*/
670: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
671: {
672:   PetscInt v;

678:   DMLabelLookupStratum(label, value, &v);
679:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
680:   return(0);
681: }

683: /*@
684:   DMLabelHasPoint - Determine whether a label assigns a value to a point

686:   Input Parameters:
687: + label - the DMLabel
688: - point - the point

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

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

695:   Level: developer

697: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
698: @*/
699: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
700: {

706:   DMLabelMakeAllValid_Private(label);
707: #if defined(PETSC_USE_DEBUG)
708:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
709:   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);
710: #endif
711:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
712:   return(0);
713: }

715: /*@
716:   DMLabelStratumHasPoint - Return true if the stratum contains a point

718:   Input Parameters:
719: + label - the DMLabel
720: . value - the stratum value
721: - point - the point

723:   Output Parameter:
724: . contains - true if the stratum contains the point

726:   Level: intermediate

728: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
729: @*/
730: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
731: {
732:   PetscInt       v;

738:   *contains = PETSC_FALSE;
739:   DMLabelLookupStratum(label, value, &v);
740:   if (v < 0) return(0);

742:   if (label->validIS[v]) {
743:     PetscInt i;

745:     ISLocate(label->points[v], point, &i);
746:     if (i >= 0) *contains = PETSC_TRUE;
747:   } else {
748:     PetscBool has;

750:     PetscHSetIHas(label->ht[v], point, &has);
751:     if (has) *contains = PETSC_TRUE;
752:   }
753:   return(0);
754: }

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

760:   Input parameter:
761: . label - a DMLabel object

763:   Output parameter:
764: . defaultValue - the default value

766:   Level: beginner

768: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
769: @*/
770: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
771: {
774:   *defaultValue = label->defaultValue;
775:   return(0);
776: }

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

782:   Input parameter:
783: . label - a DMLabel object

785:   Output parameter:
786: . defaultValue - the default value

788:   Level: beginner

790: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
791: @*/
792: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
793: {
796:   label->defaultValue = defaultValue;
797:   return(0);
798: }

800: /*@
801:   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())

803:   Input Parameters:
804: + label - the DMLabel
805: - point - the point

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

810:   Level: intermediate

812: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
813: @*/
814: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
815: {
816:   PetscInt       v;

822:   *value = label->defaultValue;
823:   for (v = 0; v < label->numStrata; ++v) {
824:     if (label->validIS[v]) {
825:       PetscInt i;

827:       ISLocate(label->points[v], point, &i);
828:       if (i >= 0) {
829:         *value = label->stratumValues[v];
830:         break;
831:       }
832:     } else {
833:       PetscBool has;

835:       PetscHSetIHas(label->ht[v], point, &has);
836:       if (has) {
837:         *value = label->stratumValues[v];
838:         break;
839:       }
840:     }
841:   }
842:   return(0);
843: }

845: /*@
846:   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.

848:   Input Parameters:
849: + label - the DMLabel
850: . point - the point
851: - value - The point value

853:   Level: intermediate

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

864:   /* Find label value, add new entry if needed */
865:   if (value == label->defaultValue) return(0);
866:   DMLabelLookupAddStratum(label, value, &v);
867:   /* Set key */
868:   DMLabelMakeInvalid_Private(label, v);
869:   PetscHSetIAdd(label->ht[v], point);
870:   return(0);
871: }

873: /*@
874:   DMLabelClearValue - Clear the value a label assigns to a point

876:   Input Parameters:
877: + label - the DMLabel
878: . point - the point
879: - value - The point value

881:   Level: intermediate

883: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
884: @*/
885: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
886: {
887:   PetscInt       v;

892:   /* Find label value */
893:   DMLabelLookupStratum(label, value, &v);
894:   if (v < 0) return(0);

896:   if (label->bt) {
897:     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);
898:     PetscBTClear(label->bt, point - label->pStart);
899:   }

901:   /* Delete key */
902:   DMLabelMakeInvalid_Private(label, v);
903:   PetscHSetIDel(label->ht[v], point);
904:   return(0);
905: }

907: /*@
908:   DMLabelInsertIS - Set all points in the IS to a value

910:   Input Parameters:
911: + label - the DMLabel
912: . is    - the point IS
913: - value - The point value

915:   Level: intermediate

917: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
918: @*/
919: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
920: {
921:   PetscInt        v, n, p;
922:   const PetscInt *points;
923:   PetscErrorCode  ierr;

928:   /* Find label value, add new entry if needed */
929:   if (value == label->defaultValue) return(0);
930:   DMLabelLookupAddStratum(label, value, &v);
931:   /* Set keys */
932:   DMLabelMakeInvalid_Private(label, v);
933:   ISGetLocalSize(is, &n);
934:   ISGetIndices(is, &points);
935:   for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
936:   ISRestoreIndices(is, &points);
937:   return(0);
938: }

940: /*@
941:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

943:   Input Parameter:
944: . label - the DMLabel

946:   Output Paramater:
947: . numValues - the number of values

949:   Level: intermediate

951: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
952: @*/
953: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
954: {
958:   *numValues = label->numStrata;
959:   return(0);
960: }

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

965:   Input Parameter:
966: . label - the DMLabel

968:   Output Paramater:
969: . is    - the value IS

971:   Level: intermediate

973: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
974: @*/
975: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
976: {

982:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
983:   return(0);
984: }

986: /*@
987:   DMLabelHasStratum - Determine whether points exist with the given value

989:   Input Parameters:
990: + label - the DMLabel
991: - value - the stratum value

993:   Output Paramater:
994: . exists - Flag saying whether points exist

996:   Level: intermediate

998: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
999: @*/
1000: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1001: {
1002:   PetscInt       v;

1008:   DMLabelLookupStratum(label, value, &v);
1009:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1010:   return(0);
1011: }

1013: /*@
1014:   DMLabelGetStratumSize - Get the size of a stratum

1016:   Input Parameters:
1017: + label - the DMLabel
1018: - value - the stratum value

1020:   Output Paramater:
1021: . size - The number of points in the stratum

1023:   Level: intermediate

1025: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1026: @*/
1027: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1028: {
1029:   PetscInt       v;

1035:   *size = 0;
1036:   DMLabelLookupStratum(label, value, &v);
1037:   if (v < 0) return(0);
1038:   DMLabelMakeValid_Private(label, v);
1039:   *size = label->stratumSizes[v];
1040:   return(0);
1041: }

1043: /*@
1044:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

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

1050:   Output Paramaters:
1051: + start - the smallest point in the stratum
1052: - end - the largest point in the stratum

1054:   Level: intermediate

1056: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1057: @*/
1058: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1059: {
1060:   PetscInt       v, min, max;

1067:   DMLabelLookupStratum(label, value, &v);
1068:   if (v < 0) return(0);
1069:   DMLabelMakeValid_Private(label, v);
1070:   if (label->stratumSizes[v] <= 0) return(0);
1071:   ISGetMinMax(label->points[v], &min, &max);
1072:   if (start) *start = min;
1073:   if (end)   *end   = max+1;
1074:   return(0);
1075: }

1077: /*@
1078:   DMLabelGetStratumIS - Get an IS with the stratum points

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

1084:   Output Paramater:
1085: . points - The stratum points

1087:   Level: intermediate

1089: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1090: @*/
1091: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1092: {
1093:   PetscInt       v;

1099:   *points = NULL;
1100:   DMLabelLookupStratum(label, value, &v);
1101:   if (v < 0) return(0);
1102:   DMLabelMakeValid_Private(label, v);
1103:   PetscObjectReference((PetscObject) label->points[v]);
1104:   *points = label->points[v];
1105:   return(0);
1106: }

1108: /*@
1109:   DMLabelSetStratumIS - Set the stratum points using an IS

1111:   Input Parameters:
1112: + label - the DMLabel
1113: . value - the stratum value
1114: - points - The stratum points

1116:   Level: intermediate

1118: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1119: @*/
1120: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1121: {
1122:   PetscInt       v;

1128:   DMLabelLookupAddStratum(label, value, &v);
1129:   if (is == label->points[v]) return(0);
1130:   DMLabelClearStratum(label, value);
1131:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1132:   PetscObjectReference((PetscObject)is);
1133:   ISDestroy(&(label->points[v]));
1134:   label->points[v] = is;
1135:   label->validIS[v] = PETSC_TRUE;
1136:   if (label->bt) {
1137:     const PetscInt *points;
1138:     PetscInt p;

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

1144:       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);
1145:       PetscBTSet(label->bt, point - label->pStart);
1146:     }
1147:   }
1148:   return(0);
1149: }

1151: /*@
1152:   DMLabelClearStratum - Remove a stratum

1154:   Input Parameters:
1155: + label - the DMLabel
1156: - value - the stratum value

1158:   Level: intermediate

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

1169:   DMLabelLookupStratum(label, value, &v);
1170:   if (v < 0) return(0);
1171:   if (label->validIS[v]) {
1172:     if (label->bt) {
1173:       PetscInt       i;
1174:       const PetscInt *points;

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

1180:         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);
1181:         PetscBTClear(label->bt, point - label->pStart);
1182:       }
1183:       ISRestoreIndices(label->points[v], &points);
1184:     }
1185:     label->stratumSizes[v] = 0;
1186:     ISDestroy(&label->points[v]);
1187:     ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_OWN_POINTER, &label->points[v]);
1188:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1189:   } else {
1190:     PetscHSetIClear(label->ht[v]);
1191:   }
1192:   return(0);
1193: }

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

1198:   Input Parameters:
1199: + label - the DMLabel
1200: . start - the first point
1201: - end - the last point

1203:   Level: intermediate

1205: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1206: @*/
1207: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1208: {
1209:   PetscInt       v;

1214:   DMLabelDestroyIndex(label);
1215:   DMLabelMakeAllValid_Private(label);
1216:   for (v = 0; v < label->numStrata; ++v) {
1217:     PetscInt off, q;
1218:     const PetscInt *points;
1219:     PetscInt numPointsNew = 0, *pointsNew = NULL;

1221:     ISGetIndices(label->points[v], &points);
1222:     for (q = 0; q < label->stratumSizes[v]; ++q)
1223:       if (points[q] >= start && points[q] < end)
1224:         numPointsNew++;
1225:     PetscMalloc1(numPointsNew, &pointsNew);
1226:     for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
1227:       if (points[q] >= start && points[q] < end)
1228:         pointsNew[off++] = points[q];
1229:     }
1230:     ISRestoreIndices(label->points[v],&points);

1232:     label->stratumSizes[v] = numPointsNew;
1233:     ISDestroy(&label->points[v]);
1234:     ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
1235:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1236:   }
1237:   DMLabelCreateIndex(label, start, end);
1238:   return(0);
1239: }

1241: /*@
1242:   DMLabelPermute - Create a new label with permuted points

1244:   Input Parameters:
1245: + label - the DMLabel
1246: - permutation - the point permutation

1248:   Output Parameter:
1249: . labelnew - the new label containing the permuted points

1251:   Level: intermediate

1253: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1254: @*/
1255: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1256: {
1257:   const PetscInt *perm;
1258:   PetscInt        numValues, numPoints, v, q;
1259:   PetscErrorCode  ierr;

1264:   DMLabelMakeAllValid_Private(label);
1265:   DMLabelDuplicate(label, labelNew);
1266:   DMLabelGetNumValues(*labelNew, &numValues);
1267:   ISGetLocalSize(permutation, &numPoints);
1268:   ISGetIndices(permutation, &perm);
1269:   for (v = 0; v < numValues; ++v) {
1270:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1271:     const PetscInt *points;
1272:     PetscInt *pointsNew;

1274:     ISGetIndices((*labelNew)->points[v],&points);
1275:     PetscMalloc1(size,&pointsNew);
1276:     for (q = 0; q < size; ++q) {
1277:       const PetscInt point = points[q];

1279:       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);
1280:       pointsNew[q] = perm[point];
1281:     }
1282:     ISRestoreIndices((*labelNew)->points[v],&points);
1283:     PetscSortInt(size, pointsNew);
1284:     ISDestroy(&((*labelNew)->points[v]));
1285:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1286:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1287:       PetscFree(pointsNew);
1288:     } else {
1289:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1290:     }
1291:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1292:   }
1293:   ISRestoreIndices(permutation, &perm);
1294:   if (label->bt) {
1295:     PetscBTDestroy(&label->bt);
1296:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1297:   }
1298:   return(0);
1299: }

1301: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1302: {
1303:   MPI_Comm       comm;
1304:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1305:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1306:   PetscSection   rootSection;
1307:   PetscSF        labelSF;

1311:   if (label) {DMLabelMakeAllValid_Private(label);}
1312:   PetscObjectGetComm((PetscObject)sf, &comm);
1313:   /* Build a section of stratum values per point, generate the according SF
1314:      and distribute point-wise stratum values to leaves. */
1315:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1316:   PetscSectionCreate(comm, &rootSection);
1317:   PetscSectionSetChart(rootSection, 0, nroots);
1318:   if (label) {
1319:     for (s = 0; s < label->numStrata; ++s) {
1320:       const PetscInt *points;

1322:       ISGetIndices(label->points[s], &points);
1323:       for (l = 0; l < label->stratumSizes[s]; l++) {
1324:         PetscSectionGetDof(rootSection, points[l], &dof);
1325:         PetscSectionSetDof(rootSection, points[l], dof+1);
1326:       }
1327:       ISRestoreIndices(label->points[s], &points);
1328:     }
1329:   }
1330:   PetscSectionSetUp(rootSection);
1331:   /* Create a point-wise array of stratum values */
1332:   PetscSectionGetStorageSize(rootSection, &size);
1333:   PetscMalloc1(size, &rootStrata);
1334:   PetscCalloc1(nroots, &rootIdx);
1335:   if (label) {
1336:     for (s = 0; s < label->numStrata; ++s) {
1337:       const PetscInt *points;

1339:       ISGetIndices(label->points[s], &points);
1340:       for (l = 0; l < label->stratumSizes[s]; l++) {
1341:         const PetscInt p = points[l];
1342:         PetscSectionGetOffset(rootSection, p, &offset);
1343:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1344:       }
1345:       ISRestoreIndices(label->points[s], &points);
1346:     }
1347:   }
1348:   /* Build SF that maps label points to remote processes */
1349:   PetscSectionCreate(comm, leafSection);
1350:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1351:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1352:   PetscFree(remoteOffsets);
1353:   /* Send the strata for each point over the derived SF */
1354:   PetscSectionGetStorageSize(*leafSection, &size);
1355:   PetscMalloc1(size, leafStrata);
1356:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1357:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1358:   /* Clean up */
1359:   PetscFree(rootStrata);
1360:   PetscFree(rootIdx);
1361:   PetscSectionDestroy(&rootSection);
1362:   PetscSFDestroy(&labelSF);
1363:   return(0);
1364: }

1366: /*@
1367:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1369:   Input Parameters:
1370: + label - the DMLabel
1371: - sf    - the map from old to new distribution

1373:   Output Parameter:
1374: . labelnew - the new redistributed label

1376:   Level: intermediate

1378: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1379: @*/
1380: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1381: {
1382:   MPI_Comm       comm;
1383:   PetscSection   leafSection;
1384:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1385:   PetscInt      *leafStrata, *strataIdx;
1386:   PetscInt     **points;
1387:   const char    *lname = NULL;
1388:   char          *name;
1389:   PetscInt       nameSize;
1390:   PetscHSetI     stratumHash;
1391:   size_t         len = 0;
1392:   PetscMPIInt    rank;

1397:   if (label) {
1399:     DMLabelMakeAllValid_Private(label);
1400:   }
1401:   PetscObjectGetComm((PetscObject)sf, &comm);
1402:   MPI_Comm_rank(comm, &rank);
1403:   /* Bcast name */
1404:   if (!rank) {
1405:     PetscObjectGetName((PetscObject) label, &lname);
1406:     PetscStrlen(lname, &len);
1407:   }
1408:   nameSize = len;
1409:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1410:   PetscMalloc1(nameSize+1, &name);
1411:   if (!rank) {PetscMemcpy(name, lname, nameSize+1);}
1412:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1413:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1414:   PetscFree(name);
1415:   /* Bcast defaultValue */
1416:   if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1417:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1418:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1419:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1420:   /* Determine received stratum values and initialise new label*/
1421:   PetscHSetICreate(&stratumHash);
1422:   PetscSectionGetStorageSize(leafSection, &size);
1423:   for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1424:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1425:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1426:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1427:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1428:   /* Turn leafStrata into indices rather than stratum values */
1429:   offset = 0;
1430:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1431:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1432:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1433:     PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1434:   }
1435:   for (p = 0; p < size; ++p) {
1436:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1437:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1438:     }
1439:   }
1440:   /* Rebuild the point strata on the receiver */
1441:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1442:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1443:   for (p=pStart; p<pEnd; p++) {
1444:     PetscSectionGetDof(leafSection, p, &dof);
1445:     PetscSectionGetOffset(leafSection, p, &offset);
1446:     for (s=0; s<dof; s++) {
1447:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1448:     }
1449:   }
1450:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1451:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1452:   PetscMalloc1((*labelNew)->numStrata,&points);
1453:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1454:     PetscHSetICreate(&(*labelNew)->ht[s]);
1455:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1456:   }
1457:   /* Insert points into new strata */
1458:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1459:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1460:   for (p=pStart; p<pEnd; p++) {
1461:     PetscSectionGetDof(leafSection, p, &dof);
1462:     PetscSectionGetOffset(leafSection, p, &offset);
1463:     for (s=0; s<dof; s++) {
1464:       stratum = leafStrata[offset+s];
1465:       points[stratum][strataIdx[stratum]++] = p;
1466:     }
1467:   }
1468:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1469:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1470:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1471:   }
1472:   PetscFree(points);
1473:   PetscHSetIDestroy(&stratumHash);
1474:   PetscFree(leafStrata);
1475:   PetscFree(strataIdx);
1476:   PetscSectionDestroy(&leafSection);
1477:   return(0);
1478: }

1480: /*@
1481:   DMLabelGather - Gather all label values from leafs into roots

1483:   Input Parameters:
1484: + label - the DMLabel
1485: - sf - the Star Forest point communication map

1487:   Output Parameters:
1488: . labelNew - the new DMLabel with localised leaf values

1490:   Level: developer

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

1494: .seealso: DMLabelDistribute()
1495: @*/
1496: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1497: {
1498:   MPI_Comm       comm;
1499:   PetscSection   rootSection;
1500:   PetscSF        sfLabel;
1501:   PetscSFNode   *rootPoints, *leafPoints;
1502:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1503:   const PetscInt *rootDegree, *ilocal;
1504:   PetscInt       *rootStrata;
1505:   const char    *lname;
1506:   char          *name;
1507:   PetscInt       nameSize;
1508:   size_t         len = 0;
1509:   PetscMPIInt    rank, size;

1515:   PetscObjectGetComm((PetscObject)sf, &comm);
1516:   MPI_Comm_rank(comm, &rank);
1517:   MPI_Comm_size(comm, &size);
1518:   /* Bcast name */
1519:   if (!rank) {
1520:     PetscObjectGetName((PetscObject) label, &lname);
1521:     PetscStrlen(lname, &len);
1522:   }
1523:   nameSize = len;
1524:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1525:   PetscMalloc1(nameSize+1, &name);
1526:   if (!rank) {PetscMemcpy(name, lname, nameSize+1);}
1527:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1528:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1529:   PetscFree(name);
1530:   /* Gather rank/index pairs of leaves into local roots to build
1531:      an inverse, multi-rooted SF. Note that this ignores local leaf
1532:      indexing due to the use of the multiSF in PetscSFGather. */
1533:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1534:   PetscMalloc1(nroots, &leafPoints);
1535:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1536:   for (p = 0; p < nleaves; p++) {
1537:     PetscInt ilp = ilocal ? ilocal[p] : p;

1539:     leafPoints[ilp].index = ilp;
1540:     leafPoints[ilp].rank  = rank;
1541:   }
1542:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1543:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1544:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1545:   PetscMalloc1(nmultiroots, &rootPoints);
1546:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1547:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1548:   PetscSFCreate(comm,& sfLabel);
1549:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1550:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1551:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1552:   /* Rebuild the point strata on the receiver */
1553:   for (p = 0, idx = 0; p < nroots; p++) {
1554:     for (d = 0; d < rootDegree[p]; d++) {
1555:       PetscSectionGetDof(rootSection, idx+d, &dof);
1556:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1557:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1558:     }
1559:     idx += rootDegree[p];
1560:   }
1561:   PetscFree(leafPoints);
1562:   PetscFree(rootStrata);
1563:   PetscSectionDestroy(&rootSection);
1564:   PetscSFDestroy(&sfLabel);
1565:   return(0);
1566: }

1568: /*@
1569:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1571:   Input Parameter:
1572: . label - the DMLabel

1574:   Output Parameters:
1575: + section - the section giving offsets for each stratum
1576: - is - An IS containing all the label points

1578:   Level: developer

1580: .seealso: DMLabelDistribute()
1581: @*/
1582: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1583: {
1584:   IS              vIS;
1585:   const PetscInt *values;
1586:   PetscInt       *points;
1587:   PetscInt        nV, vS = 0, vE = 0, v, N;
1588:   PetscErrorCode  ierr;

1592:   DMLabelGetNumValues(label, &nV);
1593:   DMLabelGetValueIS(label, &vIS);
1594:   ISGetIndices(vIS, &values);
1595:   if (nV) {vS = values[0]; vE = values[0]+1;}
1596:   for (v = 1; v < nV; ++v) {
1597:     vS = PetscMin(vS, values[v]);
1598:     vE = PetscMax(vE, values[v]+1);
1599:   }
1600:   PetscSectionCreate(PETSC_COMM_SELF, section);
1601:   PetscSectionSetChart(*section, vS, vE);
1602:   for (v = 0; v < nV; ++v) {
1603:     PetscInt n;

1605:     DMLabelGetStratumSize(label, values[v], &n);
1606:     PetscSectionSetDof(*section, values[v], n);
1607:   }
1608:   PetscSectionSetUp(*section);
1609:   PetscSectionGetStorageSize(*section, &N);
1610:   PetscMalloc1(N, &points);
1611:   for (v = 0; v < nV; ++v) {
1612:     IS              is;
1613:     const PetscInt *spoints;
1614:     PetscInt        dof, off, p;

1616:     PetscSectionGetDof(*section, values[v], &dof);
1617:     PetscSectionGetOffset(*section, values[v], &off);
1618:     DMLabelGetStratumIS(label, values[v], &is);
1619:     ISGetIndices(is, &spoints);
1620:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1621:     ISRestoreIndices(is, &spoints);
1622:     ISDestroy(&is);
1623:   }
1624:   ISRestoreIndices(vIS, &values);
1625:   ISDestroy(&vIS);
1626:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1627:   return(0);
1628: }

1630: /*@
1631:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1632:   the local section and an SF describing the section point overlap.

1634:   Input Parameters:
1635:   + s - The PetscSection for the local field layout
1636:   . sf - The SF describing parallel layout of the section points
1637:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1638:   . label - The label specifying the points
1639:   - labelValue - The label stratum specifying the points

1641:   Output Parameter:
1642:   . gsection - The PetscSection for the global field layout

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

1646:   Level: developer

1648: .seealso: PetscSectionCreate()
1649: @*/
1650: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1651: {
1652:   PetscInt      *neg = NULL, *tmpOff = NULL;
1653:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1660:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1661:   PetscSectionGetChart(s, &pStart, &pEnd);
1662:   PetscSectionSetChart(*gsection, pStart, pEnd);
1663:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1664:   if (nroots >= 0) {
1665:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1666:     PetscCalloc1(nroots, &neg);
1667:     if (nroots > pEnd-pStart) {
1668:       PetscCalloc1(nroots, &tmpOff);
1669:     } else {
1670:       tmpOff = &(*gsection)->atlasDof[-pStart];
1671:     }
1672:   }
1673:   /* Mark ghost points with negative dof */
1674:   for (p = pStart; p < pEnd; ++p) {
1675:     PetscInt value;

1677:     DMLabelGetValue(label, p, &value);
1678:     if (value != labelValue) continue;
1679:     PetscSectionGetDof(s, p, &dof);
1680:     PetscSectionSetDof(*gsection, p, dof);
1681:     PetscSectionGetConstraintDof(s, p, &cdof);
1682:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1683:     if (neg) neg[p] = -(dof+1);
1684:   }
1685:   PetscSectionSetUpBC(*gsection);
1686:   if (nroots >= 0) {
1687:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1688:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1689:     if (nroots > pEnd-pStart) {
1690:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1691:     }
1692:   }
1693:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1694:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1695:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1696:     (*gsection)->atlasOff[p] = off;
1697:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1698:   }
1699:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1700:   globalOff -= off;
1701:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1702:     (*gsection)->atlasOff[p] += globalOff;
1703:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1704:   }
1705:   /* Put in negative offsets for ghost points */
1706:   if (nroots >= 0) {
1707:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1708:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1709:     if (nroots > pEnd-pStart) {
1710:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1711:     }
1712:   }
1713:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1714:   PetscFree(neg);
1715:   return(0);
1716: }

1718: typedef struct _n_PetscSectionSym_Label
1719: {
1720:   DMLabel           label;
1721:   PetscCopyMode     *modes;
1722:   PetscInt          *sizes;
1723:   const PetscInt    ***perms;
1724:   const PetscScalar ***rots;
1725:   PetscInt          (*minMaxOrients)[2];
1726:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1727: } PetscSectionSym_Label;

1729: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1730: {
1731:   PetscInt              i, j;
1732:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1733:   PetscErrorCode        ierr;

1736:   for (i = 0; i <= sl->numStrata; i++) {
1737:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1738:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1739:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1740:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1741:       }
1742:       if (sl->perms[i]) {
1743:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1745:         PetscFree(perms);
1746:       }
1747:       if (sl->rots[i]) {
1748:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1750:         PetscFree(rots);
1751:       }
1752:     }
1753:   }
1754:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1755:   DMLabelDestroy(&sl->label);
1756:   sl->numStrata = 0;
1757:   return(0);
1758: }

1760: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1761: {

1765:   PetscSectionSymLabelReset(sym);
1766:   PetscFree(sym->data);
1767:   return(0);
1768: }

1770: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1771: {
1772:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1773:   PetscBool             isAscii;
1774:   DMLabel               label = sl->label;
1775:   const char           *name;
1776:   PetscErrorCode        ierr;

1779:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1780:   if (isAscii) {
1781:     PetscInt          i, j, k;
1782:     PetscViewerFormat format;

1784:     PetscViewerGetFormat(viewer,&format);
1785:     if (label) {
1786:       PetscViewerGetFormat(viewer,&format);
1787:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1788:         PetscViewerASCIIPushTab(viewer);
1789:         DMLabelView(label, viewer);
1790:         PetscViewerASCIIPopTab(viewer);
1791:       } else {
1792:         PetscObjectGetName((PetscObject) sl->label, &name);
1793:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
1794:       }
1795:     } else {
1796:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1797:     }
1798:     PetscViewerASCIIPushTab(viewer);
1799:     for (i = 0; i <= sl->numStrata; i++) {
1800:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1802:       if (!(sl->perms[i] || sl->rots[i])) {
1803:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1804:       } else {
1805:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1806:         PetscViewerASCIIPushTab(viewer);
1807:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1808:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1809:           PetscViewerASCIIPushTab(viewer);
1810:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1811:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1812:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1813:             } else {
1814:               PetscInt tab;

1816:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1817:               PetscViewerASCIIPushTab(viewer);
1818:               PetscViewerASCIIGetTab(viewer,&tab);
1819:               if (sl->perms[i] && sl->perms[i][j]) {
1820:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1821:                 PetscViewerASCIISetTab(viewer,0);
1822:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1823:                 PetscViewerASCIIPrintf(viewer,"\n");
1824:                 PetscViewerASCIISetTab(viewer,tab);
1825:               }
1826:               if (sl->rots[i] && sl->rots[i][j]) {
1827:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
1828:                 PetscViewerASCIISetTab(viewer,0);
1829: #if defined(PETSC_USE_COMPLEX)
1830:                 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]));}
1831: #else
1832:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1833: #endif
1834:                 PetscViewerASCIIPrintf(viewer,"\n");
1835:                 PetscViewerASCIISetTab(viewer,tab);
1836:               }
1837:               PetscViewerASCIIPopTab(viewer);
1838:             }
1839:           }
1840:           PetscViewerASCIIPopTab(viewer);
1841:         }
1842:         PetscViewerASCIIPopTab(viewer);
1843:       }
1844:     }
1845:     PetscViewerASCIIPopTab(viewer);
1846:   }
1847:   return(0);
1848: }

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

1853:   Logically collective on sym

1855:   Input parameters:
1856: + sym - the section symmetries
1857: - label - the DMLabel describing the types of points

1859:   Level: developer:

1861: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1862: @*/
1863: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1864: {
1865:   PetscSectionSym_Label *sl;
1866:   PetscErrorCode        ierr;

1870:   sl = (PetscSectionSym_Label *) sym->data;
1871:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1872:   if (label) {
1873:     sl->label = label;
1874:     PetscObjectReference((PetscObject) label);
1875:     DMLabelGetNumValues(label,&sl->numStrata);
1876:     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);
1877:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1878:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1879:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1880:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1881:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1882:   }
1883:   return(0);
1884: }

1886: /*@C
1887:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

1889:   Logically collective on PetscSectionSym

1891:   InputParameters:
1892: + sys       - the section symmetries
1893: . stratum   - the stratum value in the label that we are assigning symmetries for
1894: . size      - the number of dofs for points in the stratum of the label
1895: . minOrient - the smallest orientation for a point in this stratum
1896: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1897: . mode      - how sym should copy the perms and rots arrays
1898: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
1899: + 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

1901:   Level: developer

1903: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1904: @*/
1905: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1906: {
1907:   PetscSectionSym_Label *sl;
1908:   const char            *name;
1909:   PetscInt               i, j, k;
1910:   PetscErrorCode         ierr;

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

1919:     if (stratum == value) break;
1920:   }
1921:   PetscObjectGetName((PetscObject) sl->label, &name);
1922:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
1923:   sl->sizes[i] = size;
1924:   sl->modes[i] = mode;
1925:   sl->minMaxOrients[i][0] = minOrient;
1926:   sl->minMaxOrients[i][1] = maxOrient;
1927:   if (mode == PETSC_COPY_VALUES) {
1928:     if (perms) {
1929:       PetscInt    **ownPerms;

1931:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
1932:       for (j = 0; j < maxOrient-minOrient; j++) {
1933:         if (perms[j]) {
1934:           PetscMalloc1(size,&ownPerms[j]);
1935:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1936:         }
1937:       }
1938:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1939:     }
1940:     if (rots) {
1941:       PetscScalar **ownRots;

1943:       PetscCalloc1(maxOrient - minOrient,&ownRots);
1944:       for (j = 0; j < maxOrient-minOrient; j++) {
1945:         if (rots[j]) {
1946:           PetscMalloc1(size,&ownRots[j]);
1947:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1948:         }
1949:       }
1950:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1951:     }
1952:   } else {
1953:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1954:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
1955:   }
1956:   return(0);
1957: }

1959: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1960: {
1961:   PetscInt              i, j, numStrata;
1962:   PetscSectionSym_Label *sl;
1963:   DMLabel               label;
1964:   PetscErrorCode        ierr;

1967:   sl = (PetscSectionSym_Label *) sym->data;
1968:   numStrata = sl->numStrata;
1969:   label     = sl->label;
1970:   for (i = 0; i < numPoints; i++) {
1971:     PetscInt point = points[2*i];
1972:     PetscInt ornt  = points[2*i+1];

1974:     for (j = 0; j < numStrata; j++) {
1975:       if (label->validIS[j]) {
1976:         PetscInt k;

1978:         ISLocate(label->points[j],point,&k);
1979:         if (k >= 0) break;
1980:       } else {
1981:         PetscBool has;

1983:         PetscHSetIHas(label->ht[j], point, &has);
1984:         if (has) break;
1985:       }
1986:     }
1987:     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);
1988:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1989:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1990:   }
1991:   return(0);
1992: }

1994: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1995: {
1996:   PetscSectionSym_Label *sl;
1997:   PetscErrorCode        ierr;

2000:   PetscNewLog(sym,&sl);
2001:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2002:   sym->ops->view      = PetscSectionSymView_Label;
2003:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
2004:   sym->data           = (void *) sl;
2005:   return(0);
2006: }

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

2011:   Collective

2013:   Input Parameters:
2014: + comm - the MPI communicator for the new symmetry
2015: - label - the label defining the strata

2017:   Output Parameters:
2018: . sym - the section symmetries

2020:   Level: developer

2022: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2023: @*/
2024: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2025: {
2026:   PetscErrorCode        ierr;

2029:   DMInitializePackage();
2030:   PetscSectionSymCreate(comm,sym);
2031:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2032:   PetscSectionSymLabelSetLabel(*sym,label);
2033:   return(0);
2034: }