Actual source code: dmlabel.c

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

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

 10:   Collective

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

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

 19:   Level: beginner

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

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

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

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

 49: /*
 50:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 52:   Not collective

 54:   Input parameter:
 55: + label - The DMLabel
 56: - v - The stratum value

 58:   Output parameter:
 59: . label - The DMLabel with stratum in sorted list format

 61:   Level: developer

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

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

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

100:   Not collective

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

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

108:   Level: developer

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

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

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

123:   Not collective

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

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

132:   Level: developer

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

141:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
143:   if (label->points[v]) {
144:     ISGetIndices(label->points[v], &points);
145:     for (p = 0; p < label->stratumSizes[v]; ++p) PetscHSetIAdd(label->ht[v], points[p]);
146:     ISRestoreIndices(label->points[v], &points);
147:     ISDestroy(&(label->points[v]));
148:   }
149:   label->validIS[v] = PETSC_FALSE;
150:   return 0;
151: }

153: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
154: {
155:   PetscInt v;

157:   for (v = 0; v < label->numStrata; v++) DMLabelMakeInvalid_Private(label, v);
158:   return 0;
159: }

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

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

169:   *index = -1;
170:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
171:     for (v = 0; v < label->numStrata; ++v)
172:       if (label->stratumValues[v] == value) {
173:         *index = v;
174:         break;
175:       }
176:   } else {
177:     PetscHMapIGet(label->hmap, value, index);
178:   }
179:   if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
180:     PetscInt len, loc = -1;
181:     PetscHMapIGetSize(label->hmap, &len);
183:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
184:       PetscHMapIGet(label->hmap, value, &loc);
185:     } else {
186:       for (v = 0; v < label->numStrata; ++v)
187:         if (label->stratumValues[v] == value) {
188:           loc = v;
189:           break;
190:         }
191:     }
193:   }
194:   return 0;
195: }

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

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

254: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
255: {
256:   DMLabelLookupStratum(label, value, index);
257:   if (*index < 0) DMLabelNewStratum(label, value, index);
258:   return 0;
259: }

261: static inline PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
262: {
263:   *size = 0;
264:   if (v < 0) return 0;
265:   if (label->validIS[v]) {
266:     *size = label->stratumSizes[v];
267:   } else {
268:     PetscHSetIGetSize(label->ht[v], size);
269:   }
270:   return 0;
271: }

273: /*@
274:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

276:   Input Parameters:
277: + label - The DMLabel
278: - value - The stratum value

280:   Level: beginner

282: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
283: @*/
284: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
285: {
286:   PetscInt v;

289:   DMLabelLookupAddStratum(label, value, &v);
290:   return 0;
291: }

293: /*@
294:   DMLabelAddStrata - Adds new stratum values in a DMLabel

296:   Not collective

298:   Input Parameters:
299: + label - The DMLabel
300: . numStrata - The number of stratum values
301: - stratumValues - The stratum values

303:   Level: beginner

305: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
306: @*/
307: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
308: {
309:   PetscInt *values, v;

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

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

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

356:   Not collective

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

362:   Level: beginner

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

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

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

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

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

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

407: /*@C
408:   DMLabelView - View the label

410:   Collective on viewer

412:   Input Parameters:
413: + label - The DMLabel
414: - viewer - The PetscViewer

416:   Level: intermediate

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

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

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

436:   Not collective

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

441:   Level: beginner

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

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

472: /*@
473:   DMLabelDestroy - Destroys a DMLabel

475:   Collective on label

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

480:   Level: beginner

482: .seealso: `DMLabelReset()`, `DMLabelCreate()`
483: @*/
484: PetscErrorCode DMLabelDestroy(DMLabel *label)
485: {
486:   if (!*label) return 0;
488:   if (--((PetscObject)(*label))->refct > 0) {
489:     *label = NULL;
490:     return 0;
491:   }
492:   DMLabelReset(*label);
493:   PetscHMapIDestroy(&(*label)->hmap);
494:   PetscHeaderDestroy(label);
495:   return 0;
496: }

498: /*@
499:   DMLabelDuplicate - Duplicates a DMLabel

501:   Collective on label

503:   Input Parameter:
504: . label - The DMLabel

506:   Output Parameter:
507: . labelnew - location to put new vector

509:   Level: intermediate

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

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

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

546: /*@C
547:   DMLabelCompare - Compare two DMLabel objects

549:   Collective on comm

551:   Input Parameters:
552: + comm - Comm over which to compare labels
553: . l0 - First DMLabel
554: - l1 - Second DMLabel

556:   Output Parameters
557: + equal   - (Optional) Flag whether the two labels are equal
558: - message - (Optional) Message describing the difference

560:   Level: intermediate

562:   Notes:
563:   The output flag equal is the same on all processes.
564:   If it is passed as NULL and difference is found, an error is thrown on all processes.
565:   Make sure to pass NULL on all processes.

567:   The output message is set independently on each rank.
568:   It is set to NULL if no difference was found on the current rank. It must be freed by user.
569:   If message is passed as NULL and difference is found, the difference description is printed to stderr in synchronized manner.
570:   Make sure to pass NULL on all processes.

572:   For the comparison, we ignore the order of stratum values, and strata with no points.

574:   The communicator needs to be specified because currently DMLabel can live on PETSC_COMM_SELF even if the underlying DM is parallel.

576:   Fortran Notes:
577:   This function is currently not available from Fortran.

579: .seealso: `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
580: @*/
581: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
582: {
583:   const char *name0, *name1;
584:   char        msg[PETSC_MAX_PATH_LEN] = "";
585:   PetscBool   eq;
586:   PetscMPIInt rank;

592:   MPI_Comm_rank(comm, &rank);
593:   PetscObjectGetName((PetscObject)l0, &name0);
594:   PetscObjectGetName((PetscObject)l1, &name1);
595:   {
596:     PetscInt v0, v1;

598:     DMLabelGetDefaultValue(l0, &v0);
599:     DMLabelGetDefaultValue(l1, &v1);
600:     eq = (PetscBool)(v0 == v1);
601:     if (!eq) PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %" PetscInt_FMT " != %" PetscInt_FMT " = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1);
602:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
603:     if (!eq) goto finish;
604:   }
605:   {
606:     IS is0, is1;

608:     DMLabelGetNonEmptyStratumValuesIS(l0, &is0);
609:     DMLabelGetNonEmptyStratumValuesIS(l1, &is1);
610:     ISEqual(is0, is1, &eq);
611:     ISDestroy(&is0);
612:     ISDestroy(&is1);
613:     if (!eq) PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1);
614:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
615:     if (!eq) goto finish;
616:   }
617:   {
618:     PetscInt i, nValues;

620:     DMLabelGetNumValues(l0, &nValues);
621:     for (i = 0; i < nValues; i++) {
622:       const PetscInt v = l0->stratumValues[i];
623:       PetscInt       n;
624:       IS             is0, is1;

626:       DMLabelGetStratumSize_Private(l0, i, &n);
627:       if (!n) continue;
628:       DMLabelGetStratumIS(l0, v, &is0);
629:       DMLabelGetStratumIS(l1, v, &is1);
630:       ISEqualUnsorted(is0, is1, &eq);
631:       ISDestroy(&is0);
632:       ISDestroy(&is1);
633:       if (!eq) {
634:         PetscSNPrintf(msg, sizeof(msg), "Stratum #%" PetscInt_FMT " with value %" PetscInt_FMT " contains different points in DMLabel l0 \"%s\" and DMLabel l1 \"%s\"", i, v, name0, name1);
635:         break;
636:       }
637:     }
638:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
639:   }
640: finish:
641:   /* If message output arg not set, print to stderr */
642:   if (message) {
643:     *message = NULL;
644:     if (msg[0]) PetscStrallocpy(msg, message);
645:   } else {
646:     if (msg[0]) PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg);
647:     PetscSynchronizedFlush(comm, PETSC_STDERR);
648:   }
649:   /* If same output arg not ser and labels are not equal, throw error */
650:   if (equal) *equal = eq;
652:   return 0;
653: }

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

658:   Not collective

660:   Input Parameter:
661: . label  - The DMLabel

663:   Level: intermediate

665: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
666: @*/
667: PetscErrorCode DMLabelComputeIndex(DMLabel label)
668: {
669:   PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;

672:   DMLabelMakeAllValid_Private(label);
673:   for (v = 0; v < label->numStrata; ++v) {
674:     const PetscInt *points;
675:     PetscInt        i;

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

681:       pStart = PetscMin(point, pStart);
682:       pEnd   = PetscMax(point + 1, pEnd);
683:     }
684:     ISRestoreIndices(label->points[v], &points);
685:   }
686:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
687:   label->pEnd   = pEnd;
688:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
689:   return 0;
690: }

692: /*@
693:   DMLabelCreateIndex - Create an index structure for membership determination

695:   Not collective

697:   Input Parameters:
698: + label  - The DMLabel
699: . pStart - The smallest point
700: - pEnd   - The largest point + 1

702:   Level: intermediate

704: .seealso: `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
705: @*/
706: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
707: {
708:   PetscInt v;

711:   DMLabelDestroyIndex(label);
712:   DMLabelMakeAllValid_Private(label);
713:   label->pStart = pStart;
714:   label->pEnd   = pEnd;
715:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
716:   PetscBTCreate(pEnd - pStart, &label->bt);
717:   for (v = 0; v < label->numStrata; ++v) {
718:     const PetscInt *points;
719:     PetscInt        i;

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

726:       PetscBTSet(label->bt, point - pStart);
727:     }
728:     ISRestoreIndices(label->points[v], &points);
729:   }
730:   return 0;
731: }

733: /*@
734:   DMLabelDestroyIndex - Destroy the index structure

736:   Not collective

738:   Input Parameter:
739: . label - the DMLabel

741:   Level: intermediate

743: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
744: @*/
745: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
746: {
748:   label->pStart = -1;
749:   label->pEnd   = -1;
750:   PetscBTDestroy(&label->bt);
751:   return 0;
752: }

754: /*@
755:   DMLabelGetBounds - Return the smallest and largest point in the label

757:   Not collective

759:   Input Parameter:
760: . label - the DMLabel

762:   Output Parameters:
763: + pStart - The smallest point
764: - pEnd   - The largest point + 1

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

768:   Level: intermediate

770: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
771: @*/
772: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
773: {
775:   if ((label->pStart == -1) && (label->pEnd == -1)) DMLabelComputeIndex(label);
776:   if (pStart) {
778:     *pStart = label->pStart;
779:   }
780:   if (pEnd) {
782:     *pEnd = label->pEnd;
783:   }
784:   return 0;
785: }

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

790:   Not collective

792:   Input Parameters:
793: + label - the DMLabel
794: - value - the value

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

799:   Level: developer

801: .seealso: `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
802: @*/
803: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
804: {
805:   PetscInt v;

809:   DMLabelLookupStratum(label, value, &v);
810:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
811:   return 0;
812: }

814: /*@
815:   DMLabelHasPoint - Determine whether a label assigns a value to a point

817:   Not collective

819:   Input Parameters:
820: + label - the DMLabel
821: - point - the point

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

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

828:   Level: developer

830: .seealso: `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
831: @*/
832: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
833: {
837:   DMLabelMakeAllValid_Private(label);
838:   if (PetscDefined(USE_DEBUG)) {
841:   }
842:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
843:   return 0;
844: }

846: /*@
847:   DMLabelStratumHasPoint - Return true if the stratum contains a point

849:   Not collective

851:   Input Parameters:
852: + label - the DMLabel
853: . value - the stratum value
854: - point - the point

856:   Output Parameter:
857: . contains - true if the stratum contains the point

859:   Level: intermediate

861: .seealso: `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
862: @*/
863: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
864: {
868:   if (value == label->defaultValue) {
869:     PetscInt pointVal;

871:     DMLabelGetValue(label, point, &pointVal);
872:     *contains = (PetscBool)(pointVal == value);
873:   } else {
874:     PetscInt v;

876:     DMLabelLookupStratum(label, value, &v);
877:     if (v >= 0) {
878:       if (label->validIS[v]) {
879:         PetscInt i;

881:         ISLocate(label->points[v], point, &i);
882:         *contains = (PetscBool)(i >= 0);
883:       } else {
884:         PetscHSetIHas(label->ht[v], point, contains);
885:       }
886:     } else { // value is not present
887:       *contains = PETSC_FALSE;
888:     }
889:   }
890:   return 0;
891: }

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

897:   Not collective

899:   Input parameter:
900: . label - a DMLabel object

902:   Output parameter:
903: . defaultValue - the default value

905:   Level: beginner

907: .seealso: `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
908: @*/
909: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
910: {
912:   *defaultValue = label->defaultValue;
913:   return 0;
914: }

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

920:   Not collective

922:   Input parameter:
923: . label - a DMLabel object

925:   Output parameter:
926: . defaultValue - the default value

928:   Level: beginner

930: .seealso: `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
931: @*/
932: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
933: {
935:   label->defaultValue = defaultValue;
936:   return 0;
937: }

939: /*@
940:   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())

942:   Not collective

944:   Input Parameters:
945: + label - the DMLabel
946: - point - the point

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

951:   Note: a label may assign multiple values to a point.  No guarantees are made about which value is returned in that case.  Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.

953:   Level: intermediate

955: .seealso: `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
956: @*/
957: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
958: {
959:   PetscInt v;

964:   *value = label->defaultValue;
965:   for (v = 0; v < label->numStrata; ++v) {
966:     if (label->validIS[v]) {
967:       PetscInt i;

969:       ISLocate(label->points[v], point, &i);
970:       if (i >= 0) {
971:         *value = label->stratumValues[v];
972:         break;
973:       }
974:     } else {
975:       PetscBool has;

977:       PetscHSetIHas(label->ht[v], point, &has);
978:       if (has) {
979:         *value = label->stratumValues[v];
980:         break;
981:       }
982:     }
983:   }
984:   return 0;
985: }

987: /*@
988:   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.

990:   Not collective

992:   Input Parameters:
993: + label - the DMLabel
994: . point - the point
995: - value - The point value

997:   Level: intermediate

999: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1000: @*/
1001: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1002: {
1003:   PetscInt v;

1006:   /* Find label value, add new entry if needed */
1007:   if (value == label->defaultValue) return 0;
1008:   DMLabelLookupAddStratum(label, value, &v);
1009:   /* Set key */
1010:   DMLabelMakeInvalid_Private(label, v);
1011:   PetscHSetIAdd(label->ht[v], point);
1012:   return 0;
1013: }

1015: /*@
1016:   DMLabelClearValue - Clear the value a label assigns to a point

1018:   Not collective

1020:   Input Parameters:
1021: + label - the DMLabel
1022: . point - the point
1023: - value - The point value

1025:   Level: intermediate

1027: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1028: @*/
1029: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1030: {
1031:   PetscInt v;

1034:   /* Find label value */
1035:   DMLabelLookupStratum(label, value, &v);
1036:   if (v < 0) return 0;

1038:   if (label->bt) {
1040:     PetscBTClear(label->bt, point - label->pStart);
1041:   }

1043:   /* Delete key */
1044:   DMLabelMakeInvalid_Private(label, v);
1045:   PetscHSetIDel(label->ht[v], point);
1046:   return 0;
1047: }

1049: /*@
1050:   DMLabelInsertIS - Set all points in the IS to a value

1052:   Not collective

1054:   Input Parameters:
1055: + label - the DMLabel
1056: . is    - the point IS
1057: - value - The point value

1059:   Level: intermediate

1061: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1062: @*/
1063: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1064: {
1065:   PetscInt        v, n, p;
1066:   const PetscInt *points;

1070:   /* Find label value, add new entry if needed */
1071:   if (value == label->defaultValue) return 0;
1072:   DMLabelLookupAddStratum(label, value, &v);
1073:   /* Set keys */
1074:   DMLabelMakeInvalid_Private(label, v);
1075:   ISGetLocalSize(is, &n);
1076:   ISGetIndices(is, &points);
1077:   for (p = 0; p < n; ++p) PetscHSetIAdd(label->ht[v], points[p]);
1078:   ISRestoreIndices(is, &points);
1079:   return 0;
1080: }

1082: /*@
1083:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

1085:   Not collective

1087:   Input Parameter:
1088: . label - the DMLabel

1090:   Output Parameter:
1091: . numValues - the number of values

1093:   Level: intermediate

1095: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1096: @*/
1097: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1098: {
1101:   *numValues = label->numStrata;
1102:   return 0;
1103: }

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

1108:   Not collective

1110:   Input Parameter:
1111: . label - the DMLabel

1113:   Output Parameter:
1114: . is    - the value IS

1116:   Level: intermediate

1118:   Notes:
1119:   The output IS should be destroyed when no longer needed.
1120:   Strata which are allocated but empty [DMLabelGetStratumSize() yields 0] are counted.
1121:   If you need to count only nonempty strata, use DMLabelGetNonEmptyStratumValuesIS().

1123: .seealso: `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1124: @*/
1125: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1126: {
1129:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1130:   return 0;
1131: }

1133: /*@
1134:   DMLabelGetNonEmptyStratumValuesIS - Get an IS of all values that the DMlabel takes

1136:   Not collective

1138:   Input Parameter:
1139: . label - the DMLabel

1141:   Output Parameter:
1142: . is    - the value IS

1144:   Level: intermediate

1146:   Notes:
1147:   The output IS should be destroyed when no longer needed.
1148:   This is similar to DMLabelGetValueIS() but counts only nonempty strata.

1150: .seealso: `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1151: @*/
1152: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1153: {
1154:   PetscInt  i, j;
1155:   PetscInt *valuesArr;

1159:   PetscMalloc1(label->numStrata, &valuesArr);
1160:   for (i = 0, j = 0; i < label->numStrata; i++) {
1161:     PetscInt n;

1163:     DMLabelGetStratumSize_Private(label, i, &n);
1164:     if (n) valuesArr[j++] = label->stratumValues[i];
1165:   }
1166:   if (j == label->numStrata) {
1167:     ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1168:   } else {
1169:     ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values);
1170:   }
1171:   PetscFree(valuesArr);
1172:   return 0;
1173: }

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

1178:   Not collective

1180:   Input Parameters:
1181: + label - the DMLabel
1182: - value - the value

1184:   Output Parameter:
1185: . index - the index of value in the list of values

1187:   Level: intermediate

1189: .seealso: `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1190: @*/
1191: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1192: {
1193:   PetscInt v;

1197:   /* Do not assume they are sorted */
1198:   for (v = 0; v < label->numStrata; ++v)
1199:     if (label->stratumValues[v] == value) break;
1200:   if (v >= label->numStrata) *index = -1;
1201:   else *index = v;
1202:   return 0;
1203: }

1205: /*@
1206:   DMLabelHasStratum - Determine whether points exist with the given value

1208:   Not collective

1210:   Input Parameters:
1211: + label - the DMLabel
1212: - value - the stratum value

1214:   Output Parameter:
1215: . exists - Flag saying whether points exist

1217:   Level: intermediate

1219: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1220: @*/
1221: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1222: {
1223:   PetscInt v;

1227:   DMLabelLookupStratum(label, value, &v);
1228:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1229:   return 0;
1230: }

1232: /*@
1233:   DMLabelGetStratumSize - Get the size of a stratum

1235:   Not collective

1237:   Input Parameters:
1238: + label - the DMLabel
1239: - value - the stratum value

1241:   Output Parameter:
1242: . size - The number of points in the stratum

1244:   Level: intermediate

1246: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1247: @*/
1248: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1249: {
1250:   PetscInt v;

1254:   DMLabelLookupStratum(label, value, &v);
1255:   DMLabelGetStratumSize_Private(label, v, size);
1256:   return 0;
1257: }

1259: /*@
1260:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1262:   Not collective

1264:   Input Parameters:
1265: + label - the DMLabel
1266: - value - the stratum value

1268:   Output Parameters:
1269: + start - the smallest point in the stratum
1270: - end - the largest point in the stratum

1272:   Level: intermediate

1274: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1275: @*/
1276: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1277: {
1278:   PetscInt v, min, max;

1281:   if (start) {
1283:     *start = -1;
1284:   }
1285:   if (end) {
1287:     *end = -1;
1288:   }
1289:   DMLabelLookupStratum(label, value, &v);
1290:   if (v < 0) return 0;
1291:   DMLabelMakeValid_Private(label, v);
1292:   if (label->stratumSizes[v] <= 0) return 0;
1293:   ISGetMinMax(label->points[v], &min, &max);
1294:   if (start) *start = min;
1295:   if (end) *end = max + 1;
1296:   return 0;
1297: }

1299: /*@
1300:   DMLabelGetStratumIS - Get an IS with the stratum points

1302:   Not collective

1304:   Input Parameters:
1305: + label - the DMLabel
1306: - value - the stratum value

1308:   Output Parameter:
1309: . points - The stratum points

1311:   Level: intermediate

1313:   Notes:
1314:   The output IS should be destroyed when no longer needed.
1315:   Returns NULL if the stratum is empty.

1317: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1318: @*/
1319: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1320: {
1321:   PetscInt v;

1325:   *points = NULL;
1326:   DMLabelLookupStratum(label, value, &v);
1327:   if (v < 0) return 0;
1328:   DMLabelMakeValid_Private(label, v);
1329:   PetscObjectReference((PetscObject)label->points[v]);
1330:   *points = label->points[v];
1331:   return 0;
1332: }

1334: /*@
1335:   DMLabelSetStratumIS - Set the stratum points using an IS

1337:   Not collective

1339:   Input Parameters:
1340: + label - the DMLabel
1341: . value - the stratum value
1342: - points - The stratum points

1344:   Level: intermediate

1346: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1347: @*/
1348: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1349: {
1350:   PetscInt v;

1354:   DMLabelLookupAddStratum(label, value, &v);
1355:   if (is == label->points[v]) return 0;
1356:   DMLabelClearStratum(label, value);
1357:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1358:   PetscObjectReference((PetscObject)is);
1359:   ISDestroy(&(label->points[v]));
1360:   label->points[v]  = is;
1361:   label->validIS[v] = PETSC_TRUE;
1362:   PetscObjectStateIncrease((PetscObject)label);
1363:   if (label->bt) {
1364:     const PetscInt *points;
1365:     PetscInt        p;

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

1372:       PetscBTSet(label->bt, point - label->pStart);
1373:     }
1374:   }
1375:   return 0;
1376: }

1378: /*@
1379:   DMLabelClearStratum - Remove a stratum

1381:   Not collective

1383:   Input Parameters:
1384: + label - the DMLabel
1385: - value - the stratum value

1387:   Level: intermediate

1389: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1390: @*/
1391: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1392: {
1393:   PetscInt v;

1396:   DMLabelLookupStratum(label, value, &v);
1397:   if (v < 0) return 0;
1398:   if (label->validIS[v]) {
1399:     if (label->bt) {
1400:       PetscInt        i;
1401:       const PetscInt *points;

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

1408:         PetscBTClear(label->bt, point - label->pStart);
1409:       }
1410:       ISRestoreIndices(label->points[v], &points);
1411:     }
1412:     label->stratumSizes[v] = 0;
1413:     ISDestroy(&label->points[v]);
1414:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1415:     PetscObjectSetName((PetscObject)label->points[v], "indices");
1416:     PetscObjectStateIncrease((PetscObject)label);
1417:   } else {
1418:     PetscHSetIClear(label->ht[v]);
1419:   }
1420:   return 0;
1421: }

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

1426:   Not collective

1428:   Input Parameters:
1429: + label  - The DMLabel
1430: . value  - The label value for all points
1431: . pStart - The first point
1432: - pEnd   - A point beyond all marked points

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

1436:   Level: intermediate

1438: .seealso: `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1439: @*/
1440: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1441: {
1442:   IS pIS;

1444:   ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1445:   DMLabelSetStratumIS(label, value, pIS);
1446:   ISDestroy(&pIS);
1447:   return 0;
1448: }

1450: /*@
1451:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

1453:   Not collective

1455:   Input Parameters:
1456: + label  - The DMLabel
1457: . value  - The label value
1458: - p      - A point with this value

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

1463:   Level: intermediate

1465: .seealso: `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1466: @*/
1467: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1468: {
1469:   const PetscInt *indices;
1470:   PetscInt        v;

1474:   *index = -1;
1475:   DMLabelLookupStratum(label, value, &v);
1476:   if (v < 0) return 0;
1477:   DMLabelMakeValid_Private(label, v);
1478:   ISGetIndices(label->points[v], &indices);
1479:   PetscFindInt(p, label->stratumSizes[v], indices, index);
1480:   ISRestoreIndices(label->points[v], &indices);
1481:   return 0;
1482: }

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

1487:   Not collective

1489:   Input Parameters:
1490: + label - the DMLabel
1491: . start - the first point kept
1492: - end - one more than the last point kept

1494:   Level: intermediate

1496: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1497: @*/
1498: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1499: {
1500:   PetscInt v;

1503:   DMLabelDestroyIndex(label);
1504:   DMLabelMakeAllValid_Private(label);
1505:   for (v = 0; v < label->numStrata; ++v) ISGeneralFilter(label->points[v], start, end);
1506:   DMLabelCreateIndex(label, start, end);
1507:   return 0;
1508: }

1510: /*@
1511:   DMLabelPermute - Create a new label with permuted points

1513:   Not collective

1515:   Input Parameters:
1516: + label - the DMLabel
1517: - permutation - the point permutation

1519:   Output Parameter:
1520: . labelnew - the new label containing the permuted points

1522:   Level: intermediate

1524: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1525: @*/
1526: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1527: {
1528:   const PetscInt *perm;
1529:   PetscInt        numValues, numPoints, v, q;

1533:   DMLabelMakeAllValid_Private(label);
1534:   DMLabelDuplicate(label, labelNew);
1535:   DMLabelGetNumValues(*labelNew, &numValues);
1536:   ISGetLocalSize(permutation, &numPoints);
1537:   ISGetIndices(permutation, &perm);
1538:   for (v = 0; v < numValues; ++v) {
1539:     const PetscInt  size = (*labelNew)->stratumSizes[v];
1540:     const PetscInt *points;
1541:     PetscInt       *pointsNew;

1543:     ISGetIndices((*labelNew)->points[v], &points);
1544:     PetscMalloc1(size, &pointsNew);
1545:     for (q = 0; q < size; ++q) {
1546:       const PetscInt point = points[q];

1549:       pointsNew[q] = perm[point];
1550:     }
1551:     ISRestoreIndices((*labelNew)->points[v], &points);
1552:     PetscSortInt(size, pointsNew);
1553:     ISDestroy(&((*labelNew)->points[v]));
1554:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1555:       ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v]));
1556:       PetscFree(pointsNew);
1557:     } else {
1558:       ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v]));
1559:     }
1560:     PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices");
1561:   }
1562:   ISRestoreIndices(permutation, &perm);
1563:   if (label->bt) {
1564:     PetscBTDestroy(&label->bt);
1565:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1566:   }
1567:   return 0;
1568: }

1570: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1571: {
1572:   MPI_Comm     comm;
1573:   PetscInt     s, l, nroots, nleaves, offset, size;
1574:   PetscInt    *remoteOffsets, *rootStrata, *rootIdx;
1575:   PetscSection rootSection;
1576:   PetscSF      labelSF;

1578:   if (label) DMLabelMakeAllValid_Private(label);
1579:   PetscObjectGetComm((PetscObject)sf, &comm);
1580:   /* Build a section of stratum values per point, generate the according SF
1581:      and distribute point-wise stratum values to leaves. */
1582:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1583:   PetscSectionCreate(comm, &rootSection);
1584:   PetscSectionSetChart(rootSection, 0, nroots);
1585:   if (label) {
1586:     for (s = 0; s < label->numStrata; ++s) {
1587:       const PetscInt *points;

1589:       ISGetIndices(label->points[s], &points);
1590:       for (l = 0; l < label->stratumSizes[s]; l++) PetscSectionAddDof(rootSection, points[l], 1);
1591:       ISRestoreIndices(label->points[s], &points);
1592:     }
1593:   }
1594:   PetscSectionSetUp(rootSection);
1595:   /* Create a point-wise array of stratum values */
1596:   PetscSectionGetStorageSize(rootSection, &size);
1597:   PetscMalloc1(size, &rootStrata);
1598:   PetscCalloc1(nroots, &rootIdx);
1599:   if (label) {
1600:     for (s = 0; s < label->numStrata; ++s) {
1601:       const PetscInt *points;

1603:       ISGetIndices(label->points[s], &points);
1604:       for (l = 0; l < label->stratumSizes[s]; l++) {
1605:         const PetscInt p = points[l];
1606:         PetscSectionGetOffset(rootSection, p, &offset);
1607:         rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1608:       }
1609:       ISRestoreIndices(label->points[s], &points);
1610:     }
1611:   }
1612:   /* Build SF that maps label points to remote processes */
1613:   PetscSectionCreate(comm, leafSection);
1614:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1615:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1616:   PetscFree(remoteOffsets);
1617:   /* Send the strata for each point over the derived SF */
1618:   PetscSectionGetStorageSize(*leafSection, &size);
1619:   PetscMalloc1(size, leafStrata);
1620:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE);
1621:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE);
1622:   /* Clean up */
1623:   PetscFree(rootStrata);
1624:   PetscFree(rootIdx);
1625:   PetscSectionDestroy(&rootSection);
1626:   PetscSFDestroy(&labelSF);
1627:   return 0;
1628: }

1630: /*@
1631:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1633:   Collective on sf

1635:   Input Parameters:
1636: + label - the DMLabel
1637: - sf    - the map from old to new distribution

1639:   Output Parameter:
1640: . labelnew - the new redistributed label

1642:   Level: intermediate

1644: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1645: @*/
1646: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1647: {
1648:   MPI_Comm     comm;
1649:   PetscSection leafSection;
1650:   PetscInt     p, pStart, pEnd, s, size, dof, offset, stratum;
1651:   PetscInt    *leafStrata, *strataIdx;
1652:   PetscInt   **points;
1653:   const char  *lname = NULL;
1654:   char        *name;
1655:   PetscInt     nameSize;
1656:   PetscHSetI   stratumHash;
1657:   size_t       len = 0;
1658:   PetscMPIInt  rank;

1661:   if (label) {
1663:     DMLabelMakeAllValid_Private(label);
1664:   }
1665:   PetscObjectGetComm((PetscObject)sf, &comm);
1666:   MPI_Comm_rank(comm, &rank);
1667:   /* Bcast name */
1668:   if (rank == 0) {
1669:     PetscObjectGetName((PetscObject)label, &lname);
1670:     PetscStrlen(lname, &len);
1671:   }
1672:   nameSize = len;
1673:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1674:   PetscMalloc1(nameSize + 1, &name);
1675:   if (rank == 0) PetscArraycpy(name, lname, nameSize + 1);
1676:   MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm);
1677:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1678:   PetscFree(name);
1679:   /* Bcast defaultValue */
1680:   if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1681:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1682:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1683:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1684:   /* Determine received stratum values and initialise new label*/
1685:   PetscHSetICreate(&stratumHash);
1686:   PetscSectionGetStorageSize(leafSection, &size);
1687:   for (p = 0; p < size; ++p) PetscHSetIAdd(stratumHash, leafStrata[p]);
1688:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1689:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1690:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1691:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1692:   /* Turn leafStrata into indices rather than stratum values */
1693:   offset = 0;
1694:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1695:   PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues);
1696:   for (s = 0; s < (*labelNew)->numStrata; ++s) PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1697:   for (p = 0; p < size; ++p) {
1698:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1699:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1700:         leafStrata[p] = s;
1701:         break;
1702:       }
1703:     }
1704:   }
1705:   /* Rebuild the point strata on the receiver */
1706:   PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes);
1707:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1708:   for (p = pStart; p < pEnd; p++) {
1709:     PetscSectionGetDof(leafSection, p, &dof);
1710:     PetscSectionGetOffset(leafSection, p, &offset);
1711:     for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1712:   }
1713:   PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht);
1714:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->points);
1715:   PetscMalloc1((*labelNew)->numStrata, &points);
1716:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1717:     PetscHSetICreate(&(*labelNew)->ht[s]);
1718:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1719:   }
1720:   /* Insert points into new strata */
1721:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1722:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1723:   for (p = pStart; p < pEnd; p++) {
1724:     PetscSectionGetDof(leafSection, p, &dof);
1725:     PetscSectionGetOffset(leafSection, p, &offset);
1726:     for (s = 0; s < dof; s++) {
1727:       stratum                               = leafStrata[offset + s];
1728:       points[stratum][strataIdx[stratum]++] = p;
1729:     }
1730:   }
1731:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1732:     ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &(points[s][0]), PETSC_OWN_POINTER, &((*labelNew)->points[s]));
1733:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices");
1734:   }
1735:   PetscFree(points);
1736:   PetscHSetIDestroy(&stratumHash);
1737:   PetscFree(leafStrata);
1738:   PetscFree(strataIdx);
1739:   PetscSectionDestroy(&leafSection);
1740:   return 0;
1741: }

1743: /*@
1744:   DMLabelGather - Gather all label values from leafs into roots

1746:   Collective on sf

1748:   Input Parameters:
1749: + label - the DMLabel
1750: - sf - the Star Forest point communication map

1752:   Output Parameters:
1753: . labelNew - the new DMLabel with localised leaf values

1755:   Level: developer

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

1759: .seealso: `DMLabelDistribute()`
1760: @*/
1761: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1762: {
1763:   MPI_Comm        comm;
1764:   PetscSection    rootSection;
1765:   PetscSF         sfLabel;
1766:   PetscSFNode    *rootPoints, *leafPoints;
1767:   PetscInt        p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1768:   const PetscInt *rootDegree, *ilocal;
1769:   PetscInt       *rootStrata;
1770:   const char     *lname;
1771:   char           *name;
1772:   PetscInt        nameSize;
1773:   size_t          len = 0;
1774:   PetscMPIInt     rank, size;

1778:   PetscObjectGetComm((PetscObject)sf, &comm);
1779:   MPI_Comm_rank(comm, &rank);
1780:   MPI_Comm_size(comm, &size);
1781:   /* Bcast name */
1782:   if (rank == 0) {
1783:     PetscObjectGetName((PetscObject)label, &lname);
1784:     PetscStrlen(lname, &len);
1785:   }
1786:   nameSize = len;
1787:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1788:   PetscMalloc1(nameSize + 1, &name);
1789:   if (rank == 0) PetscArraycpy(name, lname, nameSize + 1);
1790:   MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm);
1791:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1792:   PetscFree(name);
1793:   /* Gather rank/index pairs of leaves into local roots to build
1794:      an inverse, multi-rooted SF. Note that this ignores local leaf
1795:      indexing due to the use of the multiSF in PetscSFGather. */
1796:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1797:   PetscMalloc1(nroots, &leafPoints);
1798:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1799:   for (p = 0; p < nleaves; p++) {
1800:     PetscInt ilp = ilocal ? ilocal[p] : p;

1802:     leafPoints[ilp].index = ilp;
1803:     leafPoints[ilp].rank  = rank;
1804:   }
1805:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1806:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1807:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1808:   PetscMalloc1(nmultiroots, &rootPoints);
1809:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1810:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1811:   PetscSFCreate(comm, &sfLabel);
1812:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1813:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1814:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1815:   /* Rebuild the point strata on the receiver */
1816:   for (p = 0, idx = 0; p < nroots; p++) {
1817:     for (d = 0; d < rootDegree[p]; d++) {
1818:       PetscSectionGetDof(rootSection, idx + d, &dof);
1819:       PetscSectionGetOffset(rootSection, idx + d, &offset);
1820:       for (s = 0; s < dof; s++) DMLabelSetValue(*labelNew, p, rootStrata[offset + s]);
1821:     }
1822:     idx += rootDegree[p];
1823:   }
1824:   PetscFree(leafPoints);
1825:   PetscFree(rootStrata);
1826:   PetscSectionDestroy(&rootSection);
1827:   PetscSFDestroy(&sfLabel);
1828:   return 0;
1829: }

1831: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
1832: {
1833:   const PetscInt *degree;
1834:   const PetscInt *points;
1835:   PetscInt        Nr, r, Nl, l, val, defVal;

1837:   DMLabelGetDefaultValue(label, &defVal);
1838:   /* Add in leaves */
1839:   PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL);
1840:   for (l = 0; l < Nl; ++l) {
1841:     DMLabelGetValue(label, points[l], &val);
1842:     if (val != defVal) valArray[points[l]] = val;
1843:   }
1844:   /* Add in shared roots */
1845:   PetscSFComputeDegreeBegin(pointSF, &degree);
1846:   PetscSFComputeDegreeEnd(pointSF, &degree);
1847:   for (r = 0; r < Nr; ++r) {
1848:     if (degree[r]) {
1849:       DMLabelGetValue(label, r, &val);
1850:       if (val != defVal) valArray[r] = val;
1851:     }
1852:   }
1853:   return 0;
1854: }

1856: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1857: {
1858:   const PetscInt *degree;
1859:   const PetscInt *points;
1860:   PetscInt        Nr, r, Nl, l, val, defVal;

1862:   DMLabelGetDefaultValue(label, &defVal);
1863:   /* Read out leaves */
1864:   PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL);
1865:   for (l = 0; l < Nl; ++l) {
1866:     const PetscInt p    = points[l];
1867:     const PetscInt cval = valArray[p];

1869:     if (cval != defVal) {
1870:       DMLabelGetValue(label, p, &val);
1871:       if (val == defVal) {
1872:         DMLabelSetValue(label, p, cval);
1873:         if (markPoint) (*markPoint)(label, p, cval, ctx);
1874:       }
1875:     }
1876:   }
1877:   /* Read out shared roots */
1878:   PetscSFComputeDegreeBegin(pointSF, &degree);
1879:   PetscSFComputeDegreeEnd(pointSF, &degree);
1880:   for (r = 0; r < Nr; ++r) {
1881:     if (degree[r]) {
1882:       const PetscInt cval = valArray[r];

1884:       if (cval != defVal) {
1885:         DMLabelGetValue(label, r, &val);
1886:         if (val == defVal) {
1887:           DMLabelSetValue(label, r, cval);
1888:           if (markPoint) (*markPoint)(label, r, cval, ctx);
1889:         }
1890:       }
1891:     }
1892:   }
1893:   return 0;
1894: }

1896: /*@
1897:   DMLabelPropagateBegin - Setup a cycle of label propagation

1899:   Collective on sf

1901:   Input Parameters:
1902: + label - The DMLabel to propagate across processes
1903: - sf    - The SF describing parallel layout of the label points

1905:   Level: intermediate

1907: .seealso: `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
1908: @*/
1909: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
1910: {
1911:   PetscInt    Nr, r, defVal;
1912:   PetscMPIInt size;

1914:   MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size);
1915:   if (size > 1) {
1916:     DMLabelGetDefaultValue(label, &defVal);
1917:     PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL);
1918:     if (Nr >= 0) PetscMalloc1(Nr, &label->propArray);
1919:     for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
1920:   }
1921:   return 0;
1922: }

1924: /*@
1925:   DMLabelPropagateEnd - Tear down a cycle of label propagation

1927:   Collective on sf

1929:   Input Parameters:
1930: + label - The DMLabel to propagate across processes
1931: - sf    - The SF describing parallel layout of the label points

1933:   Level: intermediate

1935: .seealso: `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
1936: @*/
1937: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
1938: {
1939:   PetscFree(label->propArray);
1940:   label->propArray = NULL;
1941:   return 0;
1942: }

1944: /*@C
1945:   DMLabelPropagatePush - Tear down a cycle of label propagation

1947:   Collective on sf

1949:   Input Parameters:
1950: + label     - The DMLabel to propagate across processes
1951: . sf        - The SF describing parallel layout of the label points
1952: . markPoint - An optional callback that is called when a point is marked, or NULL
1953: - ctx       - An optional user context for the callback, or NULL

1955:   Calling sequence of markPoint:
1956: $ markPoint(DMLabel label, PetscInt p, PetscInt val, void *ctx);

1958: + label - The DMLabel
1959: . p     - The point being marked
1960: . val   - The label value for p
1961: - ctx   - An optional user context

1963:   Level: intermediate

1965: .seealso: `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
1966: @*/
1967: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1968: {
1969:   PetscInt   *valArray = label->propArray, Nr;
1970:   PetscMPIInt size;

1972:   MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size);
1973:   PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL);
1974:   if (size > 1 && Nr >= 0) {
1975:     /* Communicate marked edges
1976:        The current implementation allocates an array the size of the number of root. We put the label values into the
1977:        array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.

1979:        TODO: We could use in-place communication with a different SF
1980:        We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
1981:        already been marked. If not, it might have been handled on the process in this round, but we add it anyway.

1983:        In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
1984:        values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new
1985:        edge to the queue.
1986:     */
1987:     DMLabelPropagateInit_Internal(label, pointSF, valArray);
1988:     PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX);
1989:     PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX);
1990:     PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE);
1991:     PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE);
1992:     DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx);
1993:   }
1994:   return 0;
1995: }

1997: /*@
1998:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

2000:   Not collective

2002:   Input Parameter:
2003: . label - the DMLabel

2005:   Output Parameters:
2006: + section - the section giving offsets for each stratum
2007: - is - An IS containing all the label points

2009:   Level: developer

2011: .seealso: `DMLabelDistribute()`
2012: @*/
2013: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2014: {
2015:   IS              vIS;
2016:   const PetscInt *values;
2017:   PetscInt       *points;
2018:   PetscInt        nV, vS = 0, vE = 0, v, N;

2021:   DMLabelGetNumValues(label, &nV);
2022:   DMLabelGetValueIS(label, &vIS);
2023:   ISGetIndices(vIS, &values);
2024:   if (nV) {
2025:     vS = values[0];
2026:     vE = values[0] + 1;
2027:   }
2028:   for (v = 1; v < nV; ++v) {
2029:     vS = PetscMin(vS, values[v]);
2030:     vE = PetscMax(vE, values[v] + 1);
2031:   }
2032:   PetscSectionCreate(PETSC_COMM_SELF, section);
2033:   PetscSectionSetChart(*section, vS, vE);
2034:   for (v = 0; v < nV; ++v) {
2035:     PetscInt n;

2037:     DMLabelGetStratumSize(label, values[v], &n);
2038:     PetscSectionSetDof(*section, values[v], n);
2039:   }
2040:   PetscSectionSetUp(*section);
2041:   PetscSectionGetStorageSize(*section, &N);
2042:   PetscMalloc1(N, &points);
2043:   for (v = 0; v < nV; ++v) {
2044:     IS              is;
2045:     const PetscInt *spoints;
2046:     PetscInt        dof, off, p;

2048:     PetscSectionGetDof(*section, values[v], &dof);
2049:     PetscSectionGetOffset(*section, values[v], &off);
2050:     DMLabelGetStratumIS(label, values[v], &is);
2051:     ISGetIndices(is, &spoints);
2052:     for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2053:     ISRestoreIndices(is, &spoints);
2054:     ISDestroy(&is);
2055:   }
2056:   ISRestoreIndices(vIS, &values);
2057:   ISDestroy(&vIS);
2058:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
2059:   return 0;
2060: }

2062: /*@
2063:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2064:   the local section and an SF describing the section point overlap.

2066:   Collective on sf

2068:   Input Parameters:
2069: + s - The PetscSection for the local field layout
2070: . sf - The SF describing parallel layout of the section points
2071: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
2072: . label - The label specifying the points
2073: - labelValue - The label stratum specifying the points

2075:   Output Parameter:
2076: . gsection - The PetscSection for the global field layout

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

2080:   Level: developer

2082: .seealso: `PetscSectionCreate()`
2083: @*/
2084: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2085: {
2086:   PetscInt *neg = NULL, *tmpOff = NULL;
2087:   PetscInt  pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

2092:   PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection);
2093:   PetscSectionGetChart(s, &pStart, &pEnd);
2094:   PetscSectionSetChart(*gsection, pStart, pEnd);
2095:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
2096:   if (nroots >= 0) {
2098:     PetscCalloc1(nroots, &neg);
2099:     if (nroots > pEnd - pStart) {
2100:       PetscCalloc1(nroots, &tmpOff);
2101:     } else {
2102:       tmpOff = &(*gsection)->atlasDof[-pStart];
2103:     }
2104:   }
2105:   /* Mark ghost points with negative dof */
2106:   for (p = pStart; p < pEnd; ++p) {
2107:     PetscInt value;

2109:     DMLabelGetValue(label, p, &value);
2110:     if (value != labelValue) continue;
2111:     PetscSectionGetDof(s, p, &dof);
2112:     PetscSectionSetDof(*gsection, p, dof);
2113:     PetscSectionGetConstraintDof(s, p, &cdof);
2114:     if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
2115:     if (neg) neg[p] = -(dof + 1);
2116:   }
2117:   PetscSectionSetUpBC(*gsection);
2118:   if (nroots >= 0) {
2119:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2120:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2121:     if (nroots > pEnd - pStart) {
2122:       for (p = pStart; p < pEnd; ++p) {
2123:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2124:       }
2125:     }
2126:   }
2127:   /* Calculate new sizes, get process offset, and calculate point offsets */
2128:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2129:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2130:     (*gsection)->atlasOff[p] = off;
2131:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2132:   }
2133:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
2134:   globalOff -= off;
2135:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2136:     (*gsection)->atlasOff[p] += globalOff;
2137:     if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2138:   }
2139:   /* Put in negative offsets for ghost points */
2140:   if (nroots >= 0) {
2141:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2142:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2143:     if (nroots > pEnd - pStart) {
2144:       for (p = pStart; p < pEnd; ++p) {
2145:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2146:       }
2147:     }
2148:   }
2149:   if (nroots >= 0 && nroots > pEnd - pStart) PetscFree(tmpOff);
2150:   PetscFree(neg);
2151:   return 0;
2152: }

2154: typedef struct _n_PetscSectionSym_Label {
2155:   DMLabel              label;
2156:   PetscCopyMode       *modes;
2157:   PetscInt            *sizes;
2158:   const PetscInt    ***perms;
2159:   const PetscScalar ***rots;
2160:   PetscInt (*minMaxOrients)[2];
2161:   PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2162: } PetscSectionSym_Label;

2164: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2165: {
2166:   PetscInt               i, j;
2167:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;

2169:   for (i = 0; i <= sl->numStrata; i++) {
2170:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2171:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2172:         if (sl->perms[i]) PetscFree(sl->perms[i][j]);
2173:         if (sl->rots[i]) PetscFree(sl->rots[i][j]);
2174:       }
2175:       if (sl->perms[i]) {
2176:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

2178:         PetscFree(perms);
2179:       }
2180:       if (sl->rots[i]) {
2181:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

2183:         PetscFree(rots);
2184:       }
2185:     }
2186:   }
2187:   PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients);
2188:   DMLabelDestroy(&sl->label);
2189:   sl->numStrata = 0;
2190:   return 0;
2191: }

2193: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2194: {
2195:   PetscSectionSymLabelReset(sym);
2196:   PetscFree(sym->data);
2197:   return 0;
2198: }

2200: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2201: {
2202:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2203:   PetscBool              isAscii;
2204:   DMLabel                label = sl->label;
2205:   const char            *name;

2207:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii);
2208:   if (isAscii) {
2209:     PetscInt          i, j, k;
2210:     PetscViewerFormat format;

2212:     PetscViewerGetFormat(viewer, &format);
2213:     if (label) {
2214:       PetscViewerGetFormat(viewer, &format);
2215:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2216:         PetscViewerASCIIPushTab(viewer);
2217:         DMLabelView(label, viewer);
2218:         PetscViewerASCIIPopTab(viewer);
2219:       } else {
2220:         PetscObjectGetName((PetscObject)sl->label, &name);
2221:         PetscViewerASCIIPrintf(viewer, "  Label '%s'\n", name);
2222:       }
2223:     } else {
2224:       PetscViewerASCIIPrintf(viewer, "No label given\n");
2225:     }
2226:     PetscViewerASCIIPushTab(viewer);
2227:     for (i = 0; i <= sl->numStrata; i++) {
2228:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

2230:       if (!(sl->perms[i] || sl->rots[i])) {
2231:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]);
2232:       } else {
2233:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]);
2234:         PetscViewerASCIIPushTab(viewer);
2235:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
2236:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2237:           PetscViewerASCIIPushTab(viewer);
2238:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2239:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2240:               PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j);
2241:             } else {
2242:               PetscInt tab;

2244:               PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j);
2245:               PetscViewerASCIIPushTab(viewer);
2246:               PetscViewerASCIIGetTab(viewer, &tab);
2247:               if (sl->perms[i] && sl->perms[i][j]) {
2248:                 PetscViewerASCIIPrintf(viewer, "Permutation:");
2249:                 PetscViewerASCIISetTab(viewer, 0);
2250:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]);
2251:                 PetscViewerASCIIPrintf(viewer, "\n");
2252:                 PetscViewerASCIISetTab(viewer, tab);
2253:               }
2254:               if (sl->rots[i] && sl->rots[i][j]) {
2255:                 PetscViewerASCIIPrintf(viewer, "Rotations:  ");
2256:                 PetscViewerASCIISetTab(viewer, 0);
2257: #if defined(PETSC_USE_COMPLEX)
2258:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %+g+i*%+g", (double)PetscRealPart(sl->rots[i][j][k]), (double)PetscImaginaryPart(sl->rots[i][j][k]));
2259: #else
2260:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]);
2261: #endif
2262:                 PetscViewerASCIIPrintf(viewer, "\n");
2263:                 PetscViewerASCIISetTab(viewer, tab);
2264:               }
2265:               PetscViewerASCIIPopTab(viewer);
2266:             }
2267:           }
2268:           PetscViewerASCIIPopTab(viewer);
2269:         }
2270:         PetscViewerASCIIPopTab(viewer);
2271:       }
2272:     }
2273:     PetscViewerASCIIPopTab(viewer);
2274:   }
2275:   return 0;
2276: }

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

2281:   Logically collective on sym

2283:   Input parameters:
2284: + sym - the section symmetries
2285: - label - the DMLabel describing the types of points

2287:   Level: developer:

2289: .seealso: `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2290: @*/
2291: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2292: {
2293:   PetscSectionSym_Label *sl;

2296:   sl = (PetscSectionSym_Label *)sym->data;
2297:   if (sl->label && sl->label != label) PetscSectionSymLabelReset(sym);
2298:   if (label) {
2299:     sl->label = label;
2300:     PetscObjectReference((PetscObject)label);
2301:     DMLabelGetNumValues(label, &sl->numStrata);
2302:     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);
2303:     PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode));
2304:     PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt));
2305:     PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **));
2306:     PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **));
2307:     PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2]));
2308:   }
2309:   return 0;
2310: }

2312: /*@C
2313:   PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum

2315:   Logically collective on sym

2317:   Input Parameters:
2318: + sym       - the section symmetries
2319: - stratum   - the stratum value in the label that we are assigning symmetries for

2321:   Output Parameters:
2322: + size      - the number of dofs for points in the stratum of the label
2323: . minOrient - the smallest orientation for a point in this stratum
2324: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2325: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2326: - 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

2328:   Level: developer

2330: .seealso: `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2331: @*/
2332: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2333: {
2334:   PetscSectionSym_Label *sl;
2335:   const char            *name;
2336:   PetscInt               i;

2339:   sl = (PetscSectionSym_Label *)sym->data;
2341:   for (i = 0; i <= sl->numStrata; i++) {
2342:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2344:     if (stratum == value) break;
2345:   }
2346:   PetscObjectGetName((PetscObject)sl->label, &name);
2348:   if (size) {
2350:     *size = sl->sizes[i];
2351:   }
2352:   if (minOrient) {
2354:     *minOrient = sl->minMaxOrients[i][0];
2355:   }
2356:   if (maxOrient) {
2358:     *maxOrient = sl->minMaxOrients[i][1];
2359:   }
2360:   if (perms) {
2362:     *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;
2363:   }
2364:   if (rots) {
2366:     *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;
2367:   }
2368:   return 0;
2369: }

2371: /*@C
2372:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

2374:   Logically collective on sym

2376:   InputParameters:
2377: + sym       - the section symmetries
2378: . stratum   - the stratum value in the label that we are assigning symmetries for
2379: . size      - the number of dofs for points in the stratum of the label
2380: . minOrient - the smallest orientation for a point in this stratum
2381: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2382: . mode      - how sym should copy the perms and rots arrays
2383: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2384: - 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

2386:   Level: developer

2388: .seealso: `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2389: @*/
2390: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2391: {
2392:   PetscSectionSym_Label *sl;
2393:   const char            *name;
2394:   PetscInt               i, j, k;

2397:   sl = (PetscSectionSym_Label *)sym->data;
2399:   for (i = 0; i <= sl->numStrata; i++) {
2400:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2402:     if (stratum == value) break;
2403:   }
2404:   PetscObjectGetName((PetscObject)sl->label, &name);
2406:   sl->sizes[i]            = size;
2407:   sl->modes[i]            = mode;
2408:   sl->minMaxOrients[i][0] = minOrient;
2409:   sl->minMaxOrients[i][1] = maxOrient;
2410:   if (mode == PETSC_COPY_VALUES) {
2411:     if (perms) {
2412:       PetscInt **ownPerms;

2414:       PetscCalloc1(maxOrient - minOrient, &ownPerms);
2415:       for (j = 0; j < maxOrient - minOrient; j++) {
2416:         if (perms[j]) {
2417:           PetscMalloc1(size, &ownPerms[j]);
2418:           for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2419:         }
2420:       }
2421:       sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2422:     }
2423:     if (rots) {
2424:       PetscScalar **ownRots;

2426:       PetscCalloc1(maxOrient - minOrient, &ownRots);
2427:       for (j = 0; j < maxOrient - minOrient; j++) {
2428:         if (rots[j]) {
2429:           PetscMalloc1(size, &ownRots[j]);
2430:           for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2431:         }
2432:       }
2433:       sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2434:     }
2435:   } else {
2436:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2437:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2438:   }
2439:   return 0;
2440: }

2442: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2443: {
2444:   PetscInt               i, j, numStrata;
2445:   PetscSectionSym_Label *sl;
2446:   DMLabel                label;

2448:   sl        = (PetscSectionSym_Label *)sym->data;
2449:   numStrata = sl->numStrata;
2450:   label     = sl->label;
2451:   for (i = 0; i < numPoints; i++) {
2452:     PetscInt point = points[2 * i];
2453:     PetscInt ornt  = points[2 * i + 1];

2455:     for (j = 0; j < numStrata; j++) {
2456:       if (label->validIS[j]) {
2457:         PetscInt k;

2459:         ISLocate(label->points[j], point, &k);
2460:         if (k >= 0) break;
2461:       } else {
2462:         PetscBool has;

2464:         PetscHSetIHas(label->ht[j], point, &has);
2465:         if (has) break;
2466:       }
2467:     }
2469:                j < numStrata ? label->stratumValues[j] : label->defaultValue);
2470:     if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2471:     if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2472:   }
2473:   return 0;
2474: }

2476: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2477: {
2478:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2479:   IS                     valIS;
2480:   const PetscInt        *values;
2481:   PetscInt               Nv, v;

2483:   DMLabelGetNumValues(sl->label, &Nv);
2484:   DMLabelGetValueIS(sl->label, &valIS);
2485:   ISGetIndices(valIS, &values);
2486:   for (v = 0; v < Nv; ++v) {
2487:     const PetscInt      val = values[v];
2488:     PetscInt            size, minOrient, maxOrient;
2489:     const PetscInt    **perms;
2490:     const PetscScalar **rots;

2492:     PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots);
2493:     PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots);
2494:   }
2495:   ISDestroy(&valIS);
2496:   return 0;
2497: }

2499: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2500: {
2501:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2502:   DMLabel                dlabel;

2504:   DMLabelDistribute(sl->label, migrationSF, &dlabel);
2505:   PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym);
2506:   DMLabelDestroy(&dlabel);
2507:   PetscSectionSymCopy(sym, *dsym);
2508:   return 0;
2509: }

2511: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2512: {
2513:   PetscSectionSym_Label *sl;

2515:   PetscNew(&sl);
2516:   sym->ops->getpoints  = PetscSectionSymGetPoints_Label;
2517:   sym->ops->distribute = PetscSectionSymDistribute_Label;
2518:   sym->ops->copy       = PetscSectionSymCopy_Label;
2519:   sym->ops->view       = PetscSectionSymView_Label;
2520:   sym->ops->destroy    = PetscSectionSymDestroy_Label;
2521:   sym->data            = (void *)sl;
2522:   return 0;
2523: }

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

2528:   Collective

2530:   Input Parameters:
2531: + comm - the MPI communicator for the new symmetry
2532: - label - the label defining the strata

2534:   Output Parameters:
2535: . sym - the section symmetries

2537:   Level: developer

2539: .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2540: @*/
2541: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2542: {
2543:   DMInitializePackage();
2544:   PetscSectionSymCreate(comm, sym);
2545:   PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL);
2546:   PetscSectionSymLabelSetLabel(*sym, label);
2547:   return 0;
2548: }