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: PetscFunctionList DMLabelList              = NULL;
      8: PetscBool         DMLabelRegisterAllCalled = PETSC_FALSE;

     10: /*@
     11:   DMLabelCreate - Create a `DMLabel` object, which is a multimap

     13:   Collective

     15:   Input Parameters:
     16: + comm - The communicator, usually `PETSC_COMM_SELF`
     17: - name - The label name

     19:   Output Parameter:
     20: . label - The `DMLabel`

     22:   Level: beginner

     24:   Notes:
     25:   The label name is actually usually the `PetscObject` name.
     26:   One can get/set it with `PetscObjectGetName()`/`PetscObjectSetName()`.

     28: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`
     29: @*/
     30: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
     31: {
     32:   PetscFunctionBegin;
     33:   PetscAssertPointer(label, 3);
     34:   PetscCall(DMInitializePackage());

     36:   PetscCall(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:   PetscCall(PetscHMapICreate(&(*label)->hmap));
     48:   PetscCall(PetscObjectSetName((PetscObject)*label, name));
     49:   PetscCall(DMLabelSetType(*label, DMLABELCONCRETE));
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: /*@
     54:   DMLabelSetUp - SetUp a `DMLabel` object

     56:   Collective

     58:   Input Parameters:
     59: . label - The `DMLabel`

     61:   Level: intermediate

     63: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
     64: @*/
     65: PetscErrorCode DMLabelSetUp(DMLabel label)
     66: {
     67:   PetscFunctionBegin;
     69:   PetscTryTypeMethod(label, setup);
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: /*
     74:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

     76:   Not collective

     78:   Input parameter:
     79: + label - The `DMLabel`
     80: - v - The stratum value

     82:   Output parameter:
     83: . label - The `DMLabel` with stratum in sorted list format

     85:   Level: developer

     87: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
     88: */
     89: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
     90: {
     91:   IS       is;
     92:   PetscInt off = 0, *pointArray, p;

     94:   PetscFunctionBegin;
     95:   if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
     96:   PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
     97:   PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
     98:   PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
     99:   PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
    100:   PetscCall(PetscHSetIClear(label->ht[v]));
    101:   PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
    102:   if (label->bt) {
    103:     for (p = 0; p < label->stratumSizes[v]; ++p) {
    104:       const PetscInt point = pointArray[p];
    105:       PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
    106:       PetscCall(PetscBTSet(label->bt, point - label->pStart));
    107:     }
    108:   }
    109:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
    110:     PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
    111:     PetscCall(PetscFree(pointArray));
    112:   } else {
    113:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
    114:   }
    115:   PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, PETSC_TRUE));
    116:   PetscCall(PetscObjectSetName((PetscObject)is, "indices"));
    117:   label->points[v]  = is;
    118:   label->validIS[v] = PETSC_TRUE;
    119:   PetscCall(PetscObjectStateIncrease((PetscObject)label));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /*
    124:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

    126:   Not Collective

    128:   Input parameter:
    129: . label - The `DMLabel`

    131:   Output parameter:
    132: . label - The `DMLabel` with all strata in sorted list format

    134:   Level: developer

    136: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
    137: */
    138: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
    139: {
    140:   PetscInt v;

    142:   PetscFunctionBegin;
    143:   for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeValid_Private(label, v));
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }

    147: /*
    148:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

    150:   Not Collective

    152:   Input parameter:
    153: + label - The `DMLabel`
    154: - v - The stratum value

    156:   Output parameter:
    157: . label - The `DMLabel` with stratum in hash format

    159:   Level: developer

    161: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
    162: */
    163: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
    164: {
    165:   PetscInt        p;
    166:   const PetscInt *points;

    168:   PetscFunctionBegin;
    169:   if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
    170:   PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeInvalid_Private", v);
    171:   if (label->points[v]) {
    172:     PetscCall(ISGetIndices(label->points[v], &points));
    173:     for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
    174:     PetscCall(ISRestoreIndices(label->points[v], &points));
    175:     PetscCall(ISDestroy(&label->points[v]));
    176:   }
    177:   label->validIS[v] = PETSC_FALSE;
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

    181: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
    182: {
    183:   PetscInt v;

    185:   PetscFunctionBegin;
    186:   for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeInvalid_Private(label, v));
    187:   PetscFunctionReturn(PETSC_SUCCESS);
    188: }

    190: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
    191:   #define DMLABEL_LOOKUP_THRESHOLD 16
    192: #endif

    194: PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
    195: {
    196:   PetscInt v;

    198:   PetscFunctionBegin;
    199:   *index = -1;
    200:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD || label->readonly) {
    201:     for (v = 0; v < label->numStrata; ++v)
    202:       if (label->stratumValues[v] == value) {
    203:         *index = v;
    204:         break;
    205:       }
    206:   } else {
    207:     PetscCall(PetscHMapIGet(label->hmap, value, index));
    208:   }
    209:   if (PetscDefined(USE_DEBUG) && !label->readonly) { /* Check strata hash map consistency */
    210:     PetscInt len, loc = -1;
    211:     PetscCall(PetscHMapIGetSize(label->hmap, &len));
    212:     PetscCheck(len == label->numStrata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
    213:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
    214:       PetscCall(PetscHMapIGet(label->hmap, value, &loc));
    215:     } else {
    216:       for (v = 0; v < label->numStrata; ++v)
    217:         if (label->stratumValues[v] == value) {
    218:           loc = v;
    219:           break;
    220:         }
    221:     }
    222:     PetscCheck(loc == *index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
    223:   }
    224:   PetscFunctionReturn(PETSC_SUCCESS);
    225: }

    227: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
    228: {
    229:   PetscInt    v;
    230:   PetscInt   *tmpV;
    231:   PetscInt   *tmpS;
    232:   PetscHSetI *tmpH, ht;
    233:   IS         *tmpP, is;
    234:   PetscBool  *tmpB;
    235:   PetscHMapI  hmap = label->hmap;

    237:   PetscFunctionBegin;
    238:   v    = label->numStrata;
    239:   tmpV = label->stratumValues;
    240:   tmpS = label->stratumSizes;
    241:   tmpH = label->ht;
    242:   tmpP = label->points;
    243:   tmpB = label->validIS;
    244:   { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free  */
    245:     PetscInt   *oldV = tmpV;
    246:     PetscInt   *oldS = tmpS;
    247:     PetscHSetI *oldH = tmpH;
    248:     IS         *oldP = tmpP;
    249:     PetscBool  *oldB = tmpB;
    250:     PetscCall(PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV));
    251:     PetscCall(PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS));
    252:     PetscCall(PetscCalloc((v + 1) * sizeof(*tmpH), &tmpH));
    253:     PetscCall(PetscCalloc((v + 1) * sizeof(*tmpP), &tmpP));
    254:     PetscCall(PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB));
    255:     PetscCall(PetscArraycpy(tmpV, oldV, v));
    256:     PetscCall(PetscArraycpy(tmpS, oldS, v));
    257:     PetscCall(PetscArraycpy(tmpH, oldH, v));
    258:     PetscCall(PetscArraycpy(tmpP, oldP, v));
    259:     PetscCall(PetscArraycpy(tmpB, oldB, v));
    260:     PetscCall(PetscFree(oldV));
    261:     PetscCall(PetscFree(oldS));
    262:     PetscCall(PetscFree(oldH));
    263:     PetscCall(PetscFree(oldP));
    264:     PetscCall(PetscFree(oldB));
    265:   }
    266:   label->numStrata     = v + 1;
    267:   label->stratumValues = tmpV;
    268:   label->stratumSizes  = tmpS;
    269:   label->ht            = tmpH;
    270:   label->points        = tmpP;
    271:   label->validIS       = tmpB;
    272:   PetscCall(PetscHSetICreate(&ht));
    273:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
    274:   PetscCall(PetscHMapISet(hmap, value, v));
    275:   tmpV[v] = value;
    276:   tmpS[v] = 0;
    277:   tmpH[v] = ht;
    278:   tmpP[v] = is;
    279:   tmpB[v] = PETSC_TRUE;
    280:   PetscCall(PetscObjectStateIncrease((PetscObject)label));
    281:   *index = v;
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
    286: {
    287:   PetscFunctionBegin;
    288:   PetscCall(DMLabelLookupStratum(label, value, index));
    289:   if (*index < 0) PetscCall(DMLabelNewStratum(label, value, index));
    290:   PetscFunctionReturn(PETSC_SUCCESS);
    291: }

    293: PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
    294: {
    295:   PetscFunctionBegin;
    296:   *size = 0;
    297:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
    298:   if (label->readonly || label->validIS[v]) {
    299:     *size = label->stratumSizes[v];
    300:   } else {
    301:     PetscCall(PetscHSetIGetSize(label->ht[v], size));
    302:   }
    303:   PetscFunctionReturn(PETSC_SUCCESS);
    304: }

    306: /*@
    307:   DMLabelAddStratum - Adds a new stratum value in a `DMLabel`

    309:   Input Parameters:
    310: + label - The `DMLabel`
    311: - value - The stratum value

    313:   Level: beginner

    315: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
    316: @*/
    317: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
    318: {
    319:   PetscInt v;

    321:   PetscFunctionBegin;
    323:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    324:   PetscCall(DMLabelLookupAddStratum(label, value, &v));
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*@
    329:   DMLabelAddStrata - Adds new stratum values in a `DMLabel`

    331:   Not Collective

    333:   Input Parameters:
    334: + label         - The `DMLabel`
    335: . numStrata     - The number of stratum values
    336: - stratumValues - The stratum values

    338:   Level: beginner

    340: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
    341: @*/
    342: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
    343: {
    344:   PetscInt *values, v;

    346:   PetscFunctionBegin;
    348:   if (numStrata) PetscAssertPointer(stratumValues, 3);
    349:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    350:   PetscCall(PetscMalloc1(numStrata, &values));
    351:   PetscCall(PetscArraycpy(values, stratumValues, numStrata));
    352:   PetscCall(PetscSortRemoveDupsInt(&numStrata, values));
    353:   if (!label->numStrata) { /* Fast preallocation */
    354:     PetscInt   *tmpV;
    355:     PetscInt   *tmpS;
    356:     PetscHSetI *tmpH, ht;
    357:     IS         *tmpP, is;
    358:     PetscBool  *tmpB;
    359:     PetscHMapI  hmap = label->hmap;

    361:     PetscCall(PetscMalloc1(numStrata, &tmpV));
    362:     PetscCall(PetscMalloc1(numStrata, &tmpS));
    363:     PetscCall(PetscCalloc1(numStrata, &tmpH));
    364:     PetscCall(PetscCalloc1(numStrata, &tmpP));
    365:     PetscCall(PetscMalloc1(numStrata, &tmpB));
    366:     label->numStrata     = numStrata;
    367:     label->stratumValues = tmpV;
    368:     label->stratumSizes  = tmpS;
    369:     label->ht            = tmpH;
    370:     label->points        = tmpP;
    371:     label->validIS       = tmpB;
    372:     for (v = 0; v < numStrata; ++v) {
    373:       PetscCall(PetscHSetICreate(&ht));
    374:       PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
    375:       PetscCall(PetscHMapISet(hmap, values[v], v));
    376:       tmpV[v] = values[v];
    377:       tmpS[v] = 0;
    378:       tmpH[v] = ht;
    379:       tmpP[v] = is;
    380:       tmpB[v] = PETSC_TRUE;
    381:     }
    382:     PetscCall(PetscObjectStateIncrease((PetscObject)label));
    383:   } else {
    384:     for (v = 0; v < numStrata; ++v) PetscCall(DMLabelAddStratum(label, values[v]));
    385:   }
    386:   PetscCall(PetscFree(values));
    387:   PetscFunctionReturn(PETSC_SUCCESS);
    388: }

    390: /*@
    391:   DMLabelAddStrataIS - Adds new stratum values in a `DMLabel`

    393:   Not Collective

    395:   Input Parameters:
    396: + label   - The `DMLabel`
    397: - valueIS - Index set with stratum values

    399:   Level: beginner

    401: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
    402: @*/
    403: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
    404: {
    405:   PetscInt        numStrata;
    406:   const PetscInt *stratumValues;

    408:   PetscFunctionBegin;
    411:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    412:   PetscCall(ISGetLocalSize(valueIS, &numStrata));
    413:   PetscCall(ISGetIndices(valueIS, &stratumValues));
    414:   PetscCall(DMLabelAddStrata(label, numStrata, stratumValues));
    415:   PetscFunctionReturn(PETSC_SUCCESS);
    416: }

    418: static PetscErrorCode DMLabelView_Concrete_Ascii(DMLabel label, PetscViewer viewer)
    419: {
    420:   PetscInt    v;
    421:   PetscMPIInt rank;

    423:   PetscFunctionBegin;
    424:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
    425:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
    426:   if (label) {
    427:     const char *name;

    429:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
    430:     PetscCall(PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name));
    431:     if (label->bt) PetscCall(PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd));
    432:     for (v = 0; v < label->numStrata; ++v) {
    433:       const PetscInt  value = label->stratumValues[v];
    434:       const PetscInt *points;
    435:       PetscInt        p;

    437:       PetscCall(ISGetIndices(label->points[v], &points));
    438:       for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
    439:       PetscCall(ISRestoreIndices(label->points[v], &points));
    440:     }
    441:   }
    442:   PetscCall(PetscViewerFlush(viewer));
    443:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    444:   PetscFunctionReturn(PETSC_SUCCESS);
    445: }

    447: static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
    448: {
    449:   PetscBool isascii;

    451:   PetscFunctionBegin;
    452:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    453:   if (isascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
    454:   PetscFunctionReturn(PETSC_SUCCESS);
    455: }

    457: /*@
    458:   DMLabelView - View the label

    460:   Collective

    462:   Input Parameters:
    463: + label  - The `DMLabel`
    464: - viewer - The `PetscViewer`

    466:   Level: intermediate

    468: .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
    469: @*/
    470: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
    471: {
    472:   PetscFunctionBegin;
    474:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer));
    476:   PetscCall(DMLabelMakeAllValid_Private(label));
    477:   PetscUseTypeMethod(label, view, viewer);
    478:   PetscFunctionReturn(PETSC_SUCCESS);
    479: }

    481: /*@
    482:   DMLabelReset - Destroys internal data structures in a `DMLabel`

    484:   Not Collective

    486:   Input Parameter:
    487: . label - The `DMLabel`

    489:   Level: beginner

    491: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
    492: @*/
    493: PetscErrorCode DMLabelReset(DMLabel label)
    494: {
    495:   PetscInt v;

    497:   PetscFunctionBegin;
    499:   for (v = 0; v < label->numStrata; ++v) {
    500:     if (label->ht[v]) PetscCall(PetscHSetIDestroy(&label->ht[v]));
    501:     PetscCall(ISDestroy(&label->points[v]));
    502:   }
    503:   label->numStrata = 0;
    504:   PetscCall(PetscFree(label->stratumValues));
    505:   PetscCall(PetscFree(label->stratumSizes));
    506:   PetscCall(PetscFree(label->ht));
    507:   PetscCall(PetscFree(label->points));
    508:   PetscCall(PetscFree(label->validIS));
    509:   PetscCall(PetscHMapIReset(label->hmap));
    510:   label->pStart = -1;
    511:   label->pEnd   = -1;
    512:   PetscCall(PetscBTDestroy(&label->bt));
    513:   PetscFunctionReturn(PETSC_SUCCESS);
    514: }

    516: /*@
    517:   DMLabelDestroy - Destroys a `DMLabel`

    519:   Collective

    521:   Input Parameter:
    522: . label - The `DMLabel`

    524:   Level: beginner

    526: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
    527: @*/
    528: PetscErrorCode DMLabelDestroy(DMLabel *label)
    529: {
    530:   PetscFunctionBegin;
    531:   if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
    533:   if (--((PetscObject)*label)->refct > 0) {
    534:     *label = NULL;
    535:     PetscFunctionReturn(PETSC_SUCCESS);
    536:   }
    537:   PetscCall(DMLabelReset(*label));
    538:   PetscCall(PetscHMapIDestroy(&(*label)->hmap));
    539:   PetscCall(PetscHeaderDestroy(label));
    540:   PetscFunctionReturn(PETSC_SUCCESS);
    541: }

    543: static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
    544: {
    545:   PetscFunctionBegin;
    546:   for (PetscInt v = 0; v < label->numStrata; ++v) {
    547:     PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
    548:     PetscCall(PetscObjectReference((PetscObject)label->points[v]));
    549:     (*labelnew)->points[v] = label->points[v];
    550:   }
    551:   PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
    552:   PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: /*@
    557:   DMLabelDuplicate - Duplicates a `DMLabel`

    559:   Collective

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

    564:   Output Parameter:
    565: . labelnew - new label

    567:   Level: intermediate

    569: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
    570: @*/
    571: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
    572: {
    573:   const char *name;

    575:   PetscFunctionBegin;
    577:   PetscCall(DMLabelMakeAllValid_Private(label));
    578:   PetscCall(PetscObjectGetName((PetscObject)label, &name));
    579:   PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));

    581:   (*labelnew)->numStrata    = label->numStrata;
    582:   (*labelnew)->defaultValue = label->defaultValue;
    583:   (*labelnew)->readonly     = label->readonly;
    584:   PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
    585:   PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
    586:   PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
    587:   PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
    588:   PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
    589:   for (PetscInt v = 0; v < label->numStrata; ++v) {
    590:     (*labelnew)->stratumValues[v] = label->stratumValues[v];
    591:     (*labelnew)->stratumSizes[v]  = label->stratumSizes[v];
    592:     (*labelnew)->validIS[v]       = PETSC_TRUE;
    593:   }
    594:   (*labelnew)->pStart = -1;
    595:   (*labelnew)->pEnd   = -1;
    596:   (*labelnew)->bt     = NULL;
    597:   PetscUseTypeMethod(label, duplicate, labelnew);
    598:   PetscFunctionReturn(PETSC_SUCCESS);
    599: }

    601: /*@C
    602:   DMLabelCompare - Compare two `DMLabel` objects

    604:   Collective; No Fortran Support

    606:   Input Parameters:
    607: + comm - Comm over which to compare labels
    608: . l0   - First `DMLabel`
    609: - l1   - Second `DMLabel`

    611:   Output Parameters:
    612: + equal   - (Optional) Flag whether the two labels are equal
    613: - message - (Optional) Message describing the difference

    615:   Level: intermediate

    617:   Notes:
    618:   The output flag equal is the same on all processes.
    619:   If it is passed as `NULL` and difference is found, an error is thrown on all processes.
    620:   Make sure to pass `NULL` on all processes.

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

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

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

    631:   Developer Note:
    632:   Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()`

    634: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
    635: @*/
    636: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
    637: {
    638:   const char *name0, *name1;
    639:   char        msg[PETSC_MAX_PATH_LEN] = "";
    640:   PetscBool   eq;
    641:   PetscMPIInt rank;

    643:   PetscFunctionBegin;
    646:   if (equal) PetscAssertPointer(equal, 4);
    647:   if (message) PetscAssertPointer(message, 5);
    648:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    649:   PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
    650:   PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
    651:   {
    652:     PetscInt v0, v1;

    654:     PetscCall(DMLabelGetDefaultValue(l0, &v0));
    655:     PetscCall(DMLabelGetDefaultValue(l1, &v1));
    656:     eq = (PetscBool)(v0 == v1);
    657:     if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %" PetscInt_FMT " != %" PetscInt_FMT " = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1));
    658:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
    659:     if (!eq) goto finish;
    660:   }
    661:   {
    662:     IS is0, is1;

    664:     PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
    665:     PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
    666:     PetscCall(ISEqual(is0, is1, &eq));
    667:     PetscCall(ISDestroy(&is0));
    668:     PetscCall(ISDestroy(&is1));
    669:     if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
    670:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
    671:     if (!eq) goto finish;
    672:   }
    673:   {
    674:     PetscInt i, nValues;

    676:     PetscCall(DMLabelGetNumValues(l0, &nValues));
    677:     for (i = 0; i < nValues; i++) {
    678:       const PetscInt v = l0->stratumValues[i];
    679:       PetscInt       n;
    680:       IS             is0, is1;

    682:       PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
    683:       if (!n) continue;
    684:       PetscCall(DMLabelGetStratumIS(l0, v, &is0));
    685:       PetscCall(DMLabelGetStratumIS(l1, v, &is1));
    686:       PetscCall(ISEqualUnsorted(is0, is1, &eq));
    687:       PetscCall(ISDestroy(&is0));
    688:       PetscCall(ISDestroy(&is1));
    689:       if (!eq) {
    690:         PetscCall(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));
    691:         break;
    692:       }
    693:     }
    694:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
    695:   }
    696: finish:
    697:   /* If message output arg not set, print to stderr */
    698:   if (message) {
    699:     *message = NULL;
    700:     if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
    701:   } else {
    702:     if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
    703:     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
    704:   }
    705:   /* If same output arg not ser and labels are not equal, throw error */
    706:   if (equal) *equal = eq;
    707:   else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
    708:   PetscFunctionReturn(PETSC_SUCCESS);
    709: }

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

    714:   Not Collective

    716:   Input Parameter:
    717: . label - The `DMLabel`

    719:   Level: intermediate

    721: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    722: @*/
    723: PetscErrorCode DMLabelComputeIndex(DMLabel label)
    724: {
    725:   PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v;

    727:   PetscFunctionBegin;
    729:   PetscCall(DMLabelMakeAllValid_Private(label));
    730:   for (v = 0; v < label->numStrata; ++v) {
    731:     const PetscInt *points;
    732:     PetscInt        i;

    734:     PetscCall(ISGetIndices(label->points[v], &points));
    735:     for (i = 0; i < label->stratumSizes[v]; ++i) {
    736:       const PetscInt point = points[i];

    738:       pStart = PetscMin(point, pStart);
    739:       pEnd   = PetscMax(point + 1, pEnd);
    740:     }
    741:     PetscCall(ISRestoreIndices(label->points[v], &points));
    742:   }
    743:   label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
    744:   label->pEnd   = pEnd;
    745:   PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
    746:   PetscFunctionReturn(PETSC_SUCCESS);
    747: }

    749: /*@
    750:   DMLabelCreateIndex - Create an index structure for membership determination

    752:   Not Collective

    754:   Input Parameters:
    755: + label  - The `DMLabel`
    756: . pStart - The smallest point
    757: - pEnd   - The largest point + 1

    759:   Level: intermediate

    761: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    762: @*/
    763: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
    764: {
    765:   PetscInt v;

    767:   PetscFunctionBegin;
    769:   PetscCall(DMLabelDestroyIndex(label));
    770:   PetscCall(DMLabelMakeAllValid_Private(label));
    771:   label->pStart = pStart;
    772:   label->pEnd   = pEnd;
    773:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
    774:   PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
    775:   for (v = 0; v < label->numStrata; ++v) {
    776:     IS              pointIS;
    777:     const PetscInt *points;
    778:     PetscInt        i;

    780:     PetscUseTypeMethod(label, getstratumis, v, &pointIS);
    781:     PetscCall(ISGetIndices(pointIS, &points));
    782:     for (i = 0; i < label->stratumSizes[v]; ++i) {
    783:       const PetscInt point = points[i];

    785:       PetscCheck(!(point < pStart) && !(point >= pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " in stratum %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->stratumValues[v], pStart, pEnd);
    786:       PetscCall(PetscBTSet(label->bt, point - pStart));
    787:     }
    788:     PetscCall(ISRestoreIndices(label->points[v], &points));
    789:     PetscCall(ISDestroy(&pointIS));
    790:   }
    791:   PetscFunctionReturn(PETSC_SUCCESS);
    792: }

    794: /*@
    795:   DMLabelDestroyIndex - Destroy the index structure

    797:   Not Collective

    799:   Input Parameter:
    800: . label - the `DMLabel`

    802:   Level: intermediate

    804: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    805: @*/
    806: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
    807: {
    808:   PetscFunctionBegin;
    810:   label->pStart = -1;
    811:   label->pEnd   = -1;
    812:   PetscCall(PetscBTDestroy(&label->bt));
    813:   PetscFunctionReturn(PETSC_SUCCESS);
    814: }

    816: /*@
    817:   DMLabelGetBounds - Return the smallest and largest point in the label

    819:   Not Collective

    821:   Input Parameter:
    822: . label - the `DMLabel`

    824:   Output Parameters:
    825: + pStart - The smallest point
    826: - pEnd   - The largest point + 1

    828:   Level: intermediate

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

    833: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    834: @*/
    835: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
    836: {
    837:   PetscFunctionBegin;
    839:   if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
    840:   if (pStart) {
    841:     PetscAssertPointer(pStart, 2);
    842:     *pStart = label->pStart;
    843:   }
    844:   if (pEnd) {
    845:     PetscAssertPointer(pEnd, 3);
    846:     *pEnd = label->pEnd;
    847:   }
    848:   PetscFunctionReturn(PETSC_SUCCESS);
    849: }

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

    854:   Not Collective

    856:   Input Parameters:
    857: + label - the `DMLabel`
    858: - value - the value

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

    863:   Level: developer

    865: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    866: @*/
    867: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
    868: {
    869:   PetscInt v;

    871:   PetscFunctionBegin;
    873:   PetscAssertPointer(contains, 3);
    874:   PetscCall(DMLabelLookupStratum(label, value, &v));
    875:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
    876:   PetscFunctionReturn(PETSC_SUCCESS);
    877: }

    879: /*@
    880:   DMLabelHasPoint - Determine whether a label assigns a value to a point

    882:   Not Collective

    884:   Input Parameters:
    885: + label - the `DMLabel`
    886: - point - the point

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

    891:   Level: developer

    893:   Note:
    894:   The user must call `DMLabelCreateIndex()` before this function.

    896: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    897: @*/
    898: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
    899: {
    900:   PetscFunctionBeginHot;
    902:   PetscAssertPointer(contains, 3);
    903:   PetscCall(DMLabelMakeAllValid_Private(label));
    904:   if (PetscDefined(USE_DEBUG)) {
    905:     PetscCheck(label->bt, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
    906:     PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
    907:   }
    908:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
    909:   PetscFunctionReturn(PETSC_SUCCESS);
    910: }

    912: /*@
    913:   DMLabelStratumHasPoint - Return true if the stratum contains a point

    915:   Not Collective

    917:   Input Parameters:
    918: + label - the `DMLabel`
    919: . value - the stratum value
    920: - point - the point

    922:   Output Parameter:
    923: . contains - true if the stratum contains the point

    925:   Level: intermediate

    927: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    928: @*/
    929: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
    930: {
    931:   PetscFunctionBeginHot;
    933:   PetscAssertPointer(contains, 4);
    934:   if (value == label->defaultValue) {
    935:     PetscInt pointVal;

    937:     PetscCall(DMLabelGetValue(label, point, &pointVal));
    938:     *contains = (PetscBool)(pointVal == value);
    939:   } else {
    940:     PetscInt v;

    942:     PetscCall(DMLabelLookupStratum(label, value, &v));
    943:     if (v >= 0) {
    944:       if (label->validIS[v] || label->readonly) {
    945:         IS       is;
    946:         PetscInt i;

    948:         PetscUseTypeMethod(label, getstratumis, v, &is);
    949:         PetscCall(ISLocate(is, point, &i));
    950:         PetscCall(ISDestroy(&is));
    951:         *contains = (PetscBool)(i >= 0);
    952:       } else {
    953:         PetscCall(PetscHSetIHas(label->ht[v], point, contains));
    954:       }
    955:     } else { // value is not present
    956:       *contains = PETSC_FALSE;
    957:     }
    958:   }
    959:   PetscFunctionReturn(PETSC_SUCCESS);
    960: }

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

    966:   Not Collective

    968:   Input Parameter:
    969: . label - a `DMLabel` object

    971:   Output Parameter:
    972: . defaultValue - the default value

    974:   Level: beginner

    976: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    977: @*/
    978: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
    979: {
    980:   PetscFunctionBegin;
    982:   *defaultValue = label->defaultValue;
    983:   PetscFunctionReturn(PETSC_SUCCESS);
    984: }

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

    990:   Not Collective

    992:   Input Parameter:
    993: . label - a `DMLabel` object

    995:   Output Parameter:
    996: . defaultValue - the default value

    998:   Level: beginner

    1000: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    1001: @*/
    1002: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
    1003: {
    1004:   PetscFunctionBegin;
    1006:   label->defaultValue = defaultValue;
    1007:   PetscFunctionReturn(PETSC_SUCCESS);
    1008: }

    1010: /*@
    1011:   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
    1012:   `DMLabelSetDefaultValue()`)

    1014:   Not Collective

    1016:   Input Parameters:
    1017: + label - the `DMLabel`
    1018: - point - the point

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

    1023:   Level: intermediate

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

    1029: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
    1030: @*/
    1031: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
    1032: {
    1033:   PetscInt v;

    1035:   PetscFunctionBeginHot;
    1037:   PetscAssertPointer(value, 3);
    1038:   *value = label->defaultValue;
    1039:   for (v = 0; v < label->numStrata; ++v) {
    1040:     if (label->validIS[v] || label->readonly) {
    1041:       IS       is;
    1042:       PetscInt i;

    1044:       PetscUseTypeMethod(label, getstratumis, v, &is);
    1045:       PetscCall(ISLocate(label->points[v], point, &i));
    1046:       PetscCall(ISDestroy(&is));
    1047:       if (i >= 0) {
    1048:         *value = label->stratumValues[v];
    1049:         break;
    1050:       }
    1051:     } else {
    1052:       PetscBool has;

    1054:       PetscCall(PetscHSetIHas(label->ht[v], point, &has));
    1055:       if (has) {
    1056:         *value = label->stratumValues[v];
    1057:         break;
    1058:       }
    1059:     }
    1060:   }
    1061:   PetscFunctionReturn(PETSC_SUCCESS);
    1062: }

    1064: /*@
    1065:   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
    1066:   be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.

    1068:   Not Collective

    1070:   Input Parameters:
    1071: + label - the `DMLabel`
    1072: . point - the point
    1073: - value - The point value

    1075:   Level: intermediate

    1077: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
    1078: @*/
    1079: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
    1080: {
    1081:   PetscInt v;

    1083:   PetscFunctionBegin;
    1085:   /* Find label value, add new entry if needed */
    1086:   if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
    1087:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1088:   PetscCall(DMLabelLookupAddStratum(label, value, &v));
    1089:   /* Set key */
    1090:   PetscCall(DMLabelMakeInvalid_Private(label, v));
    1091:   PetscCall(PetscHSetIAdd(label->ht[v], point));
    1092:   PetscFunctionReturn(PETSC_SUCCESS);
    1093: }

    1095: /*@
    1096:   DMLabelClearValue - Clear the value a label assigns to a point

    1098:   Not Collective

    1100:   Input Parameters:
    1101: + label - the `DMLabel`
    1102: . point - the point
    1103: - value - The point value

    1105:   Level: intermediate

    1107: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    1108: @*/
    1109: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
    1110: {
    1111:   PetscInt v;

    1113:   PetscFunctionBegin;
    1115:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1116:   /* Find label value */
    1117:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1118:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);

    1120:   if (label->bt) {
    1121:     PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
    1122:     PetscCall(PetscBTClear(label->bt, point - label->pStart));
    1123:   }

    1125:   /* Delete key */
    1126:   PetscCall(DMLabelMakeInvalid_Private(label, v));
    1127:   PetscCall(PetscHSetIDel(label->ht[v], point));
    1128:   PetscFunctionReturn(PETSC_SUCCESS);
    1129: }

    1131: /*@
    1132:   DMLabelInsertIS - Set all points in the `IS` to a value

    1134:   Not Collective

    1136:   Input Parameters:
    1137: + label - the `DMLabel`
    1138: . is    - the point `IS`
    1139: - value - The point value

    1141:   Level: intermediate

    1143: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1144: @*/
    1145: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
    1146: {
    1147:   PetscInt        v, n, p;
    1148:   const PetscInt *points;

    1150:   PetscFunctionBegin;
    1153:   /* Find label value, add new entry if needed */
    1154:   if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
    1155:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1156:   PetscCall(DMLabelLookupAddStratum(label, value, &v));
    1157:   /* Set keys */
    1158:   PetscCall(DMLabelMakeInvalid_Private(label, v));
    1159:   PetscCall(ISGetLocalSize(is, &n));
    1160:   PetscCall(ISGetIndices(is, &points));
    1161:   for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
    1162:   PetscCall(ISRestoreIndices(is, &points));
    1163:   PetscFunctionReturn(PETSC_SUCCESS);
    1164: }

    1166: /*@
    1167:   DMLabelGetNumValues - Get the number of values that the `DMLabel` takes

    1169:   Not Collective

    1171:   Input Parameter:
    1172: . label - the `DMLabel`

    1174:   Output Parameter:
    1175: . numValues - the number of values

    1177:   Level: intermediate

    1179: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1180: @*/
    1181: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
    1182: {
    1183:   PetscFunctionBegin;
    1185:   PetscAssertPointer(numValues, 2);
    1186:   *numValues = label->numStrata;
    1187:   PetscFunctionReturn(PETSC_SUCCESS);
    1188: }

    1190: /*@
    1191:   DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes

    1193:   Not Collective

    1195:   Input Parameter:
    1196: . label - the `DMLabel`

    1198:   Output Parameter:
    1199: . values - the value `IS`

    1201:   Level: intermediate

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

    1208: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1209: @*/
    1210: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
    1211: {
    1212:   PetscFunctionBegin;
    1214:   PetscAssertPointer(values, 2);
    1215:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
    1216:   PetscFunctionReturn(PETSC_SUCCESS);
    1217: }

    1219: /*@
    1220:   DMLabelGetValueBounds - Return the smallest and largest value in the label

    1222:   Not Collective

    1224:   Input Parameter:
    1225: . label - the `DMLabel`

    1227:   Output Parameters:
    1228: + minValue - The smallest value
    1229: - maxValue - The largest value

    1231:   Level: intermediate

    1233: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
    1234: @*/
    1235: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
    1236: {
    1237:   PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;

    1239:   PetscFunctionBegin;
    1241:   for (PetscInt v = 0; v < label->numStrata; ++v) {
    1242:     min = PetscMin(min, label->stratumValues[v]);
    1243:     max = PetscMax(max, label->stratumValues[v]);
    1244:   }
    1245:   if (minValue) {
    1246:     PetscAssertPointer(minValue, 2);
    1247:     *minValue = min;
    1248:   }
    1249:   if (maxValue) {
    1250:     PetscAssertPointer(maxValue, 3);
    1251:     *maxValue = max;
    1252:   }
    1253:   PetscFunctionReturn(PETSC_SUCCESS);
    1254: }

    1256: /*@
    1257:   DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes

    1259:   Not Collective

    1261:   Input Parameter:
    1262: . label - the `DMLabel`

    1264:   Output Parameter:
    1265: . values - the value `IS`

    1267:   Level: intermediate

    1269:   Notes:
    1270:   The output `IS` should be destroyed when no longer needed.
    1271:   This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.

    1273: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1274: @*/
    1275: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
    1276: {
    1277:   PetscInt  i, j;
    1278:   PetscInt *valuesArr;

    1280:   PetscFunctionBegin;
    1282:   PetscAssertPointer(values, 2);
    1283:   PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
    1284:   for (i = 0, j = 0; i < label->numStrata; i++) {
    1285:     PetscInt n;

    1287:     PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
    1288:     if (n) valuesArr[j++] = label->stratumValues[i];
    1289:   }
    1290:   if (j == label->numStrata) {
    1291:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
    1292:   } else {
    1293:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
    1294:   }
    1295:   PetscCall(PetscFree(valuesArr));
    1296:   PetscFunctionReturn(PETSC_SUCCESS);
    1297: }

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

    1302:   Not Collective

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

    1308:   Output Parameter:
    1309: . index - the index of value in the list of values

    1311:   Level: intermediate

    1313: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1314: @*/
    1315: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
    1316: {
    1317:   PetscInt v;

    1319:   PetscFunctionBegin;
    1321:   PetscAssertPointer(index, 3);
    1322:   /* Do not assume they are sorted */
    1323:   for (v = 0; v < label->numStrata; ++v)
    1324:     if (label->stratumValues[v] == value) break;
    1325:   if (v >= label->numStrata) *index = -1;
    1326:   else *index = v;
    1327:   PetscFunctionReturn(PETSC_SUCCESS);
    1328: }

    1330: /*@
    1331:   DMLabelHasStratum - Determine whether points exist with the given value

    1333:   Not Collective

    1335:   Input Parameters:
    1336: + label - the `DMLabel`
    1337: - value - the stratum value

    1339:   Output Parameter:
    1340: . exists - Flag saying whether points exist

    1342:   Level: intermediate

    1344: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1345: @*/
    1346: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
    1347: {
    1348:   PetscInt v;

    1350:   PetscFunctionBegin;
    1352:   PetscAssertPointer(exists, 3);
    1353:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1354:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
    1355:   PetscFunctionReturn(PETSC_SUCCESS);
    1356: }

    1358: /*@
    1359:   DMLabelGetStratumSize - Get the size of a stratum

    1361:   Not Collective

    1363:   Input Parameters:
    1364: + label - the `DMLabel`
    1365: - value - the stratum value

    1367:   Output Parameter:
    1368: . size - The number of points in the stratum

    1370:   Level: intermediate

    1372: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1373: @*/
    1374: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
    1375: {
    1376:   PetscInt v;

    1378:   PetscFunctionBegin;
    1380:   PetscAssertPointer(size, 3);
    1381:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1382:   PetscCall(DMLabelGetStratumSize_Private(label, v, size));
    1383:   PetscFunctionReturn(PETSC_SUCCESS);
    1384: }

    1386: /*@
    1387:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

    1389:   Not Collective

    1391:   Input Parameters:
    1392: + label - the `DMLabel`
    1393: - value - the stratum value

    1395:   Output Parameters:
    1396: + start - the smallest point in the stratum
    1397: - end   - the largest point in the stratum

    1399:   Level: intermediate

    1401: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1402: @*/
    1403: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
    1404: {
    1405:   IS       is;
    1406:   PetscInt v, min, max;

    1408:   PetscFunctionBegin;
    1410:   if (start) {
    1411:     PetscAssertPointer(start, 3);
    1412:     *start = -1;
    1413:   }
    1414:   if (end) {
    1415:     PetscAssertPointer(end, 4);
    1416:     *end = -1;
    1417:   }
    1418:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1419:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
    1420:   PetscCall(DMLabelMakeValid_Private(label, v));
    1421:   if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
    1422:   PetscUseTypeMethod(label, getstratumis, v, &is);
    1423:   PetscCall(ISGetMinMax(is, &min, &max));
    1424:   PetscCall(ISDestroy(&is));
    1425:   if (start) *start = min;
    1426:   if (end) *end = max + 1;
    1427:   PetscFunctionReturn(PETSC_SUCCESS);
    1428: }

    1430: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
    1431: {
    1432:   PetscFunctionBegin;
    1433:   PetscCall(PetscObjectReference((PetscObject)label->points[v]));
    1434:   *pointIS = label->points[v];
    1435:   PetscFunctionReturn(PETSC_SUCCESS);
    1436: }

    1438: /*@
    1439:   DMLabelGetStratumIS - Get an `IS` with the stratum points

    1441:   Not Collective

    1443:   Input Parameters:
    1444: + label - the `DMLabel`
    1445: - value - the stratum value

    1447:   Output Parameter:
    1448: . points - The stratum points

    1450:   Level: intermediate

    1452:   Notes:
    1453:   The output `IS` should be destroyed when no longer needed.
    1454:   Returns `NULL` if the stratum is empty.

    1456: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1457: @*/
    1458: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
    1459: {
    1460:   PetscInt v;

    1462:   PetscFunctionBegin;
    1464:   PetscAssertPointer(points, 3);
    1465:   *points = NULL;
    1466:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1467:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
    1468:   PetscCall(DMLabelMakeValid_Private(label, v));
    1469:   PetscUseTypeMethod(label, getstratumis, v, points);
    1470:   PetscFunctionReturn(PETSC_SUCCESS);
    1471: }

    1473: /*@
    1474:   DMLabelSetStratumIS - Set the stratum points using an `IS`

    1476:   Not Collective

    1478:   Input Parameters:
    1479: + label - the `DMLabel`
    1480: . value - the stratum value
    1481: - is    - The stratum points

    1483:   Level: intermediate

    1485: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1486: @*/
    1487: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
    1488: {
    1489:   PetscInt v;

    1491:   PetscFunctionBegin;
    1494:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1495:   PetscCall(DMLabelLookupAddStratum(label, value, &v));
    1496:   if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
    1497:   PetscCall(DMLabelClearStratum(label, value));
    1498:   PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
    1499:   PetscCall(PetscObjectReference((PetscObject)is));
    1500:   PetscCall(ISDestroy(&label->points[v]));
    1501:   label->points[v]  = is;
    1502:   label->validIS[v] = PETSC_TRUE;
    1503:   PetscCall(PetscObjectStateIncrease((PetscObject)label));
    1504:   if (label->bt) {
    1505:     const PetscInt *points;
    1506:     PetscInt        p;

    1508:     PetscCall(ISGetIndices(is, &points));
    1509:     for (p = 0; p < label->stratumSizes[v]; ++p) {
    1510:       const PetscInt point = points[p];

    1512:       PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
    1513:       PetscCall(PetscBTSet(label->bt, point - label->pStart));
    1514:     }
    1515:   }
    1516:   PetscFunctionReturn(PETSC_SUCCESS);
    1517: }

    1519: /*@
    1520:   DMLabelClearStratum - Remove a stratum

    1522:   Not Collective

    1524:   Input Parameters:
    1525: + label - the `DMLabel`
    1526: - value - the stratum value

    1528:   Level: intermediate

    1530: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1531: @*/
    1532: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
    1533: {
    1534:   PetscInt v;

    1536:   PetscFunctionBegin;
    1538:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1539:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1540:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
    1541:   if (label->validIS[v]) {
    1542:     if (label->bt) {
    1543:       PetscInt        i;
    1544:       const PetscInt *points;

    1546:       PetscCall(ISGetIndices(label->points[v], &points));
    1547:       for (i = 0; i < label->stratumSizes[v]; ++i) {
    1548:         const PetscInt point = points[i];

    1550:         PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
    1551:         PetscCall(PetscBTClear(label->bt, point - label->pStart));
    1552:       }
    1553:       PetscCall(ISRestoreIndices(label->points[v], &points));
    1554:     }
    1555:     label->stratumSizes[v] = 0;
    1556:     PetscCall(ISDestroy(&label->points[v]));
    1557:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
    1558:     PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
    1559:     PetscCall(PetscObjectStateIncrease((PetscObject)label));
    1560:   } else {
    1561:     PetscCall(PetscHSetIClear(label->ht[v]));
    1562:   }
    1563:   PetscFunctionReturn(PETSC_SUCCESS);
    1564: }

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

    1569:   Not Collective

    1571:   Input Parameters:
    1572: + label  - The `DMLabel`
    1573: . value  - The label value for all points
    1574: . pStart - The first point
    1575: - pEnd   - A point beyond all marked points

    1577:   Level: intermediate

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

    1582: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
    1583: @*/
    1584: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
    1585: {
    1586:   IS pIS;

    1588:   PetscFunctionBegin;
    1589:   PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
    1590:   PetscCall(DMLabelSetStratumIS(label, value, pIS));
    1591:   PetscCall(ISDestroy(&pIS));
    1592:   PetscFunctionReturn(PETSC_SUCCESS);
    1593: }

    1595: /*@
    1596:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

    1598:   Not Collective

    1600:   Input Parameters:
    1601: + label - The `DMLabel`
    1602: . value - The label value
    1603: - p     - A point with this value

    1605:   Output Parameter:
    1606: . 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

    1608:   Level: intermediate

    1610: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
    1611: @*/
    1612: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
    1613: {
    1614:   IS       pointIS;
    1615:   PetscInt v;

    1617:   PetscFunctionBegin;
    1619:   PetscAssertPointer(index, 4);
    1620:   *index = -1;
    1621:   PetscCall(DMLabelLookupStratum(label, value, &v));
    1622:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
    1623:   PetscCall(DMLabelMakeValid_Private(label, v));
    1624:   PetscUseTypeMethod(label, getstratumis, v, &pointIS);
    1625:   PetscCall(ISLocate(pointIS, p, index));
    1626:   PetscCall(ISDestroy(&pointIS));
    1627:   PetscFunctionReturn(PETSC_SUCCESS);
    1628: }

    1630: /*@
    1631:   DMLabelFilter - Remove all points outside of [`start`, `end`)

    1633:   Not Collective

    1635:   Input Parameters:
    1636: + label - the `DMLabel`
    1637: . start - the first point kept
    1638: - end   - one more than the last point kept

    1640:   Level: intermediate

    1642: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1643: @*/
    1644: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
    1645: {
    1646:   PetscInt v;

    1648:   PetscFunctionBegin;
    1650:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1651:   PetscCall(DMLabelDestroyIndex(label));
    1652:   PetscCall(DMLabelMakeAllValid_Private(label));
    1653:   for (v = 0; v < label->numStrata; ++v) {
    1654:     PetscCall(ISGeneralFilter(label->points[v], start, end));
    1655:     PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
    1656:   }
    1657:   PetscCall(DMLabelCreateIndex(label, start, end));
    1658:   PetscFunctionReturn(PETSC_SUCCESS);
    1659: }

    1661: /*@
    1662:   DMLabelPermute - Create a new label with permuted points

    1664:   Not Collective

    1666:   Input Parameters:
    1667: + label       - the `DMLabel`
    1668: - permutation - the point permutation

    1670:   Output Parameter:
    1671: . labelNew - the new label containing the permuted points

    1673:   Level: intermediate

    1675: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1676: @*/
    1677: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
    1678: {
    1679:   const PetscInt *perm;
    1680:   PetscInt        numValues, numPoints, v, q;

    1682:   PetscFunctionBegin;
    1685:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1686:   PetscCall(DMLabelMakeAllValid_Private(label));
    1687:   PetscCall(DMLabelDuplicate(label, labelNew));
    1688:   PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
    1689:   PetscCall(ISGetLocalSize(permutation, &numPoints));
    1690:   PetscCall(ISGetIndices(permutation, &perm));
    1691:   for (v = 0; v < numValues; ++v) {
    1692:     const PetscInt  size = (*labelNew)->stratumSizes[v];
    1693:     const PetscInt *points;
    1694:     PetscInt       *pointsNew;

    1696:     PetscCall(ISGetIndices((*labelNew)->points[v], &points));
    1697:     PetscCall(PetscCalloc1(size, &pointsNew));
    1698:     for (q = 0; q < size; ++q) {
    1699:       const PetscInt point = points[q];

    1701:       PetscCheck(!(point < 0) && !(point >= numPoints), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ") for the remapping", point, numPoints);
    1702:       pointsNew[q] = perm[point];
    1703:     }
    1704:     PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
    1705:     PetscCall(PetscSortInt(size, pointsNew));
    1706:     PetscCall(ISDestroy(&(*labelNew)->points[v]));
    1707:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
    1708:       PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
    1709:       PetscCall(PetscFree(pointsNew));
    1710:     } else {
    1711:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
    1712:     }
    1713:     PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
    1714:   }
    1715:   PetscCall(ISRestoreIndices(permutation, &perm));
    1716:   if (label->bt) {
    1717:     PetscCall(PetscBTDestroy(&label->bt));
    1718:     PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
    1719:   }
    1720:   PetscFunctionReturn(PETSC_SUCCESS);
    1721: }

    1723: /*@
    1724:   DMLabelPermuteValues - Permute the values in a label

    1726:   Not collective

    1728:   Input Parameters:
    1729: + label       - the `DMLabel`
    1730: - permutation - the value permutation, permutation[old value] = new value

    1732:   Output Parameter:
    1733: . label - the `DMLabel` now with permuted values

    1735:   Note:
    1736:   The modification is done in-place

    1738:   Level: intermediate

    1740: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1741: @*/
    1742: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
    1743: {
    1744:   PetscInt Nv, Np;

    1746:   PetscFunctionBegin;
    1749:   PetscCall(DMLabelGetNumValues(label, &Nv));
    1750:   PetscCall(ISGetLocalSize(permutation, &Np));
    1751:   PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
    1752:   if (PetscDefined(USE_DEBUG)) {
    1753:     PetscBool flg;
    1754:     PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
    1755:     PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
    1756:   }
    1757:   PetscCall(DMLabelRewriteValues(label, permutation));
    1758:   PetscFunctionReturn(PETSC_SUCCESS);
    1759: }

    1761: /*@
    1762:   DMLabelRewriteValues - Permute the values in a label, but some may be omitted

    1764:   Not collective

    1766:   Input Parameters:
    1767: + label       - the `DMLabel`
    1768: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted

    1770:   Output Parameter:
    1771: . label - the `DMLabel` now with permuted values

    1773:   Note:
    1774:   The modification is done in-place

    1776:   Level: intermediate

    1778: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1779: @*/
    1780: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
    1781: {
    1782:   const PetscInt *perm;
    1783:   PetscInt        Nv, Np;

    1785:   PetscFunctionBegin;
    1788:   PetscCall(DMLabelMakeAllValid_Private(label));
    1789:   PetscCall(DMLabelGetNumValues(label, &Nv));
    1790:   PetscCall(ISGetLocalSize(permutation, &Np));
    1791:   PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
    1792:   PetscCall(ISGetIndices(permutation, &perm));
    1793:   for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
    1794:   PetscCall(ISRestoreIndices(permutation, &perm));
    1795:   PetscFunctionReturn(PETSC_SUCCESS);
    1796: }

    1798: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
    1799: {
    1800:   MPI_Comm     comm;
    1801:   PetscInt     s, l, nroots, nleaves, offset, size;
    1802:   PetscInt    *remoteOffsets, *rootStrata, *rootIdx;
    1803:   PetscSection rootSection;
    1804:   PetscSF      labelSF;

    1806:   PetscFunctionBegin;
    1807:   if (label) PetscCall(DMLabelMakeAllValid_Private(label));
    1808:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
    1809:   /* Build a section of stratum values per point, generate the according SF
    1810:      and distribute point-wise stratum values to leaves. */
    1811:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
    1812:   PetscCall(PetscSectionCreate(comm, &rootSection));
    1813:   PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
    1814:   if (label) {
    1815:     for (s = 0; s < label->numStrata; ++s) {
    1816:       const PetscInt *points;

    1818:       PetscCall(ISGetIndices(label->points[s], &points));
    1819:       for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
    1820:       PetscCall(ISRestoreIndices(label->points[s], &points));
    1821:     }
    1822:   }
    1823:   PetscCall(PetscSectionSetUp(rootSection));
    1824:   /* Create a point-wise array of stratum values */
    1825:   PetscCall(PetscSectionGetStorageSize(rootSection, &size));
    1826:   PetscCall(PetscMalloc1(size, &rootStrata));
    1827:   PetscCall(PetscCalloc1(nroots, &rootIdx));
    1828:   if (label) {
    1829:     for (s = 0; s < label->numStrata; ++s) {
    1830:       const PetscInt *points;

    1832:       PetscCall(ISGetIndices(label->points[s], &points));
    1833:       for (l = 0; l < label->stratumSizes[s]; l++) {
    1834:         const PetscInt p = points[l];
    1835:         PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
    1836:         rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
    1837:       }
    1838:       PetscCall(ISRestoreIndices(label->points[s], &points));
    1839:     }
    1840:   }
    1841:   /* Build SF that maps label points to remote processes */
    1842:   PetscCall(PetscSectionCreate(comm, leafSection));
    1843:   PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
    1844:   PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
    1845:   PetscCall(PetscFree(remoteOffsets));
    1846:   /* Send the strata for each point over the derived SF */
    1847:   PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
    1848:   PetscCall(PetscMalloc1(size, leafStrata));
    1849:   PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
    1850:   PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
    1851:   /* Clean up */
    1852:   PetscCall(PetscFree(rootStrata));
    1853:   PetscCall(PetscFree(rootIdx));
    1854:   PetscCall(PetscSectionDestroy(&rootSection));
    1855:   PetscCall(PetscSFDestroy(&labelSF));
    1856:   PetscFunctionReturn(PETSC_SUCCESS);
    1857: }

    1859: /*@
    1860:   DMLabelDistribute - Create a new label pushed forward over the `PetscSF`

    1862:   Collective

    1864:   Input Parameters:
    1865: + label - the `DMLabel`
    1866: - sf    - the map from old to new distribution

    1868:   Output Parameter:
    1869: . labelNew - the new redistributed label

    1871:   Level: intermediate

    1873: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
    1874: @*/
    1875: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
    1876: {
    1877:   MPI_Comm     comm;
    1878:   PetscSection leafSection;
    1879:   PetscInt     p, pStart, pEnd, s, size, dof, offset, stratum;
    1880:   PetscInt    *leafStrata, *strataIdx;
    1881:   PetscInt   **points;
    1882:   const char  *lname = NULL;
    1883:   char        *name;
    1884:   PetscMPIInt  nameSize;
    1885:   PetscHSetI   stratumHash;
    1886:   size_t       len = 0;
    1887:   PetscMPIInt  rank;

    1889:   PetscFunctionBegin;
    1891:   if (label) {
    1893:     PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    1894:     PetscCall(DMLabelMakeAllValid_Private(label));
    1895:   }
    1896:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
    1897:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1898:   /* Bcast name */
    1899:   if (rank == 0) {
    1900:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
    1901:     PetscCall(PetscStrlen(lname, &len));
    1902:   }
    1903:   PetscCall(PetscMPIIntCast(len, &nameSize));
    1904:   PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
    1905:   PetscCall(PetscMalloc1(nameSize + 1, &name));
    1906:   if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
    1907:   PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
    1908:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
    1909:   PetscCall(PetscFree(name));
    1910:   /* Bcast defaultValue */
    1911:   if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
    1912:   PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
    1913:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
    1914:   PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
    1915:   /* Determine received stratum values and initialise new label*/
    1916:   PetscCall(PetscHSetICreate(&stratumHash));
    1917:   PetscCall(PetscSectionGetStorageSize(leafSection, &size));
    1918:   for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
    1919:   PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
    1920:   PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
    1921:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
    1922:   PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
    1923:   /* Turn leafStrata into indices rather than stratum values */
    1924:   offset = 0;
    1925:   PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
    1926:   PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
    1927:   for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
    1928:   for (p = 0; p < size; ++p) {
    1929:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
    1930:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
    1931:         leafStrata[p] = s;
    1932:         break;
    1933:       }
    1934:     }
    1935:   }
    1936:   /* Rebuild the point strata on the receiver */
    1937:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
    1938:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
    1939:   for (p = pStart; p < pEnd; p++) {
    1940:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
    1941:     PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
    1942:     for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
    1943:   }
    1944:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
    1945:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
    1946:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
    1947:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
    1948:     PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
    1949:     PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
    1950:   }
    1951:   /* Insert points into new strata */
    1952:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
    1953:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
    1954:   for (p = pStart; p < pEnd; p++) {
    1955:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
    1956:     PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
    1957:     for (s = 0; s < dof; s++) {
    1958:       stratum                               = leafStrata[offset + s];
    1959:       points[stratum][strataIdx[stratum]++] = p;
    1960:     }
    1961:   }
    1962:   for (s = 0; s < (*labelNew)->numStrata; s++) {
    1963:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
    1964:     PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
    1965:   }
    1966:   PetscCall(PetscFree(points));
    1967:   PetscCall(PetscHSetIDestroy(&stratumHash));
    1968:   PetscCall(PetscFree(leafStrata));
    1969:   PetscCall(PetscFree(strataIdx));
    1970:   PetscCall(PetscSectionDestroy(&leafSection));
    1971:   PetscFunctionReturn(PETSC_SUCCESS);
    1972: }

    1974: /*@
    1975:   DMLabelGather - Gather all label values from leafs into roots

    1977:   Collective

    1979:   Input Parameters:
    1980: + label - the `DMLabel`
    1981: - sf    - the `PetscSF` communication map

    1983:   Output Parameter:
    1984: . labelNew - the new `DMLabel` with localised leaf values

    1986:   Level: developer

    1988:   Note:
    1989:   This is the inverse operation to `DMLabelDistribute()`.

    1991: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
    1992: @*/
    1993: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
    1994: {
    1995:   MPI_Comm        comm;
    1996:   PetscSection    rootSection;
    1997:   PetscSF         sfLabel;
    1998:   PetscSFNode    *rootPoints, *leafPoints;
    1999:   PetscInt        p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
    2000:   const PetscInt *rootDegree, *ilocal;
    2001:   PetscInt       *rootStrata;
    2002:   const char     *lname;
    2003:   char           *name;
    2004:   PetscMPIInt     nameSize;
    2005:   size_t          len = 0;
    2006:   PetscMPIInt     rank, size;

    2008:   PetscFunctionBegin;
    2011:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    2012:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
    2013:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    2014:   PetscCallMPI(MPI_Comm_size(comm, &size));
    2015:   /* Bcast name */
    2016:   if (rank == 0) {
    2017:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
    2018:     PetscCall(PetscStrlen(lname, &len));
    2019:   }
    2020:   PetscCall(PetscMPIIntCast(len, &nameSize));
    2021:   PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
    2022:   PetscCall(PetscMalloc1(nameSize + 1, &name));
    2023:   if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
    2024:   PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
    2025:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
    2026:   PetscCall(PetscFree(name));
    2027:   /* Gather rank/index pairs of leaves into local roots to build
    2028:      an inverse, multi-rooted SF. Note that this ignores local leaf
    2029:      indexing due to the use of the multiSF in PetscSFGather. */
    2030:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
    2031:   PetscCall(PetscMalloc1(nroots, &leafPoints));
    2032:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
    2033:   for (p = 0; p < nleaves; p++) {
    2034:     PetscInt ilp = ilocal ? ilocal[p] : p;

    2036:     leafPoints[ilp].index = ilp;
    2037:     leafPoints[ilp].rank  = rank;
    2038:   }
    2039:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
    2040:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
    2041:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
    2042:   PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
    2043:   PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
    2044:   PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
    2045:   PetscCall(PetscSFCreate(comm, &sfLabel));
    2046:   PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
    2047:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
    2048:   PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
    2049:   /* Rebuild the point strata on the receiver */
    2050:   for (p = 0, idx = 0; p < nroots; p++) {
    2051:     for (d = 0; d < rootDegree[p]; d++) {
    2052:       PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
    2053:       PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
    2054:       for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
    2055:     }
    2056:     idx += rootDegree[p];
    2057:   }
    2058:   PetscCall(PetscFree(leafPoints));
    2059:   PetscCall(PetscFree(rootStrata));
    2060:   PetscCall(PetscSectionDestroy(&rootSection));
    2061:   PetscCall(PetscSFDestroy(&sfLabel));
    2062:   PetscFunctionReturn(PETSC_SUCCESS);
    2063: }

    2065: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
    2066: {
    2067:   const PetscInt *degree;
    2068:   const PetscInt *points;
    2069:   PetscInt        Nr, r, Nl, l, val, defVal;

    2071:   PetscFunctionBegin;
    2072:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
    2073:   /* Add in leaves */
    2074:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
    2075:   for (l = 0; l < Nl; ++l) {
    2076:     PetscCall(DMLabelGetValue(label, points[l], &val));
    2077:     if (val != defVal) valArray[points[l]] = val;
    2078:   }
    2079:   /* Add in shared roots */
    2080:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &degree));
    2081:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &degree));
    2082:   for (r = 0; r < Nr; ++r) {
    2083:     if (degree[r]) {
    2084:       PetscCall(DMLabelGetValue(label, r, &val));
    2085:       if (val != defVal) valArray[r] = val;
    2086:     }
    2087:   }
    2088:   PetscFunctionReturn(PETSC_SUCCESS);
    2089: }

    2091: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
    2092: {
    2093:   const PetscInt *degree;
    2094:   const PetscInt *points;
    2095:   PetscInt        Nr, r, Nl, l, val, defVal;

    2097:   PetscFunctionBegin;
    2098:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
    2099:   /* Read out leaves */
    2100:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
    2101:   for (l = 0; l < Nl; ++l) {
    2102:     const PetscInt p    = points[l];
    2103:     const PetscInt cval = valArray[p];

    2105:     if (cval != defVal) {
    2106:       PetscCall(DMLabelGetValue(label, p, &val));
    2107:       if (val == defVal) {
    2108:         PetscCall(DMLabelSetValue(label, p, cval));
    2109:         if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
    2110:       }
    2111:     }
    2112:   }
    2113:   /* Read out shared roots */
    2114:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &degree));
    2115:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &degree));
    2116:   for (r = 0; r < Nr; ++r) {
    2117:     if (degree[r]) {
    2118:       const PetscInt cval = valArray[r];

    2120:       if (cval != defVal) {
    2121:         PetscCall(DMLabelGetValue(label, r, &val));
    2122:         if (val == defVal) {
    2123:           PetscCall(DMLabelSetValue(label, r, cval));
    2124:           if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
    2125:         }
    2126:       }
    2127:     }
    2128:   }
    2129:   PetscFunctionReturn(PETSC_SUCCESS);
    2130: }

    2132: /*@
    2133:   DMLabelPropagateBegin - Setup a cycle of label propagation

    2135:   Collective

    2137:   Input Parameters:
    2138: + label - The `DMLabel` to propagate across processes
    2139: - sf    - The `PetscSF` describing parallel layout of the label points

    2141:   Level: intermediate

    2143: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
    2144: @*/
    2145: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
    2146: {
    2147:   PetscInt    Nr, r, defVal;
    2148:   PetscMPIInt size;

    2150:   PetscFunctionBegin;
    2151:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    2152:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
    2153:   if (size > 1) {
    2154:     PetscCall(DMLabelGetDefaultValue(label, &defVal));
    2155:     PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
    2156:     if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
    2157:     for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
    2158:   }
    2159:   PetscFunctionReturn(PETSC_SUCCESS);
    2160: }

    2162: /*@
    2163:   DMLabelPropagateEnd - Tear down a cycle of label propagation

    2165:   Collective

    2167:   Input Parameters:
    2168: + label   - The `DMLabel` to propagate across processes
    2169: - pointSF - The `PetscSF` describing parallel layout of the label points

    2171:   Level: intermediate

    2173: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
    2174: @*/
    2175: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
    2176: {
    2177:   PetscFunctionBegin;
    2178:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    2179:   PetscCall(PetscFree(label->propArray));
    2180:   label->propArray = NULL;
    2181:   PetscFunctionReturn(PETSC_SUCCESS);
    2182: }

    2184: /*@C
    2185:   DMLabelPropagatePush - Tear down a cycle of label propagation

    2187:   Collective

    2189:   Input Parameters:
    2190: + label     - The `DMLabel` to propagate across processes
    2191: . pointSF   - The `PetscSF` describing parallel layout of the label points
    2192: . markPoint - An optional callback that is called when a point is marked, or `NULL`
    2193: - ctx       - An optional user context for the callback, or `NULL`

    2195:   Calling sequence of `markPoint`:
    2196: + label - The `DMLabel`
    2197: . p     - The point being marked
    2198: . val   - The label value for `p`
    2199: - ctx   - An optional user context

    2201:   Level: intermediate

    2203: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
    2204: @*/
    2205: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx)
    2206: {
    2207:   PetscInt   *valArray = label->propArray, Nr;
    2208:   PetscMPIInt size;

    2210:   PetscFunctionBegin;
    2211:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
    2212:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
    2213:   PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
    2214:   if (size > 1 && Nr >= 0) {
    2215:     /* Communicate marked edges
    2216:        The current implementation allocates an array the size of the number of root. We put the label values into the
    2217:        array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.

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

    2223:        In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
    2224:        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
    2225:        edge to the queue.
    2226:     */
    2227:     PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
    2228:     PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
    2229:     PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
    2230:     PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
    2231:     PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
    2232:     PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
    2233:   }
    2234:   PetscFunctionReturn(PETSC_SUCCESS);
    2235: }

    2237: /*@
    2238:   DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label

    2240:   Not Collective

    2242:   Input Parameter:
    2243: . label - the `DMLabel`

    2245:   Output Parameters:
    2246: + section - the section giving offsets for each stratum
    2247: - is      - An `IS` containing all the label points

    2249:   Level: developer

    2251: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
    2252: @*/
    2253: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
    2254: {
    2255:   IS              vIS;
    2256:   const PetscInt *values;
    2257:   PetscInt       *points;
    2258:   PetscInt        nV, vS = 0, vE = 0, v, N;

    2260:   PetscFunctionBegin;
    2262:   PetscCall(DMLabelGetNumValues(label, &nV));
    2263:   PetscCall(DMLabelGetValueIS(label, &vIS));
    2264:   PetscCall(ISGetIndices(vIS, &values));
    2265:   if (nV) {
    2266:     vS = values[0];
    2267:     vE = values[0] + 1;
    2268:   }
    2269:   for (v = 1; v < nV; ++v) {
    2270:     vS = PetscMin(vS, values[v]);
    2271:     vE = PetscMax(vE, values[v] + 1);
    2272:   }
    2273:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
    2274:   PetscCall(PetscSectionSetChart(*section, vS, vE));
    2275:   for (v = 0; v < nV; ++v) {
    2276:     PetscInt n;

    2278:     PetscCall(DMLabelGetStratumSize(label, values[v], &n));
    2279:     PetscCall(PetscSectionSetDof(*section, values[v], n));
    2280:   }
    2281:   PetscCall(PetscSectionSetUp(*section));
    2282:   PetscCall(PetscSectionGetStorageSize(*section, &N));
    2283:   PetscCall(PetscMalloc1(N, &points));
    2284:   for (v = 0; v < nV; ++v) {
    2285:     IS              is;
    2286:     const PetscInt *spoints;
    2287:     PetscInt        dof, off, p;

    2289:     PetscCall(PetscSectionGetDof(*section, values[v], &dof));
    2290:     PetscCall(PetscSectionGetOffset(*section, values[v], &off));
    2291:     PetscCall(DMLabelGetStratumIS(label, values[v], &is));
    2292:     PetscCall(ISGetIndices(is, &spoints));
    2293:     for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
    2294:     PetscCall(ISRestoreIndices(is, &spoints));
    2295:     PetscCall(ISDestroy(&is));
    2296:   }
    2297:   PetscCall(ISRestoreIndices(vIS, &values));
    2298:   PetscCall(ISDestroy(&vIS));
    2299:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
    2300:   PetscFunctionReturn(PETSC_SUCCESS);
    2301: }

    2303: /*@C
    2304:   DMLabelRegister - Adds a new label component implementation

    2306:   Not Collective

    2308:   Input Parameters:
    2309: + name        - The name of a new user-defined creation routine
    2310: - create_func - The creation routine itself

    2312:   Notes:
    2313:   `DMLabelRegister()` may be called multiple times to add several user-defined labels

    2315:   Example Usage:
    2316: .vb
    2317:   DMLabelRegister("my_label", MyLabelCreate);
    2318: .ve

    2320:   Then, your label type can be chosen with the procedural interface via
    2321: .vb
    2322:   DMLabelCreate(MPI_Comm, DMLabel *);
    2323:   DMLabelSetType(DMLabel, "my_label");
    2324: .ve
    2325:   or at runtime via the option
    2326: .vb
    2327:   -dm_label_type my_label
    2328: .ve

    2330:   Level: advanced

    2332: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
    2333: @*/
    2334: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
    2335: {
    2336:   PetscFunctionBegin;
    2337:   PetscCall(DMInitializePackage());
    2338:   PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
    2339:   PetscFunctionReturn(PETSC_SUCCESS);
    2340: }

    2342: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
    2343: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);

    2345: /*@C
    2346:   DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.

    2348:   Not Collective

    2350:   Level: advanced

    2352: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
    2353: @*/
    2354: PetscErrorCode DMLabelRegisterAll(void)
    2355: {
    2356:   PetscFunctionBegin;
    2357:   if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    2358:   DMLabelRegisterAllCalled = PETSC_TRUE;

    2360:   PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
    2361:   PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
    2362:   PetscFunctionReturn(PETSC_SUCCESS);
    2363: }

    2365: /*@C
    2366:   DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.

    2368:   Level: developer

    2370: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
    2371: @*/
    2372: PetscErrorCode DMLabelRegisterDestroy(void)
    2373: {
    2374:   PetscFunctionBegin;
    2375:   PetscCall(PetscFunctionListDestroy(&DMLabelList));
    2376:   DMLabelRegisterAllCalled = PETSC_FALSE;
    2377:   PetscFunctionReturn(PETSC_SUCCESS);
    2378: }

    2380: /*@
    2381:   DMLabelSetType - Sets the particular implementation for a label.

    2383:   Collective

    2385:   Input Parameters:
    2386: + label  - The label
    2387: - method - The name of the label type

    2389:   Options Database Key:
    2390: . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType`

    2392:   Level: intermediate

    2394: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
    2395: @*/
    2396: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
    2397: {
    2398:   PetscErrorCode (*r)(DMLabel);
    2399:   PetscBool match;

    2401:   PetscFunctionBegin;
    2403:   PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
    2404:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    2406:   PetscCall(DMLabelRegisterAll());
    2407:   PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
    2408:   PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);

    2410:   PetscTryTypeMethod(label, destroy);
    2411:   PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
    2412:   PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
    2413:   PetscCall((*r)(label));
    2414:   PetscFunctionReturn(PETSC_SUCCESS);
    2415: }

    2417: /*@
    2418:   DMLabelGetType - Gets the type name (as a string) from the label.

    2420:   Not Collective

    2422:   Input Parameter:
    2423: . label - The `DMLabel`

    2425:   Output Parameter:
    2426: . type - The `DMLabel` type name

    2428:   Level: intermediate

    2430: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
    2431: @*/
    2432: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
    2433: {
    2434:   PetscFunctionBegin;
    2436:   PetscAssertPointer(type, 2);
    2437:   PetscCall(DMLabelRegisterAll());
    2438:   *type = ((PetscObject)label)->type_name;
    2439:   PetscFunctionReturn(PETSC_SUCCESS);
    2440: }

    2442: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
    2443: {
    2444:   PetscFunctionBegin;
    2445:   label->ops->view         = DMLabelView_Concrete;
    2446:   label->ops->setup        = NULL;
    2447:   label->ops->duplicate    = DMLabelDuplicate_Concrete;
    2448:   label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
    2449:   PetscFunctionReturn(PETSC_SUCCESS);
    2450: }

    2452: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
    2453: {
    2454:   PetscFunctionBegin;
    2456:   PetscCall(DMLabelInitialize_Concrete(label));
    2457:   PetscFunctionReturn(PETSC_SUCCESS);
    2458: }

    2460: /*@
    2461:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
    2462:   the local section and an `PetscSF` describing the section point overlap.

    2464:   Collective

    2466:   Input Parameters:
    2467: + s                  - The `PetscSection` for the local field layout
    2468: . sf                 - The `PetscSF` describing parallel layout of the section points
    2469: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
    2470: . label              - The label specifying the points
    2471: - labelValue         - The label stratum specifying the points

    2473:   Output Parameter:
    2474: . gsection - The `PetscSection` for the global field layout

    2476:   Level: developer

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

    2481: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
    2482: @*/
    2483: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
    2484: {
    2485:   PetscInt *neg = NULL, *tmpOff = NULL;
    2486:   PetscInt  pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

    2488:   PetscFunctionBegin;
    2492:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
    2493:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    2494:   PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
    2495:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
    2496:   if (nroots >= 0) {
    2497:     PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
    2498:     PetscCall(PetscCalloc1(nroots, &neg));
    2499:     if (nroots > pEnd - pStart) {
    2500:       PetscCall(PetscCalloc1(nroots, &tmpOff));
    2501:     } else {
    2502:       tmpOff = &(*gsection)->atlasDof[-pStart];
    2503:     }
    2504:   }
    2505:   /* Mark ghost points with negative dof */
    2506:   for (p = pStart; p < pEnd; ++p) {
    2507:     PetscInt value;

    2509:     PetscCall(DMLabelGetValue(label, p, &value));
    2510:     if (value != labelValue) continue;
    2511:     PetscCall(PetscSectionGetDof(s, p, &dof));
    2512:     PetscCall(PetscSectionSetDof(*gsection, p, dof));
    2513:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
    2514:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
    2515:     if (neg) neg[p] = -(dof + 1);
    2516:   }
    2517:   PetscCall(PetscSectionSetUpBC(*gsection));
    2518:   if (nroots >= 0) {
    2519:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
    2520:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
    2521:     if (nroots > pEnd - pStart) {
    2522:       for (p = pStart; p < pEnd; ++p) {
    2523:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
    2524:       }
    2525:     }
    2526:   }
    2527:   /* Calculate new sizes, get process offset, and calculate point offsets */
    2528:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
    2529:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
    2530:     (*gsection)->atlasOff[p] = off;
    2531:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
    2532:   }
    2533:   PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
    2534:   globalOff -= off;
    2535:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
    2536:     (*gsection)->atlasOff[p] += globalOff;
    2537:     if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
    2538:   }
    2539:   /* Put in negative offsets for ghost points */
    2540:   if (nroots >= 0) {
    2541:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
    2542:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
    2543:     if (nroots > pEnd - pStart) {
    2544:       for (p = pStart; p < pEnd; ++p) {
    2545:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
    2546:       }
    2547:     }
    2548:   }
    2549:   if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
    2550:   PetscCall(PetscFree(neg));
    2551:   PetscFunctionReturn(PETSC_SUCCESS);
    2552: }

    2554: typedef struct _n_PetscSectionSym_Label {
    2555:   DMLabel              label;
    2556:   PetscCopyMode       *modes;
    2557:   PetscInt            *sizes;
    2558:   const PetscInt    ***perms;
    2559:   const PetscScalar ***rots;
    2560:   PetscInt (*minMaxOrients)[2];
    2561:   PetscInt numStrata; /* numStrata is only increasing, functions as a state */
    2562: } PetscSectionSym_Label;

    2564: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
    2565: {
    2566:   PetscInt               i, j;
    2567:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;

    2569:   PetscFunctionBegin;
    2570:   for (i = 0; i <= sl->numStrata; i++) {
    2571:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
    2572:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
    2573:         if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
    2574:         if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
    2575:       }
    2576:       if (sl->perms[i]) {
    2577:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

    2579:         PetscCall(PetscFree(perms));
    2580:       }
    2581:       if (sl->rots[i]) {
    2582:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

    2584:         PetscCall(PetscFree(rots));
    2585:       }
    2586:     }
    2587:   }
    2588:   PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
    2589:   PetscCall(DMLabelDestroy(&sl->label));
    2590:   sl->numStrata = 0;
    2591:   PetscFunctionReturn(PETSC_SUCCESS);
    2592: }

    2594: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
    2595: {
    2596:   PetscFunctionBegin;
    2597:   PetscCall(PetscSectionSymLabelReset(sym));
    2598:   PetscCall(PetscFree(sym->data));
    2599:   PetscFunctionReturn(PETSC_SUCCESS);
    2600: }

    2602: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
    2603: {
    2604:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
    2605:   PetscBool              isAscii;
    2606:   DMLabel                label = sl->label;
    2607:   const char            *name;

    2609:   PetscFunctionBegin;
    2610:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
    2611:   if (isAscii) {
    2612:     PetscInt          i, j, k;
    2613:     PetscViewerFormat format;

    2615:     PetscCall(PetscViewerGetFormat(viewer, &format));
    2616:     if (label) {
    2617:       PetscCall(PetscViewerGetFormat(viewer, &format));
    2618:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    2619:         PetscCall(PetscViewerASCIIPushTab(viewer));
    2620:         PetscCall(DMLabelView(label, viewer));
    2621:         PetscCall(PetscViewerASCIIPopTab(viewer));
    2622:       } else {
    2623:         PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
    2624:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Label '%s'\n", name));
    2625:       }
    2626:     } else {
    2627:       PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
    2628:     }
    2629:     PetscCall(PetscViewerASCIIPushTab(viewer));
    2630:     for (i = 0; i <= sl->numStrata; i++) {
    2631:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

    2633:       if (!(sl->perms[i] || sl->rots[i])) {
    2634:         PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
    2635:       } else {
    2636:         PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
    2637:         PetscCall(PetscViewerASCIIPushTab(viewer));
    2638:         PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
    2639:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    2640:           PetscCall(PetscViewerASCIIPushTab(viewer));
    2641:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
    2642:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
    2643:               PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
    2644:             } else {
    2645:               PetscInt tab;

    2647:               PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
    2648:               PetscCall(PetscViewerASCIIPushTab(viewer));
    2649:               PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
    2650:               if (sl->perms[i] && sl->perms[i][j]) {
    2651:                 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
    2652:                 PetscCall(PetscViewerASCIISetTab(viewer, 0));
    2653:                 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
    2654:                 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    2655:                 PetscCall(PetscViewerASCIISetTab(viewer, tab));
    2656:               }
    2657:               if (sl->rots[i] && sl->rots[i][j]) {
    2658:                 PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations:  "));
    2659:                 PetscCall(PetscViewerASCIISetTab(viewer, 0));
    2660: #if defined(PETSC_USE_COMPLEX)
    2661:                 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g+i*%+g", (double)PetscRealPart(sl->rots[i][j][k]), (double)PetscImaginaryPart(sl->rots[i][j][k])));
    2662: #else
    2663:                 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
    2664: #endif
    2665:                 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    2666:                 PetscCall(PetscViewerASCIISetTab(viewer, tab));
    2667:               }
    2668:               PetscCall(PetscViewerASCIIPopTab(viewer));
    2669:             }
    2670:           }
    2671:           PetscCall(PetscViewerASCIIPopTab(viewer));
    2672:         }
    2673:         PetscCall(PetscViewerASCIIPopTab(viewer));
    2674:       }
    2675:     }
    2676:     PetscCall(PetscViewerASCIIPopTab(viewer));
    2677:   }
    2678:   PetscFunctionReturn(PETSC_SUCCESS);
    2679: }

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

    2684:   Logically

    2686:   Input Parameters:
    2687: + sym   - the section symmetries
    2688: - label - the `DMLabel` describing the types of points

    2690:   Level: developer:

    2692: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
    2693: @*/
    2694: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
    2695: {
    2696:   PetscSectionSym_Label *sl;

    2698:   PetscFunctionBegin;
    2700:   sl = (PetscSectionSym_Label *)sym->data;
    2701:   if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
    2702:   if (label) {
    2703:     sl->label = label;
    2704:     PetscCall(PetscObjectReference((PetscObject)label));
    2705:     PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
    2706:     PetscCall(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));
    2707:     PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
    2708:     PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
    2709:     PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
    2710:     PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
    2711:     PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
    2712:   }
    2713:   PetscFunctionReturn(PETSC_SUCCESS);
    2714: }

    2716: /*@C
    2717:   PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum

    2719:   Logically Collective

    2721:   Input Parameters:
    2722: + sym     - the section symmetries
    2723: - stratum - the stratum value in the label that we are assigning symmetries for

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

    2732:   Level: developer

    2734: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
    2735: @*/
    2736: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
    2737: {
    2738:   PetscSectionSym_Label *sl;
    2739:   const char            *name;
    2740:   PetscInt               i;

    2742:   PetscFunctionBegin;
    2744:   sl = (PetscSectionSym_Label *)sym->data;
    2745:   PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
    2746:   for (i = 0; i <= sl->numStrata; i++) {
    2747:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

    2749:     if (stratum == value) break;
    2750:   }
    2751:   PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
    2752:   PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
    2753:   if (size) {
    2754:     PetscAssertPointer(size, 3);
    2755:     *size = sl->sizes[i];
    2756:   }
    2757:   if (minOrient) {
    2758:     PetscAssertPointer(minOrient, 4);
    2759:     *minOrient = sl->minMaxOrients[i][0];
    2760:   }
    2761:   if (maxOrient) {
    2762:     PetscAssertPointer(maxOrient, 5);
    2763:     *maxOrient = sl->minMaxOrients[i][1];
    2764:   }
    2765:   if (perms) {
    2766:     PetscAssertPointer(perms, 6);
    2767:     *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
    2768:   }
    2769:   if (rots) {
    2770:     PetscAssertPointer(rots, 7);
    2771:     *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
    2772:   }
    2773:   PetscFunctionReturn(PETSC_SUCCESS);
    2774: }

    2776: /*@C
    2777:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

    2779:   Logically

    2781:   Input Parameters:
    2782: + sym       - the section symmetries
    2783: . stratum   - the stratum value in the label that we are assigning symmetries for
    2784: . size      - the number of dofs for points in the `stratum` of the label
    2785: . minOrient - the smallest orientation for a point in this `stratum`
    2786: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
    2787: . mode      - how `sym` should copy the `perms` and `rots` arrays
    2788: . perms     - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation.  A `NULL` permutation is the identity
    2789: - 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

    2791:   Level: developer

    2793: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
    2794: @*/
    2795: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
    2796: {
    2797:   PetscSectionSym_Label *sl;
    2798:   const char            *name;
    2799:   PetscInt               i, j, k;

    2801:   PetscFunctionBegin;
    2803:   sl = (PetscSectionSym_Label *)sym->data;
    2804:   PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
    2805:   for (i = 0; i <= sl->numStrata; i++) {
    2806:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

    2808:     if (stratum == value) break;
    2809:   }
    2810:   PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
    2811:   PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
    2812:   sl->sizes[i]            = size;
    2813:   sl->modes[i]            = mode;
    2814:   sl->minMaxOrients[i][0] = minOrient;
    2815:   sl->minMaxOrients[i][1] = maxOrient;
    2816:   if (mode == PETSC_COPY_VALUES) {
    2817:     if (perms) {
    2818:       PetscInt **ownPerms;

    2820:       PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
    2821:       for (j = 0; j < maxOrient - minOrient; j++) {
    2822:         if (perms[j]) {
    2823:           PetscCall(PetscMalloc1(size, &ownPerms[j]));
    2824:           for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
    2825:         }
    2826:       }
    2827:       sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
    2828:     }
    2829:     if (rots) {
    2830:       PetscScalar **ownRots;

    2832:       PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
    2833:       for (j = 0; j < maxOrient - minOrient; j++) {
    2834:         if (rots[j]) {
    2835:           PetscCall(PetscMalloc1(size, &ownRots[j]));
    2836:           for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
    2837:         }
    2838:       }
    2839:       sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
    2840:     }
    2841:   } else {
    2842:     sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
    2843:     sl->rots[i]  = PetscSafePointerPlusOffset(rots, -minOrient);
    2844:   }
    2845:   PetscFunctionReturn(PETSC_SUCCESS);
    2846: }

    2848: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
    2849: {
    2850:   PetscInt               i, j, numStrata;
    2851:   PetscSectionSym_Label *sl;
    2852:   DMLabel                label;

    2854:   PetscFunctionBegin;
    2855:   sl        = (PetscSectionSym_Label *)sym->data;
    2856:   numStrata = sl->numStrata;
    2857:   label     = sl->label;
    2858:   for (i = 0; i < numPoints; i++) {
    2859:     PetscInt point = points[2 * i];
    2860:     PetscInt ornt  = points[2 * i + 1];

    2862:     for (j = 0; j < numStrata; j++) {
    2863:       if (label->validIS[j]) {
    2864:         PetscInt k;

    2866:         PetscCall(ISLocate(label->points[j], point, &k));
    2867:         if (k >= 0) break;
    2868:       } else {
    2869:         PetscBool has;

    2871:         PetscCall(PetscHSetIHas(label->ht[j], point, &has));
    2872:         if (has) break;
    2873:       }
    2874:     }
    2875:     PetscCheck(!(sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) || !(ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " orientation %" PetscInt_FMT " not in range [%" PetscInt_FMT ", %" PetscInt_FMT ") for stratum %" PetscInt_FMT, point, ornt, sl->minMaxOrients[j][0], sl->minMaxOrients[j][1],
    2876:                j < numStrata ? label->stratumValues[j] : label->defaultValue);
    2877:     if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
    2878:     if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
    2879:   }
    2880:   PetscFunctionReturn(PETSC_SUCCESS);
    2881: }

    2883: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
    2884: {
    2885:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
    2886:   IS                     valIS;
    2887:   const PetscInt        *values;
    2888:   PetscInt               Nv, v;

    2890:   PetscFunctionBegin;
    2891:   PetscCall(DMLabelGetNumValues(sl->label, &Nv));
    2892:   PetscCall(DMLabelGetValueIS(sl->label, &valIS));
    2893:   PetscCall(ISGetIndices(valIS, &values));
    2894:   for (v = 0; v < Nv; ++v) {
    2895:     const PetscInt      val = values[v];
    2896:     PetscInt            size, minOrient, maxOrient;
    2897:     const PetscInt    **perms;
    2898:     const PetscScalar **rots;

    2900:     PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
    2901:     PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
    2902:   }
    2903:   PetscCall(ISDestroy(&valIS));
    2904:   PetscFunctionReturn(PETSC_SUCCESS);
    2905: }

    2907: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
    2908: {
    2909:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
    2910:   DMLabel                dlabel;

    2912:   PetscFunctionBegin;
    2913:   PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
    2914:   PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
    2915:   PetscCall(DMLabelDestroy(&dlabel));
    2916:   PetscCall(PetscSectionSymCopy(sym, *dsym));
    2917:   PetscFunctionReturn(PETSC_SUCCESS);
    2918: }

    2920: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
    2921: {
    2922:   PetscSectionSym_Label *sl;

    2924:   PetscFunctionBegin;
    2925:   PetscCall(PetscNew(&sl));
    2926:   sym->ops->getpoints  = PetscSectionSymGetPoints_Label;
    2927:   sym->ops->distribute = PetscSectionSymDistribute_Label;
    2928:   sym->ops->copy       = PetscSectionSymCopy_Label;
    2929:   sym->ops->view       = PetscSectionSymView_Label;
    2930:   sym->ops->destroy    = PetscSectionSymDestroy_Label;
    2931:   sym->data            = (void *)sl;
    2932:   PetscFunctionReturn(PETSC_SUCCESS);
    2933: }

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

    2938:   Collective

    2940:   Input Parameters:
    2941: + comm  - the MPI communicator for the new symmetry
    2942: - label - the label defining the strata

    2944:   Output Parameter:
    2945: . sym - the section symmetries

    2947:   Level: developer

    2949: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
    2950: @*/
    2951: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
    2952: {
    2953:   PetscFunctionBegin;
    2954:   PetscCall(DMInitializePackage());
    2955:   PetscCall(PetscSectionSymCreate(comm, sym));
    2956:   PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
    2957:   PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
    2958:   PetscFunctionReturn(PETSC_SUCCESS);
    2959: }