Actual source code: dmlabel.c

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

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

     10:   Collective

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

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

     19:   Level: beginner

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

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

     33:   DMInitializePackage();

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

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

     52: /*
     53:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

     55:   Not collective

     57:   Input parameter:
     58: + label - The DMLabel
     59: - v - The stratum value

     61:   Output parameter:
     62: . label - The DMLabel with stratum in sorted list format

     64:   Level: developer

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

     74:   if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
     76:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
     77:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
     78:   PetscMalloc1(label->stratumSizes[v], &pointArray);
     79:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
     80:   PetscHSetIClear(label->ht[v]);
     81:   PetscSortInt(label->stratumSizes[v], pointArray);
     82:   if (label->bt) {
     83:     for (p = 0; p < label->stratumSizes[v]; ++p) {
     84:       const PetscInt point = pointArray[p];
     85:       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
     86:       PetscBTSet(label->bt, point - label->pStart);
     87:     }
     88:   }
     89:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
     90:     ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
     91:     PetscFree(pointArray);
     92:   } else {
     93:     ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
     94:   }
     95:   PetscObjectSetName((PetscObject) is, "indices");
     96:   label->points[v]  = is;
     97:   label->validIS[v] = PETSC_TRUE;
     98:   PetscObjectStateIncrease((PetscObject) label);
     99:   return(0);
    100: }

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

    105:   Not collective

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

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

    113:   Level: developer

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

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

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

    132:   Not collective

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

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

    141:   Level: developer

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

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

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

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

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

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

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

    257: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
    258: {
    261:   DMLabelLookupStratum(label, value, index);
    262:   if (*index < 0) {DMLabelNewStratum(label, value, index);}
    263:   return(0);
    264: }

    266: /*@
    267:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

    269:   Input Parameter:
    270: + label - The DMLabel
    271: - value - The stratum value

    273:   Level: beginner

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

    284:   DMLabelLookupAddStratum(label, value, &v);
    285:   return(0);
    286: }

    288: /*@
    289:   DMLabelAddStrata - Adds new stratum values in a DMLabel

    291:   Not collective

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

    298:   Level: beginner

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

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

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

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

    355:   Not collective

    357:   Input Parameter:
    358: + label - The DMLabel
    359: - valueIS - Index set with stratum values

    361:   Level: beginner

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

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

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

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

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

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

    412: /*@C
    413:   DMLabelView - View the label

    415:   Collective on viewer

    417:   Input Parameters:
    418: + label - The DMLabel
    419: - viewer - The PetscViewer

    421:   Level: intermediate

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

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

    442: /*@
    443:   DMLabelReset - Destroys internal data structures in a DMLabel

    445:   Not collective

    447:   Input Parameter:
    448: . label - The DMLabel

    450:   Level: beginner

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

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

    483: /*@
    484:   DMLabelDestroy - Destroys a DMLabel

    486:   Collective on label

    488:   Input Parameter:
    489: . label - The DMLabel

    491:   Level: beginner

    493: .seealso: DMLabelReset(), DMLabelCreate()
    494: @*/
    495: PetscErrorCode DMLabelDestroy(DMLabel *label)
    496: {

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

    509: /*@
    510:   DMLabelDuplicate - Duplicates a DMLabel

    512:   Collective on label

    514:   Input Parameter:
    515: . label - The DMLabel

    517:   Output Parameter:
    518: . labelnew - location to put new vector

    520:   Level: intermediate

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

    532:   DMLabelMakeAllValid_Private(label);
    533:   PetscObjectGetName((PetscObject) label, &name);
    534:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

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

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

    562:   Not collective

    564:   Input Parameter:
    565: . label  - The DMLabel

    567:   Level: intermediate

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

    578:   DMLabelMakeAllValid_Private(label);
    579:   for (v = 0; v < label->numStrata; ++v) {
    580:     const PetscInt *points;
    581:     PetscInt       i;

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

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

    598: /*@
    599:   DMLabelCreateIndex - Create an index structure for membership determination

    601:   Not collective

    603:   Input Parameters:
    604: + label  - The DMLabel
    605: . pStart - The smallest point
    606: - pEnd   - The largest point + 1

    608:   Level: intermediate

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

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

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

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

    641: /*@
    642:   DMLabelDestroyIndex - Destroy the index structure

    644:   Not collective

    646:   Input Parameter:
    647: . label - the DMLabel

    649:   Level: intermediate

    651: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
    652: @*/
    653: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
    654: {

    659:   label->pStart = -1;
    660:   label->pEnd   = -1;
    661:   PetscBTDestroy(&label->bt);
    662:   return(0);
    663: }

    665: /*@
    666:   DMLabelGetBounds - Return the smallest and largest point in the label

    668:   Not collective

    670:   Input Parameter:
    671: . label - the DMLabel

    673:   Output Parameters:
    674: + pStart - The smallest point
    675: - pEnd   - The largest point + 1

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

    679:   Level: intermediate

    681: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
    682: @*/
    683: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
    684: {

    689:   if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
    690:   if (pStart) {
    692:     *pStart = label->pStart;
    693:   }
    694:   if (pEnd) {
    696:     *pEnd = label->pEnd;
    697:   }
    698:   return(0);
    699: }

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

    704:   Not collective

    706:   Input Parameters:
    707: + label - the DMLabel
    708: - value - the value

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

    713:   Level: developer

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

    725:   DMLabelLookupStratum(label, value, &v);
    726:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
    727:   return(0);
    728: }

    730: /*@
    731:   DMLabelHasPoint - Determine whether a label assigns a value to a point

    733:   Not collective

    735:   Input Parameters:
    736: + label - the DMLabel
    737: - point - the point

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

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

    744:   Level: developer

    746: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
    747: @*/
    748: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
    749: {

    755:   DMLabelMakeAllValid_Private(label);
    756:   if (PetscDefined(USE_DEBUG)) {
    757:     if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
    758:     if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
    759:   }
    760:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
    761:   return(0);
    762: }

    764: /*@
    765:   DMLabelStratumHasPoint - Return true if the stratum contains a point

    767:   Not collective

    769:   Input Parameters:
    770: + label - the DMLabel
    771: . value - the stratum value
    772: - point - the point

    774:   Output Parameter:
    775: . contains - true if the stratum contains the point

    777:   Level: intermediate

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

    789:   *contains = PETSC_FALSE;
    790:   DMLabelLookupStratum(label, value, &v);
    791:   if (v < 0) return(0);

    793:   if (label->validIS[v]) {
    794:     PetscInt i;

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

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

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

    811:   Not collective

    813:   Input parameter:
    814: . label - a DMLabel object

    816:   Output parameter:
    817: . defaultValue - the default value

    819:   Level: beginner

    821: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
    822: @*/
    823: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
    824: {
    827:   *defaultValue = label->defaultValue;
    828:   return(0);
    829: }

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

    835:   Not collective

    837:   Input parameter:
    838: . label - a DMLabel object

    840:   Output parameter:
    841: . defaultValue - the default value

    843:   Level: beginner

    845: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
    846: @*/
    847: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
    848: {
    851:   label->defaultValue = defaultValue;
    852:   return(0);
    853: }

    855: /*@
    856:   DMLabelGetValue - Return the value a label assigns to a point, or the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue())

    858:   Not collective

    860:   Input Parameters:
    861: + label - the DMLabel
    862: - point - the point

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

    867:   Level: intermediate

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

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

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

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

    902: /*@
    903:   DMLabelSetValue - Set the value a label assigns to a point.  If the value is the same as the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue() to something different), then this function will do nothing.

    905:   Not collective

    907:   Input Parameters:
    908: + label - the DMLabel
    909: . point - the point
    910: - value - The point value

    912:   Level: intermediate

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

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

    932: /*@
    933:   DMLabelClearValue - Clear the value a label assigns to a point

    935:   Not collective

    937:   Input Parameters:
    938: + label - the DMLabel
    939: . point - the point
    940: - value - The point value

    942:   Level: intermediate

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

    953:   /* Find label value */
    954:   DMLabelLookupStratum(label, value, &v);
    955:   if (v < 0) return(0);

    957:   if (label->bt) {
    958:     if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
    959:     PetscBTClear(label->bt, point - label->pStart);
    960:   }

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

    968: /*@
    969:   DMLabelInsertIS - Set all points in the IS to a value

    971:   Not collective

    973:   Input Parameters:
    974: + label - the DMLabel
    975: . is    - the point IS
    976: - value - The point value

    978:   Level: intermediate

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

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

    1003: /*@
    1004:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

    1006:   Not collective

    1008:   Input Parameter:
    1009: . label - the DMLabel

    1011:   Output Paramater:
    1012: . numValues - the number of values

    1014:   Level: intermediate

    1016: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
    1017: @*/
    1018: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
    1019: {
    1023:   *numValues = label->numStrata;
    1024:   return(0);
    1025: }

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

    1030:   Not collective

    1032:   Input Parameter:
    1033: . label - the DMLabel

    1035:   Output Paramater:
    1036: . is    - the value IS

    1038:   Level: intermediate

    1040: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
    1041: @*/
    1042: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
    1043: {

    1049:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
    1050:   return(0);
    1051: }

    1053: /*@
    1054:   DMLabelHasStratum - Determine whether points exist with the given value

    1056:   Not collective

    1058:   Input Parameters:
    1059: + label - the DMLabel
    1060: - value - the stratum value

    1062:   Output Paramater:
    1063: . exists - Flag saying whether points exist

    1065:   Level: intermediate

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

    1077:   DMLabelLookupStratum(label, value, &v);
    1078:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
    1079:   return(0);
    1080: }

    1082: /*@
    1083:   DMLabelGetStratumSize - Get the size of a stratum

    1085:   Not collective

    1087:   Input Parameters:
    1088: + label - the DMLabel
    1089: - value - the stratum value

    1091:   Output Paramater:
    1092: . size - The number of points in the stratum

    1094:   Level: intermediate

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

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

    1114: /*@
    1115:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

    1117:   Not collective

    1119:   Input Parameters:
    1120: + label - the DMLabel
    1121: - value - the stratum value

    1123:   Output Paramaters:
    1124: + start - the smallest point in the stratum
    1125: - end - the largest point in the stratum

    1127:   Level: intermediate

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

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

    1150: /*@
    1151:   DMLabelGetStratumIS - Get an IS with the stratum points

    1153:   Not collective

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

    1159:   Output Paramater:
    1160: . points - The stratum points

    1162:   Level: intermediate

    1164:   Notes:
    1165:   The output IS should be destroyed when no longer needed.
    1166:   Returns NULL if the stratum is empty.

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

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

    1187: /*@
    1188:   DMLabelSetStratumIS - Set the stratum points using an IS

    1190:   Not collective

    1192:   Input Parameters:
    1193: + label - the DMLabel
    1194: . value - the stratum value
    1195: - points - The stratum points

    1197:   Level: intermediate

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

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

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

    1226:       if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
    1227:       PetscBTSet(label->bt, point - label->pStart);
    1228:     }
    1229:   }
    1230:   return(0);
    1231: }

    1233: /*@
    1234:   DMLabelClearStratum - Remove a stratum

    1236:   Not collective

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

    1242:   Level: intermediate

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

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

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

    1264:         if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
    1265:         PetscBTClear(label->bt, point - label->pStart);
    1266:       }
    1267:       ISRestoreIndices(label->points[v], &points);
    1268:     }
    1269:     label->stratumSizes[v] = 0;
    1270:     ISDestroy(&label->points[v]);
    1271:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
    1272:     PetscObjectSetName((PetscObject) label->points[v], "indices");
    1273:     PetscObjectStateIncrease((PetscObject) label);
    1274:   } else {
    1275:     PetscHSetIClear(label->ht[v]);
    1276:   }
    1277:   return(0);
    1278: }

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

    1283:   Not collective

    1285:   Input Parameters:
    1286: + label  - The DMLabel
    1287: . value  - The label value for all points
    1288: . pStart - The first point
    1289: - pEnd   - A point beyond all marked points

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

    1293:   Level: intermediate

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

    1303:   ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
    1304:   DMLabelSetStratumIS(label, value, pIS);
    1305:   ISDestroy(&pIS);
    1306:   return(0);
    1307: }

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

    1312:   Not collective

    1314:   Input Parameters:
    1315: + label - the DMLabel
    1316: . start - the first point kept
    1317: - end - one more than the last point kept

    1319:   Level: intermediate

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

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

    1339: /*@
    1340:   DMLabelPermute - Create a new label with permuted points

    1342:   Not collective

    1344:   Input Parameters:
    1345: + label - the DMLabel
    1346: - permutation - the point permutation

    1348:   Output Parameter:
    1349: . labelnew - the new label containing the permuted points

    1351:   Level: intermediate

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

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

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

    1379:       if ((point < 0) || (point >= numPoints)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [0, %D) for the remapping", point, numPoints);
    1380:       pointsNew[q] = perm[point];
    1381:     }
    1382:     ISRestoreIndices((*labelNew)->points[v],&points);
    1383:     PetscSortInt(size, pointsNew);
    1384:     ISDestroy(&((*labelNew)->points[v]));
    1385:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
    1386:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
    1387:       PetscFree(pointsNew);
    1388:     } else {
    1389:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
    1390:     }
    1391:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
    1392:   }
    1393:   ISRestoreIndices(permutation, &perm);
    1394:   if (label->bt) {
    1395:     PetscBTDestroy(&label->bt);
    1396:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
    1397:   }
    1398:   return(0);
    1399: }

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

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

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

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

    1466: /*@
    1467:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

    1469:   Collective on sf

    1471:   Input Parameters:
    1472: + label - the DMLabel
    1473: - sf    - the map from old to new distribution

    1475:   Output Parameter:
    1476: . labelnew - the new redistributed label

    1478:   Level: intermediate

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

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

    1582: /*@
    1583:   DMLabelGather - Gather all label values from leafs into roots

    1585:   Collective on sf

    1587:   Input Parameters:
    1588: + label - the DMLabel
    1589: - sf - the Star Forest point communication map

    1591:   Output Parameters:
    1592: . labelNew - the new DMLabel with localised leaf values

    1594:   Level: developer

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

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

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

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

    1672: /*@
    1673:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

    1675:   Not collective

    1677:   Input Parameter:
    1678: . label - the DMLabel

    1680:   Output Parameters:
    1681: + section - the section giving offsets for each stratum
    1682: - is - An IS containing all the label points

    1684:   Level: developer

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

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

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

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

    1736: /*@
    1737:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
    1738:   the local section and an SF describing the section point overlap.

    1740:   Collective on sf

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

    1749:   Output Parameter:
    1750:   . gsection - The PetscSection for the global field layout

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

    1754:   Level: developer

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

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

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

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

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

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

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

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

    1868: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
    1869: {

    1873:   PetscSectionSymLabelReset(sym);
    1874:   PetscFree(sym->data);
    1875:   return(0);
    1876: }

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

    1887:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
    1888:   if (isAscii) {
    1889:     PetscInt          i, j, k;
    1890:     PetscViewerFormat format;

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

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

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

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

    1961:   Logically collective on sym

    1963:   Input parameters:
    1964: + sym - the section symmetries
    1965: - label - the DMLabel describing the types of points

    1967:   Level: developer:

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

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

    1994: /*@C
    1995:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

    1997:   Logically collective on sym

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

    2009:   Level: developer

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

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

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

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

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

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

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

    2082:     for (j = 0; j < numStrata; j++) {
    2083:       if (label->validIS[j]) {
    2084:         PetscInt k;

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

    2091:         PetscHSetIHas(label->ht[j], point, &has);
    2092:         if (has) break;
    2093:       }
    2094:     }
    2095:     if ((sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) && (ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1])) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D orientation %D not in range [%D, %D) for stratum %D",point,ornt,sl->minMaxOrients[j][0],sl->minMaxOrients[j][1],j < numStrata ? label->stratumValues[j] : label->defaultValue);
    2096:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
    2097:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
    2098:   }
    2099:   return(0);
    2100: }

    2102: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
    2103: {
    2104:   PetscSectionSym_Label *sl;
    2105:   PetscErrorCode        ierr;

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

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

    2119:   Collective

    2121:   Input Parameters:
    2122: + comm - the MPI communicator for the new symmetry
    2123: - label - the label defining the strata

    2125:   Output Parameters:
    2126: . sym - the section symmetries

    2128:   Level: developer

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

    2137:   DMInitializePackage();
    2138:   PetscSectionSymCreate(comm,sym);
    2139:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
    2140:   PetscSectionSymLabelSetLabel(*sym,label);
    2141:   return(0);
    2142: }