Actual source code: dmlabel.c

petsc-3.13.6 2020-09-29
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:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
 90:     ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
 91:     PetscFree(pointArray);
 92:   } else {
 93:     ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
 94:   }
 95:   PetscObjectSetName((PetscObject) is, "indices");
 96:   label->points[v]  = is;
 97:   label->validIS[v] = PETSC_TRUE;
 98:   PetscObjectStateIncrease((PetscObject) label);
 99:   return(0);
100: }

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

105:   Not collective

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

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

113:   Level: developer

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

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

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

132:   Not collective

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

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

141:   Level: developer

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

151:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
153:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
154:   if (label->points[v]) {
155:     ISGetIndices(label->points[v], &points);
156:     for (p = 0; p < label->stratumSizes[v]; ++p) {
157:       PetscHSetIAdd(label->ht[v], points[p]);
158:     }
159:     ISRestoreIndices(label->points[v],&points);
160:     ISDestroy(&(label->points[v]));
161:   }
162:   label->validIS[v] = PETSC_FALSE;
163:   return(0);
164: }

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

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

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

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

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

259: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
260: {
263:   DMLabelLookupStratum(label, value, index);
264:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
265:   return(0);
266: }

268: /*@
269:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

271:   Input Parameter:
272: + label - The DMLabel
273: - value - The stratum value

275:   Level: beginner

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

286:   DMLabelLookupAddStratum(label, value, &v);
287:   return(0);
288: }

290: /*@
291:   DMLabelAddStrata - Adds new stratum values in a DMLabel

293:   Not collective

295:   Input Parameter:
296: + label - The DMLabel
297: . numStrata - The number of stratum values
298: - stratumValues - The stratum values

300:   Level: beginner

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

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

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

354: /*@
355:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

357:   Not collective

359:   Input Parameter:
360: + label - The DMLabel
361: - valueIS - Index set with stratum values

363:   Level: beginner

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

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

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

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

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

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

414: /*@C
415:   DMLabelView - View the label

417:   Collective on viewer

419:   Input Parameters:
420: + label - The DMLabel
421: - viewer - The PetscViewer

423:   Level: intermediate

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

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

444: /*@
445:   DMLabelReset - Destroys internal data structures in a DMLabel

447:   Not collective

449:   Input Parameter:
450: . label - The DMLabel

452:   Level: beginner

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

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

480: /*@
481:   DMLabelDestroy - Destroys a DMLabel

483:   Collective on label

485:   Input Parameter:
486: . label - The DMLabel

488:   Level: beginner

490: .seealso: DMLabelReset(), DMLabelCreate()
491: @*/
492: PetscErrorCode DMLabelDestroy(DMLabel *label)
493: {

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

506: /*@
507:   DMLabelDuplicate - Duplicates a DMLabel

509:   Collective on label

511:   Input Parameter:
512: . label - The DMLabel

514:   Output Parameter:
515: . labelnew - location to put new vector

517:   Level: intermediate

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

529:   DMLabelMakeAllValid_Private(label);
530:   PetscObjectGetName((PetscObject) label, &name);
531:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

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

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

559:   Not collective

561:   Input Parameter:
562: . label  - The DMLabel

564:   Level: intermediate

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

575:   DMLabelMakeAllValid_Private(label);
576:   for (v = 0; v < label->numStrata; ++v) {
577:     const PetscInt *points;
578:     PetscInt       i;

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

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

595: /*@
596:   DMLabelCreateIndex - Create an index structure for membership determination

598:   Not collective

600:   Input Parameters:
601: + label  - The DMLabel
602: . pStart - The smallest point
603: - pEnd   - The largest point + 1

605:   Level: intermediate

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

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

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

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

638: /*@
639:   DMLabelDestroyIndex - Destroy the index structure

641:   Not collective

643:   Input Parameter:
644: . label - the DMLabel

646:   Level: intermediate

648: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
649: @*/
650: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
651: {

656:   label->pStart = -1;
657:   label->pEnd   = -1;
658:   PetscBTDestroy(&label->bt);
659:   return(0);
660: }

662: /*@
663:   DMLabelGetBounds - Return the smallest and largest point in the label

665:   Not collective

667:   Input Parameter:
668: . label - the DMLabel

670:   Output Parameters:
671: + pStart - The smallest point
672: - pEnd   - The largest point + 1

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

676:   Level: intermediate

678: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
679: @*/
680: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
681: {

686:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
687:   if (pStart) {
689:     *pStart = label->pStart;
690:   }
691:   if (pEnd) {
693:     *pEnd = label->pEnd;
694:   }
695:   return(0);
696: }

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

701:   Not collective

703:   Input Parameters:
704: + label - the DMLabel
705: - value - the value

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

710:   Level: developer

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

722:   DMLabelLookupStratum(label, value, &v);
723:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
724:   return(0);
725: }

727: /*@
728:   DMLabelHasPoint - Determine whether a label assigns a value to a point

730:   Not collective

732:   Input Parameters:
733: + label - the DMLabel
734: - point - the point

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

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

741:   Level: developer

743: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
744: @*/
745: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
746: {

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

761: /*@
762:   DMLabelStratumHasPoint - Return true if the stratum contains a point

764:   Not collective

766:   Input Parameters:
767: + label - the DMLabel
768: . value - the stratum value
769: - point - the point

771:   Output Parameter:
772: . contains - true if the stratum contains the point

774:   Level: intermediate

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

786:   *contains = PETSC_FALSE;
787:   DMLabelLookupStratum(label, value, &v);
788:   if (v < 0) return(0);

790:   if (label->validIS[v]) {
791:     PetscInt i;

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

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

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

808:   Not collective

810:   Input parameter:
811: . label - a DMLabel object

813:   Output parameter:
814: . defaultValue - the default value

816:   Level: beginner

818: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
819: @*/
820: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
821: {
824:   *defaultValue = label->defaultValue;
825:   return(0);
826: }

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

832:   Not collective

834:   Input parameter:
835: . label - a DMLabel object

837:   Output parameter:
838: . defaultValue - the default value

840:   Level: beginner

842: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
843: @*/
844: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
845: {
848:   label->defaultValue = defaultValue;
849:   return(0);
850: }

852: /*@
853:   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())

855:   Not collective

857:   Input Parameters:
858: + label - the DMLabel
859: - point - the point

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

864:   Level: intermediate

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

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

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

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

899: /*@
900:   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.

902:   Not collective

904:   Input Parameters:
905: + label - the DMLabel
906: . point - the point
907: - value - The point value

909:   Level: intermediate

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

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

929: /*@
930:   DMLabelClearValue - Clear the value a label assigns to a point

932:   Not collective

934:   Input Parameters:
935: + label - the DMLabel
936: . point - the point
937: - value - The point value

939:   Level: intermediate

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

950:   /* Find label value */
951:   DMLabelLookupStratum(label, value, &v);
952:   if (v < 0) return(0);

954:   if (label->bt) {
955:     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);
956:     PetscBTClear(label->bt, point - label->pStart);
957:   }

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

965: /*@
966:   DMLabelInsertIS - Set all points in the IS to a value

968:   Not collective

970:   Input Parameters:
971: + label - the DMLabel
972: . is    - the point IS
973: - value - The point value

975:   Level: intermediate

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

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

1000: /*@
1001:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

1003:   Not collective

1005:   Input Parameter:
1006: . label - the DMLabel

1008:   Output Paramater:
1009: . numValues - the number of values

1011:   Level: intermediate

1013: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1014: @*/
1015: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1016: {
1020:   *numValues = label->numStrata;
1021:   return(0);
1022: }

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

1027:   Not collective

1029:   Input Parameter:
1030: . label - the DMLabel

1032:   Output Paramater:
1033: . is    - the value IS

1035:   Level: intermediate

1037: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1038: @*/
1039: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1040: {

1046:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1047:   return(0);
1048: }

1050: /*@
1051:   DMLabelHasStratum - Determine whether points exist with the given value

1053:   Not collective

1055:   Input Parameters:
1056: + label - the DMLabel
1057: - value - the stratum value

1059:   Output Paramater:
1060: . exists - Flag saying whether points exist

1062:   Level: intermediate

1064: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1065: @*/
1066: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1067: {
1068:   PetscInt       v;

1074:   DMLabelLookupStratum(label, value, &v);
1075:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1076:   return(0);
1077: }

1079: /*@
1080:   DMLabelGetStratumSize - Get the size of a stratum

1082:   Not collective

1084:   Input Parameters:
1085: + label - the DMLabel
1086: - value - the stratum value

1088:   Output Paramater:
1089: . size - The number of points in the stratum

1091:   Level: intermediate

1093: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1094: @*/
1095: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1096: {
1097:   PetscInt       v;

1103:   *size = 0;
1104:   DMLabelLookupStratum(label, value, &v);
1105:   if (v < 0) return(0);
1106:   DMLabelMakeValid_Private(label, v);
1107:   *size = label->stratumSizes[v];
1108:   return(0);
1109: }

1111: /*@
1112:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1114:   Not collective

1116:   Input Parameters:
1117: + label - the DMLabel
1118: - value - the stratum value

1120:   Output Paramaters:
1121: + start - the smallest point in the stratum
1122: - end - the largest point in the stratum

1124:   Level: intermediate

1126: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1127: @*/
1128: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1129: {
1130:   PetscInt       v, min, max;

1137:   DMLabelLookupStratum(label, value, &v);
1138:   if (v < 0) return(0);
1139:   DMLabelMakeValid_Private(label, v);
1140:   if (label->stratumSizes[v] <= 0) return(0);
1141:   ISGetMinMax(label->points[v], &min, &max);
1142:   if (start) *start = min;
1143:   if (end)   *end   = max+1;
1144:   return(0);
1145: }

1147: /*@
1148:   DMLabelGetStratumIS - Get an IS with the stratum points

1150:   Not collective

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

1156:   Output Paramater:
1157: . points - The stratum points

1159:   Level: intermediate

1161:   Notes:
1162:   The output IS should be destroyed when no longer needed.
1163:   Returns NULL if the stratum is empty.

1165: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1166: @*/
1167: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1168: {
1169:   PetscInt       v;

1175:   *points = NULL;
1176:   DMLabelLookupStratum(label, value, &v);
1177:   if (v < 0) return(0);
1178:   DMLabelMakeValid_Private(label, v);
1179:   PetscObjectReference((PetscObject) label->points[v]);
1180:   *points = label->points[v];
1181:   return(0);
1182: }

1184: /*@
1185:   DMLabelSetStratumIS - Set the stratum points using an IS

1187:   Not collective

1189:   Input Parameters:
1190: + label - the DMLabel
1191: . value - the stratum value
1192: - points - The stratum points

1194:   Level: intermediate

1196: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1197: @*/
1198: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1199: {
1200:   PetscInt       v;

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

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

1223:       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);
1224:       PetscBTSet(label->bt, point - label->pStart);
1225:     }
1226:   }
1227:   return(0);
1228: }

1230: /*@
1231:   DMLabelClearStratum - Remove a stratum

1233:   Not collective

1235:   Input Parameters:
1236: + label - the DMLabel
1237: - value - the stratum value

1239:   Level: intermediate

1241: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1242: @*/
1243: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1244: {
1245:   PetscInt       v;

1250:   DMLabelLookupStratum(label, value, &v);
1251:   if (v < 0) return(0);
1252:   if (label->validIS[v]) {
1253:     if (label->bt) {
1254:       PetscInt       i;
1255:       const PetscInt *points;

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

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

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

1280:   Not collective

1282:   Input Parameters:
1283: + label  - The DMLabel
1284: . value  - The label value for all points
1285: . pStart - The first point
1286: - pEnd   - A point beyond all marked points

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

1290:   Level: intermediate

1292: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1293: @*/
1294: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1295: {
1296:   IS             pIS;

1300:   ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1301:   DMLabelSetStratumIS(label, value, pIS);
1302:   ISDestroy(&pIS);
1303:   return(0);
1304: }

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

1309:   Not collective

1311:   Input Parameters:
1312: + label - the DMLabel
1313: . start - the first point kept
1314: - end - one more than the last point kept

1316:   Level: intermediate

1318: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1319: @*/
1320: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1321: {
1322:   PetscInt       v;

1327:   DMLabelDestroyIndex(label);
1328:   DMLabelMakeAllValid_Private(label);
1329:   for (v = 0; v < label->numStrata; ++v) {
1330:     ISGeneralFilter(label->points[v], start, end);
1331:   }
1332:   DMLabelCreateIndex(label, start, end);
1333:   return(0);
1334: }

1336: /*@
1337:   DMLabelPermute - Create a new label with permuted points

1339:   Not collective

1341:   Input Parameters:
1342: + label - the DMLabel
1343: - permutation - the point permutation

1345:   Output Parameter:
1346: . labelnew - the new label containing the permuted points

1348:   Level: intermediate

1350: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1351: @*/
1352: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1353: {
1354:   const PetscInt *perm;
1355:   PetscInt        numValues, numPoints, v, q;
1356:   PetscErrorCode  ierr;

1361:   DMLabelMakeAllValid_Private(label);
1362:   DMLabelDuplicate(label, labelNew);
1363:   DMLabelGetNumValues(*labelNew, &numValues);
1364:   ISGetLocalSize(permutation, &numPoints);
1365:   ISGetIndices(permutation, &perm);
1366:   for (v = 0; v < numValues; ++v) {
1367:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1368:     const PetscInt *points;
1369:     PetscInt *pointsNew;

1371:     ISGetIndices((*labelNew)->points[v],&points);
1372:     PetscMalloc1(size,&pointsNew);
1373:     for (q = 0; q < size; ++q) {
1374:       const PetscInt point = points[q];

1376:       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);
1377:       pointsNew[q] = perm[point];
1378:     }
1379:     ISRestoreIndices((*labelNew)->points[v],&points);
1380:     PetscSortInt(size, pointsNew);
1381:     ISDestroy(&((*labelNew)->points[v]));
1382:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1383:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1384:       PetscFree(pointsNew);
1385:     } else {
1386:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1387:     }
1388:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1389:   }
1390:   ISRestoreIndices(permutation, &perm);
1391:   if (label->bt) {
1392:     PetscBTDestroy(&label->bt);
1393:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1394:   }
1395:   return(0);
1396: }

1398: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1399: {
1400:   MPI_Comm       comm;
1401:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1402:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1403:   PetscSection   rootSection;
1404:   PetscSF        labelSF;

1408:   if (label) {DMLabelMakeAllValid_Private(label);}
1409:   PetscObjectGetComm((PetscObject)sf, &comm);
1410:   /* Build a section of stratum values per point, generate the according SF
1411:      and distribute point-wise stratum values to leaves. */
1412:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1413:   PetscSectionCreate(comm, &rootSection);
1414:   PetscSectionSetChart(rootSection, 0, nroots);
1415:   if (label) {
1416:     for (s = 0; s < label->numStrata; ++s) {
1417:       const PetscInt *points;

1419:       ISGetIndices(label->points[s], &points);
1420:       for (l = 0; l < label->stratumSizes[s]; l++) {
1421:         PetscSectionGetDof(rootSection, points[l], &dof);
1422:         PetscSectionSetDof(rootSection, points[l], dof+1);
1423:       }
1424:       ISRestoreIndices(label->points[s], &points);
1425:     }
1426:   }
1427:   PetscSectionSetUp(rootSection);
1428:   /* Create a point-wise array of stratum values */
1429:   PetscSectionGetStorageSize(rootSection, &size);
1430:   PetscMalloc1(size, &rootStrata);
1431:   PetscCalloc1(nroots, &rootIdx);
1432:   if (label) {
1433:     for (s = 0; s < label->numStrata; ++s) {
1434:       const PetscInt *points;

1436:       ISGetIndices(label->points[s], &points);
1437:       for (l = 0; l < label->stratumSizes[s]; l++) {
1438:         const PetscInt p = points[l];
1439:         PetscSectionGetOffset(rootSection, p, &offset);
1440:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1441:       }
1442:       ISRestoreIndices(label->points[s], &points);
1443:     }
1444:   }
1445:   /* Build SF that maps label points to remote processes */
1446:   PetscSectionCreate(comm, leafSection);
1447:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1448:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1449:   PetscFree(remoteOffsets);
1450:   /* Send the strata for each point over the derived SF */
1451:   PetscSectionGetStorageSize(*leafSection, &size);
1452:   PetscMalloc1(size, leafStrata);
1453:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1454:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1455:   /* Clean up */
1456:   PetscFree(rootStrata);
1457:   PetscFree(rootIdx);
1458:   PetscSectionDestroy(&rootSection);
1459:   PetscSFDestroy(&labelSF);
1460:   return(0);
1461: }

1463: /*@
1464:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1466:   Collective on sf

1468:   Input Parameters:
1469: + label - the DMLabel
1470: - sf    - the map from old to new distribution

1472:   Output Parameter:
1473: . labelnew - the new redistributed label

1475:   Level: intermediate

1477: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1478: @*/
1479: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1480: {
1481:   MPI_Comm       comm;
1482:   PetscSection   leafSection;
1483:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1484:   PetscInt      *leafStrata, *strataIdx;
1485:   PetscInt     **points;
1486:   const char    *lname = NULL;
1487:   char          *name;
1488:   PetscInt       nameSize;
1489:   PetscHSetI     stratumHash;
1490:   size_t         len = 0;
1491:   PetscMPIInt    rank;

1496:   if (label) {
1498:     DMLabelMakeAllValid_Private(label);
1499:   }
1500:   PetscObjectGetComm((PetscObject)sf, &comm);
1501:   MPI_Comm_rank(comm, &rank);
1502:   /* Bcast name */
1503:   if (!rank) {
1504:     PetscObjectGetName((PetscObject) label, &lname);
1505:     PetscStrlen(lname, &len);
1506:   }
1507:   nameSize = len;
1508:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1509:   PetscMalloc1(nameSize+1, &name);
1510:   if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1511:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1512:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1513:   PetscFree(name);
1514:   /* Bcast defaultValue */
1515:   if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1516:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1517:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1518:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1519:   /* Determine received stratum values and initialise new label*/
1520:   PetscHSetICreate(&stratumHash);
1521:   PetscSectionGetStorageSize(leafSection, &size);
1522:   for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1523:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1524:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1525:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1526:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1527:   /* Turn leafStrata into indices rather than stratum values */
1528:   offset = 0;
1529:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1530:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1531:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1532:     PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1533:   }
1534:   for (p = 0; p < size; ++p) {
1535:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1536:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1537:     }
1538:   }
1539:   /* Rebuild the point strata on the receiver */
1540:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1541:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1542:   for (p=pStart; p<pEnd; p++) {
1543:     PetscSectionGetDof(leafSection, p, &dof);
1544:     PetscSectionGetOffset(leafSection, p, &offset);
1545:     for (s=0; s<dof; s++) {
1546:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1547:     }
1548:   }
1549:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1550:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1551:   PetscMalloc1((*labelNew)->numStrata,&points);
1552:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1553:     PetscHSetICreate(&(*labelNew)->ht[s]);
1554:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1555:   }
1556:   /* Insert points into new strata */
1557:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1558:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1559:   for (p=pStart; p<pEnd; p++) {
1560:     PetscSectionGetDof(leafSection, p, &dof);
1561:     PetscSectionGetOffset(leafSection, p, &offset);
1562:     for (s=0; s<dof; s++) {
1563:       stratum = leafStrata[offset+s];
1564:       points[stratum][strataIdx[stratum]++] = p;
1565:     }
1566:   }
1567:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1568:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1569:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1570:   }
1571:   PetscFree(points);
1572:   PetscHSetIDestroy(&stratumHash);
1573:   PetscFree(leafStrata);
1574:   PetscFree(strataIdx);
1575:   PetscSectionDestroy(&leafSection);
1576:   return(0);
1577: }

1579: /*@
1580:   DMLabelGather - Gather all label values from leafs into roots

1582:   Collective on sf

1584:   Input Parameters:
1585: + label - the DMLabel
1586: - sf - the Star Forest point communication map

1588:   Output Parameters:
1589: . labelNew - the new DMLabel with localised leaf values

1591:   Level: developer

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

1595: .seealso: DMLabelDistribute()
1596: @*/
1597: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1598: {
1599:   MPI_Comm       comm;
1600:   PetscSection   rootSection;
1601:   PetscSF        sfLabel;
1602:   PetscSFNode   *rootPoints, *leafPoints;
1603:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1604:   const PetscInt *rootDegree, *ilocal;
1605:   PetscInt       *rootStrata;
1606:   const char    *lname;
1607:   char          *name;
1608:   PetscInt       nameSize;
1609:   size_t         len = 0;
1610:   PetscMPIInt    rank, size;

1616:   PetscObjectGetComm((PetscObject)sf, &comm);
1617:   MPI_Comm_rank(comm, &rank);
1618:   MPI_Comm_size(comm, &size);
1619:   /* Bcast name */
1620:   if (!rank) {
1621:     PetscObjectGetName((PetscObject) label, &lname);
1622:     PetscStrlen(lname, &len);
1623:   }
1624:   nameSize = len;
1625:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1626:   PetscMalloc1(nameSize+1, &name);
1627:   if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1628:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1629:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1630:   PetscFree(name);
1631:   /* Gather rank/index pairs of leaves into local roots to build
1632:      an inverse, multi-rooted SF. Note that this ignores local leaf
1633:      indexing due to the use of the multiSF in PetscSFGather. */
1634:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1635:   PetscMalloc1(nroots, &leafPoints);
1636:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1637:   for (p = 0; p < nleaves; p++) {
1638:     PetscInt ilp = ilocal ? ilocal[p] : p;

1640:     leafPoints[ilp].index = ilp;
1641:     leafPoints[ilp].rank  = rank;
1642:   }
1643:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1644:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1645:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1646:   PetscMalloc1(nmultiroots, &rootPoints);
1647:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1648:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1649:   PetscSFCreate(comm,& sfLabel);
1650:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1651:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1652:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1653:   /* Rebuild the point strata on the receiver */
1654:   for (p = 0, idx = 0; p < nroots; p++) {
1655:     for (d = 0; d < rootDegree[p]; d++) {
1656:       PetscSectionGetDof(rootSection, idx+d, &dof);
1657:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1658:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1659:     }
1660:     idx += rootDegree[p];
1661:   }
1662:   PetscFree(leafPoints);
1663:   PetscFree(rootStrata);
1664:   PetscSectionDestroy(&rootSection);
1665:   PetscSFDestroy(&sfLabel);
1666:   return(0);
1667: }

1669: /*@
1670:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1672:   Not collective

1674:   Input Parameter:
1675: . label - the DMLabel

1677:   Output Parameters:
1678: + section - the section giving offsets for each stratum
1679: - is - An IS containing all the label points

1681:   Level: developer

1683: .seealso: DMLabelDistribute()
1684: @*/
1685: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1686: {
1687:   IS              vIS;
1688:   const PetscInt *values;
1689:   PetscInt       *points;
1690:   PetscInt        nV, vS = 0, vE = 0, v, N;
1691:   PetscErrorCode  ierr;

1695:   DMLabelGetNumValues(label, &nV);
1696:   DMLabelGetValueIS(label, &vIS);
1697:   ISGetIndices(vIS, &values);
1698:   if (nV) {vS = values[0]; vE = values[0]+1;}
1699:   for (v = 1; v < nV; ++v) {
1700:     vS = PetscMin(vS, values[v]);
1701:     vE = PetscMax(vE, values[v]+1);
1702:   }
1703:   PetscSectionCreate(PETSC_COMM_SELF, section);
1704:   PetscSectionSetChart(*section, vS, vE);
1705:   for (v = 0; v < nV; ++v) {
1706:     PetscInt n;

1708:     DMLabelGetStratumSize(label, values[v], &n);
1709:     PetscSectionSetDof(*section, values[v], n);
1710:   }
1711:   PetscSectionSetUp(*section);
1712:   PetscSectionGetStorageSize(*section, &N);
1713:   PetscMalloc1(N, &points);
1714:   for (v = 0; v < nV; ++v) {
1715:     IS              is;
1716:     const PetscInt *spoints;
1717:     PetscInt        dof, off, p;

1719:     PetscSectionGetDof(*section, values[v], &dof);
1720:     PetscSectionGetOffset(*section, values[v], &off);
1721:     DMLabelGetStratumIS(label, values[v], &is);
1722:     ISGetIndices(is, &spoints);
1723:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1724:     ISRestoreIndices(is, &spoints);
1725:     ISDestroy(&is);
1726:   }
1727:   ISRestoreIndices(vIS, &values);
1728:   ISDestroy(&vIS);
1729:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1730:   return(0);
1731: }

1733: /*@
1734:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1735:   the local section and an SF describing the section point overlap.

1737:   Collective on sf

1739:   Input Parameters:
1740:   + s - The PetscSection for the local field layout
1741:   . sf - The SF describing parallel layout of the section points
1742:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1743:   . label - The label specifying the points
1744:   - labelValue - The label stratum specifying the points

1746:   Output Parameter:
1747:   . gsection - The PetscSection for the global field layout

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

1751:   Level: developer

1753: .seealso: PetscSectionCreate()
1754: @*/
1755: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1756: {
1757:   PetscInt      *neg = NULL, *tmpOff = NULL;
1758:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1765:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1766:   PetscSectionGetChart(s, &pStart, &pEnd);
1767:   PetscSectionSetChart(*gsection, pStart, pEnd);
1768:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1769:   if (nroots >= 0) {
1770:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1771:     PetscCalloc1(nroots, &neg);
1772:     if (nroots > pEnd-pStart) {
1773:       PetscCalloc1(nroots, &tmpOff);
1774:     } else {
1775:       tmpOff = &(*gsection)->atlasDof[-pStart];
1776:     }
1777:   }
1778:   /* Mark ghost points with negative dof */
1779:   for (p = pStart; p < pEnd; ++p) {
1780:     PetscInt value;

1782:     DMLabelGetValue(label, p, &value);
1783:     if (value != labelValue) continue;
1784:     PetscSectionGetDof(s, p, &dof);
1785:     PetscSectionSetDof(*gsection, p, dof);
1786:     PetscSectionGetConstraintDof(s, p, &cdof);
1787:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1788:     if (neg) neg[p] = -(dof+1);
1789:   }
1790:   PetscSectionSetUpBC(*gsection);
1791:   if (nroots >= 0) {
1792:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1793:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1794:     if (nroots > pEnd-pStart) {
1795:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1796:     }
1797:   }
1798:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1799:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1800:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1801:     (*gsection)->atlasOff[p] = off;
1802:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1803:   }
1804:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1805:   globalOff -= off;
1806:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1807:     (*gsection)->atlasOff[p] += globalOff;
1808:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1809:   }
1810:   /* Put in negative offsets for ghost points */
1811:   if (nroots >= 0) {
1812:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1813:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1814:     if (nroots > pEnd-pStart) {
1815:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1816:     }
1817:   }
1818:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1819:   PetscFree(neg);
1820:   return(0);
1821: }

1823: typedef struct _n_PetscSectionSym_Label
1824: {
1825:   DMLabel           label;
1826:   PetscCopyMode     *modes;
1827:   PetscInt          *sizes;
1828:   const PetscInt    ***perms;
1829:   const PetscScalar ***rots;
1830:   PetscInt          (*minMaxOrients)[2];
1831:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1832: } PetscSectionSym_Label;

1834: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1835: {
1836:   PetscInt              i, j;
1837:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1838:   PetscErrorCode        ierr;

1841:   for (i = 0; i <= sl->numStrata; i++) {
1842:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1843:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1844:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1845:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1846:       }
1847:       if (sl->perms[i]) {
1848:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1850:         PetscFree(perms);
1851:       }
1852:       if (sl->rots[i]) {
1853:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1855:         PetscFree(rots);
1856:       }
1857:     }
1858:   }
1859:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1860:   DMLabelDestroy(&sl->label);
1861:   sl->numStrata = 0;
1862:   return(0);
1863: }

1865: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1866: {

1870:   PetscSectionSymLabelReset(sym);
1871:   PetscFree(sym->data);
1872:   return(0);
1873: }

1875: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1876: {
1877:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1878:   PetscBool             isAscii;
1879:   DMLabel               label = sl->label;
1880:   const char           *name;
1881:   PetscErrorCode        ierr;

1884:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1885:   if (isAscii) {
1886:     PetscInt          i, j, k;
1887:     PetscViewerFormat format;

1889:     PetscViewerGetFormat(viewer,&format);
1890:     if (label) {
1891:       PetscViewerGetFormat(viewer,&format);
1892:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1893:         PetscViewerASCIIPushTab(viewer);
1894:         DMLabelView(label, viewer);
1895:         PetscViewerASCIIPopTab(viewer);
1896:       } else {
1897:         PetscObjectGetName((PetscObject) sl->label, &name);
1898:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
1899:       }
1900:     } else {
1901:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1902:     }
1903:     PetscViewerASCIIPushTab(viewer);
1904:     for (i = 0; i <= sl->numStrata; i++) {
1905:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1907:       if (!(sl->perms[i] || sl->rots[i])) {
1908:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1909:       } else {
1910:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1911:         PetscViewerASCIIPushTab(viewer);
1912:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1913:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1914:           PetscViewerASCIIPushTab(viewer);
1915:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1916:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1917:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1918:             } else {
1919:               PetscInt tab;

1921:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1922:               PetscViewerASCIIPushTab(viewer);
1923:               PetscViewerASCIIGetTab(viewer,&tab);
1924:               if (sl->perms[i] && sl->perms[i][j]) {
1925:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1926:                 PetscViewerASCIISetTab(viewer,0);
1927:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1928:                 PetscViewerASCIIPrintf(viewer,"\n");
1929:                 PetscViewerASCIISetTab(viewer,tab);
1930:               }
1931:               if (sl->rots[i] && sl->rots[i][j]) {
1932:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
1933:                 PetscViewerASCIISetTab(viewer,0);
1934: #if defined(PETSC_USE_COMPLEX)
1935:                 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]));}
1936: #else
1937:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1938: #endif
1939:                 PetscViewerASCIIPrintf(viewer,"\n");
1940:                 PetscViewerASCIISetTab(viewer,tab);
1941:               }
1942:               PetscViewerASCIIPopTab(viewer);
1943:             }
1944:           }
1945:           PetscViewerASCIIPopTab(viewer);
1946:         }
1947:         PetscViewerASCIIPopTab(viewer);
1948:       }
1949:     }
1950:     PetscViewerASCIIPopTab(viewer);
1951:   }
1952:   return(0);
1953: }

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

1958:   Logically collective on sym

1960:   Input parameters:
1961: + sym - the section symmetries
1962: - label - the DMLabel describing the types of points

1964:   Level: developer:

1966: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1967: @*/
1968: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1969: {
1970:   PetscSectionSym_Label *sl;
1971:   PetscErrorCode        ierr;

1975:   sl = (PetscSectionSym_Label *) sym->data;
1976:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1977:   if (label) {
1978:     sl->label = label;
1979:     PetscObjectReference((PetscObject) label);
1980:     DMLabelGetNumValues(label,&sl->numStrata);
1981:     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);
1982:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1983:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1984:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1985:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1986:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1987:   }
1988:   return(0);
1989: }

1991: /*@C
1992:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

1994:   Logically collective on sym

1996:   InputParameters:
1997: + sym       - the section symmetries
1998: . stratum   - the stratum value in the label that we are assigning symmetries for
1999: . size      - the number of dofs for points in the stratum of the label
2000: . minOrient - the smallest orientation for a point in this stratum
2001: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2002: . mode      - how sym should copy the perms and rots arrays
2003: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2004: - 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

2006:   Level: developer

2008: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2009: @*/
2010: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2011: {
2012:   PetscSectionSym_Label *sl;
2013:   const char            *name;
2014:   PetscInt               i, j, k;
2015:   PetscErrorCode         ierr;

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

2024:     if (stratum == value) break;
2025:   }
2026:   PetscObjectGetName((PetscObject) sl->label, &name);
2027:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2028:   sl->sizes[i] = size;
2029:   sl->modes[i] = mode;
2030:   sl->minMaxOrients[i][0] = minOrient;
2031:   sl->minMaxOrients[i][1] = maxOrient;
2032:   if (mode == PETSC_COPY_VALUES) {
2033:     if (perms) {
2034:       PetscInt    **ownPerms;

2036:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
2037:       for (j = 0; j < maxOrient-minOrient; j++) {
2038:         if (perms[j]) {
2039:           PetscMalloc1(size,&ownPerms[j]);
2040:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2041:         }
2042:       }
2043:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2044:     }
2045:     if (rots) {
2046:       PetscScalar **ownRots;

2048:       PetscCalloc1(maxOrient - minOrient,&ownRots);
2049:       for (j = 0; j < maxOrient-minOrient; j++) {
2050:         if (rots[j]) {
2051:           PetscMalloc1(size,&ownRots[j]);
2052:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2053:         }
2054:       }
2055:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2056:     }
2057:   } else {
2058:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2059:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2060:   }
2061:   return(0);
2062: }

2064: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2065: {
2066:   PetscInt              i, j, numStrata;
2067:   PetscSectionSym_Label *sl;
2068:   DMLabel               label;
2069:   PetscErrorCode        ierr;

2072:   sl = (PetscSectionSym_Label *) sym->data;
2073:   numStrata = sl->numStrata;
2074:   label     = sl->label;
2075:   for (i = 0; i < numPoints; i++) {
2076:     PetscInt point = points[2*i];
2077:     PetscInt ornt  = points[2*i+1];

2079:     for (j = 0; j < numStrata; j++) {
2080:       if (label->validIS[j]) {
2081:         PetscInt k;

2083:         ISLocate(label->points[j],point,&k);
2084:         if (k >= 0) break;
2085:       } else {
2086:         PetscBool has;

2088:         PetscHSetIHas(label->ht[j], point, &has);
2089:         if (has) break;
2090:       }
2091:     }
2092:     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);
2093:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2094:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2095:   }
2096:   return(0);
2097: }

2099: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2100: {
2101:   PetscSectionSym_Label *sl;
2102:   PetscErrorCode        ierr;

2105:   PetscNewLog(sym,&sl);
2106:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2107:   sym->ops->view      = PetscSectionSymView_Label;
2108:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
2109:   sym->data           = (void *) sl;
2110:   return(0);
2111: }

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

2116:   Collective

2118:   Input Parameters:
2119: + comm - the MPI communicator for the new symmetry
2120: - label - the label defining the strata

2122:   Output Parameters:
2123: . sym - the section symmetries

2125:   Level: developer

2127: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2128: @*/
2129: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2130: {
2131:   PetscErrorCode        ierr;

2134:   DMInitializePackage();
2135:   PetscSectionSymCreate(comm,sym);
2136:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2137:   PetscSectionSymLabelSetLabel(*sym,label);
2138:   return(0);
2139: }