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: /*@C
 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));

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

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

 57:   Collective

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

 62:   Level: intermediate

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

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

 77:   Not collective

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

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

 86:   Level: developer

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

 95:   if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) return PETSC_SUCCESS;
 96:   PetscFunctionBegin;
 97:   PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
 98:   PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
 99:   PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
100:   PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
101:   PetscCall(PetscHSetIClear(label->ht[v]));
102:   PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
103:   if (label->bt) {
104:     for (p = 0; p < label->stratumSizes[v]; ++p) {
105:       const PetscInt point = pointArray[p];
106:       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);
107:       PetscCall(PetscBTSet(label->bt, point - label->pStart));
108:     }
109:   }
110:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
111:     PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
112:     PetscCall(PetscFree(pointArray));
113:   } else {
114:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
115:   }
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:   if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) return PETSC_SUCCESS;
169:   PetscFunctionBegin;
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 iascii;

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

457: /*@C
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: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
632: @*/
633: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
634: {
635:   const char *name0, *name1;
636:   char        msg[PETSC_MAX_PATH_LEN] = "";
637:   PetscBool   eq;
638:   PetscMPIInt rank;

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

651:     PetscCall(DMLabelGetDefaultValue(l0, &v0));
652:     PetscCall(DMLabelGetDefaultValue(l1, &v1));
653:     eq = (PetscBool)(v0 == v1);
654:     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));
655:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
656:     if (!eq) goto finish;
657:   }
658:   {
659:     IS is0, is1;

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

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

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

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

711:   Not Collective

713:   Input Parameter:
714: . label - The `DMLabel`

716:   Level: intermediate

718: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
719: @*/
720: PetscErrorCode DMLabelComputeIndex(DMLabel label)
721: {
722:   PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;

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

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

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

746: /*@
747:   DMLabelCreateIndex - Create an index structure for membership determination

749:   Not Collective

751:   Input Parameters:
752: + label  - The `DMLabel`
753: . pStart - The smallest point
754: - pEnd   - The largest point + 1

756:   Level: intermediate

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

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

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

782:       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);
783:       PetscCall(PetscBTSet(label->bt, point - pStart));
784:     }
785:     PetscCall(ISRestoreIndices(label->points[v], &points));
786:     PetscCall(ISDestroy(&pointIS));
787:   }
788:   PetscFunctionReturn(PETSC_SUCCESS);
789: }

791: /*@
792:   DMLabelDestroyIndex - Destroy the index structure

794:   Not Collective

796:   Input Parameter:
797: . label - the `DMLabel`

799:   Level: intermediate

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

813: /*@
814:   DMLabelGetBounds - Return the smallest and largest point in the label

816:   Not Collective

818:   Input Parameter:
819: . label - the `DMLabel`

821:   Output Parameters:
822: + pStart - The smallest point
823: - pEnd   - The largest point + 1

825:   Level: intermediate

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

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

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

851:   Not Collective

853:   Input Parameters:
854: + label - the `DMLabel`
855: - value - the value

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

860:   Level: developer

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

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

876: /*@
877:   DMLabelHasPoint - Determine whether a label assigns a value to a point

879:   Not Collective

881:   Input Parameters:
882: + label - the `DMLabel`
883: - point - the point

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

888:   Level: developer

890:   Note:
891:   The user must call `DMLabelCreateIndex()` before this function.

893: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
894: @*/
895: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
896: {
897:   PetscFunctionBeginHot;
899:   PetscAssertPointer(contains, 3);
900:   PetscCall(DMLabelMakeAllValid_Private(label));
901:   if (PetscDefined(USE_DEBUG)) {
902:     PetscCheck(label->bt, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
903:     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);
904:   }
905:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
906:   PetscFunctionReturn(PETSC_SUCCESS);
907: }

909: /*@
910:   DMLabelStratumHasPoint - Return true if the stratum contains a point

912:   Not Collective

914:   Input Parameters:
915: + label - the `DMLabel`
916: . value - the stratum value
917: - point - the point

919:   Output Parameter:
920: . contains - true if the stratum contains the point

922:   Level: intermediate

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

934:     PetscCall(DMLabelGetValue(label, point, &pointVal));
935:     *contains = (PetscBool)(pointVal == value);
936:   } else {
937:     PetscInt v;

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

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

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

963:   Not Collective

965:   Input Parameter:
966: . label - a `DMLabel` object

968:   Output Parameter:
969: . defaultValue - the default value

971:   Level: beginner

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

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

987:   Not Collective

989:   Input Parameter:
990: . label - a `DMLabel` object

992:   Output Parameter:
993: . defaultValue - the default value

995:   Level: beginner

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

1007: /*@
1008:   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
1009:   `DMLabelSetDefaultValue()`)

1011:   Not Collective

1013:   Input Parameters:
1014: + label - the `DMLabel`
1015: - point - the point

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

1020:   Level: intermediate

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

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

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

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

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

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

1065:   Not Collective

1067:   Input Parameters:
1068: + label - the `DMLabel`
1069: . point - the point
1070: - value - The point value

1072:   Level: intermediate

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

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

1092: /*@
1093:   DMLabelClearValue - Clear the value a label assigns to a point

1095:   Not Collective

1097:   Input Parameters:
1098: + label - the `DMLabel`
1099: . point - the point
1100: - value - The point value

1102:   Level: intermediate

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

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

1117:   if (label->bt) {
1118:     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);
1119:     PetscCall(PetscBTClear(label->bt, point - label->pStart));
1120:   }

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

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

1131:   Not Collective

1133:   Input Parameters:
1134: + label - the `DMLabel`
1135: . is    - the point `IS`
1136: - value - The point value

1138:   Level: intermediate

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

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

1163: /*@
1164:   DMLabelGetNumValues - Get the number of values that the `DMLabel` takes

1166:   Not Collective

1168:   Input Parameter:
1169: . label - the `DMLabel`

1171:   Output Parameter:
1172: . numValues - the number of values

1174:   Level: intermediate

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

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

1190:   Not Collective

1192:   Input Parameter:
1193: . label - the `DMLabel`

1195:   Output Parameter:
1196: . values - the value `IS`

1198:   Level: intermediate

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

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

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

1219:   Not Collective

1221:   Input Parameter:
1222: . label - the `DMLabel`

1224:   Output Parameter:
1225: . values - the value `IS`

1227:   Level: intermediate

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

1233: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1234: @*/
1235: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1236: {
1237:   PetscInt  i, j;
1238:   PetscInt *valuesArr;

1240:   PetscFunctionBegin;
1242:   PetscAssertPointer(values, 2);
1243:   PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1244:   for (i = 0, j = 0; i < label->numStrata; i++) {
1245:     PetscInt n;

1247:     PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1248:     if (n) valuesArr[j++] = label->stratumValues[i];
1249:   }
1250:   if (j == label->numStrata) {
1251:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1252:   } else {
1253:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1254:   }
1255:   PetscCall(PetscFree(valuesArr));
1256:   PetscFunctionReturn(PETSC_SUCCESS);
1257: }

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

1262:   Not Collective

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

1268:   Output Parameter:
1269: . index - the index of value in the list of values

1271:   Level: intermediate

1273: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1274: @*/
1275: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1276: {
1277:   PetscInt v;

1279:   PetscFunctionBegin;
1281:   PetscAssertPointer(index, 3);
1282:   /* Do not assume they are sorted */
1283:   for (v = 0; v < label->numStrata; ++v)
1284:     if (label->stratumValues[v] == value) break;
1285:   if (v >= label->numStrata) *index = -1;
1286:   else *index = v;
1287:   PetscFunctionReturn(PETSC_SUCCESS);
1288: }

1290: /*@
1291:   DMLabelHasStratum - Determine whether points exist with the given value

1293:   Not Collective

1295:   Input Parameters:
1296: + label - the `DMLabel`
1297: - value - the stratum value

1299:   Output Parameter:
1300: . exists - Flag saying whether points exist

1302:   Level: intermediate

1304: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1305: @*/
1306: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1307: {
1308:   PetscInt v;

1310:   PetscFunctionBegin;
1312:   PetscAssertPointer(exists, 3);
1313:   PetscCall(DMLabelLookupStratum(label, value, &v));
1314:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1315:   PetscFunctionReturn(PETSC_SUCCESS);
1316: }

1318: /*@
1319:   DMLabelGetStratumSize - Get the size of a stratum

1321:   Not Collective

1323:   Input Parameters:
1324: + label - the `DMLabel`
1325: - value - the stratum value

1327:   Output Parameter:
1328: . size - The number of points in the stratum

1330:   Level: intermediate

1332: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1333: @*/
1334: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1335: {
1336:   PetscInt v;

1338:   PetscFunctionBegin;
1340:   PetscAssertPointer(size, 3);
1341:   PetscCall(DMLabelLookupStratum(label, value, &v));
1342:   PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1343:   PetscFunctionReturn(PETSC_SUCCESS);
1344: }

1346: /*@
1347:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1349:   Not Collective

1351:   Input Parameters:
1352: + label - the `DMLabel`
1353: - value - the stratum value

1355:   Output Parameters:
1356: + start - the smallest point in the stratum
1357: - end   - the largest point in the stratum

1359:   Level: intermediate

1361: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1362: @*/
1363: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1364: {
1365:   IS       is;
1366:   PetscInt v, min, max;

1368:   PetscFunctionBegin;
1370:   if (start) {
1371:     PetscAssertPointer(start, 3);
1372:     *start = -1;
1373:   }
1374:   if (end) {
1375:     PetscAssertPointer(end, 4);
1376:     *end = -1;
1377:   }
1378:   PetscCall(DMLabelLookupStratum(label, value, &v));
1379:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1380:   PetscCall(DMLabelMakeValid_Private(label, v));
1381:   if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1382:   PetscUseTypeMethod(label, getstratumis, v, &is);
1383:   PetscCall(ISGetMinMax(is, &min, &max));
1384:   PetscCall(ISDestroy(&is));
1385:   if (start) *start = min;
1386:   if (end) *end = max + 1;
1387:   PetscFunctionReturn(PETSC_SUCCESS);
1388: }

1390: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1391: {
1392:   PetscFunctionBegin;
1393:   PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1394:   *pointIS = label->points[v];
1395:   PetscFunctionReturn(PETSC_SUCCESS);
1396: }

1398: /*@
1399:   DMLabelGetStratumIS - Get an `IS` with the stratum points

1401:   Not Collective

1403:   Input Parameters:
1404: + label - the `DMLabel`
1405: - value - the stratum value

1407:   Output Parameter:
1408: . points - The stratum points

1410:   Level: intermediate

1412:   Notes:
1413:   The output `IS` should be destroyed when no longer needed.
1414:   Returns `NULL` if the stratum is empty.

1416: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1417: @*/
1418: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1419: {
1420:   PetscInt v;

1422:   PetscFunctionBegin;
1424:   PetscAssertPointer(points, 3);
1425:   *points = NULL;
1426:   PetscCall(DMLabelLookupStratum(label, value, &v));
1427:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1428:   PetscCall(DMLabelMakeValid_Private(label, v));
1429:   PetscUseTypeMethod(label, getstratumis, v, points);
1430:   PetscFunctionReturn(PETSC_SUCCESS);
1431: }

1433: /*@
1434:   DMLabelSetStratumIS - Set the stratum points using an `IS`

1436:   Not Collective

1438:   Input Parameters:
1439: + label - the `DMLabel`
1440: . value - the stratum value
1441: - is    - The stratum points

1443:   Level: intermediate

1445: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1446: @*/
1447: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1448: {
1449:   PetscInt v;

1451:   PetscFunctionBegin;
1454:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1455:   PetscCall(DMLabelLookupAddStratum(label, value, &v));
1456:   if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1457:   PetscCall(DMLabelClearStratum(label, value));
1458:   PetscCall(ISGetLocalSize(is, &(label->stratumSizes[v])));
1459:   PetscCall(PetscObjectReference((PetscObject)is));
1460:   PetscCall(ISDestroy(&(label->points[v])));
1461:   label->points[v]  = is;
1462:   label->validIS[v] = PETSC_TRUE;
1463:   PetscCall(PetscObjectStateIncrease((PetscObject)label));
1464:   if (label->bt) {
1465:     const PetscInt *points;
1466:     PetscInt        p;

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

1472:       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);
1473:       PetscCall(PetscBTSet(label->bt, point - label->pStart));
1474:     }
1475:   }
1476:   PetscFunctionReturn(PETSC_SUCCESS);
1477: }

1479: /*@
1480:   DMLabelClearStratum - Remove a stratum

1482:   Not Collective

1484:   Input Parameters:
1485: + label - the `DMLabel`
1486: - value - the stratum value

1488:   Level: intermediate

1490: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1491: @*/
1492: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1493: {
1494:   PetscInt v;

1496:   PetscFunctionBegin;
1498:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1499:   PetscCall(DMLabelLookupStratum(label, value, &v));
1500:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1501:   if (label->validIS[v]) {
1502:     if (label->bt) {
1503:       PetscInt        i;
1504:       const PetscInt *points;

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

1510:         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);
1511:         PetscCall(PetscBTClear(label->bt, point - label->pStart));
1512:       }
1513:       PetscCall(ISRestoreIndices(label->points[v], &points));
1514:     }
1515:     label->stratumSizes[v] = 0;
1516:     PetscCall(ISDestroy(&label->points[v]));
1517:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1518:     PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1519:     PetscCall(PetscObjectStateIncrease((PetscObject)label));
1520:   } else {
1521:     PetscCall(PetscHSetIClear(label->ht[v]));
1522:   }
1523:   PetscFunctionReturn(PETSC_SUCCESS);
1524: }

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

1529:   Not Collective

1531:   Input Parameters:
1532: + label  - The `DMLabel`
1533: . value  - The label value for all points
1534: . pStart - The first point
1535: - pEnd   - A point beyond all marked points

1537:   Level: intermediate

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

1542: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1543: @*/
1544: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1545: {
1546:   IS pIS;

1548:   PetscFunctionBegin;
1549:   PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1550:   PetscCall(DMLabelSetStratumIS(label, value, pIS));
1551:   PetscCall(ISDestroy(&pIS));
1552:   PetscFunctionReturn(PETSC_SUCCESS);
1553: }

1555: /*@
1556:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

1558:   Not Collective

1560:   Input Parameters:
1561: + label - The `DMLabel`
1562: . value - The label value
1563: - p     - A point with this value

1565:   Output Parameter:
1566: . 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

1568:   Level: intermediate

1570: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1571: @*/
1572: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1573: {
1574:   IS              pointIS;
1575:   const PetscInt *indices;
1576:   PetscInt        v;

1578:   PetscFunctionBegin;
1580:   PetscAssertPointer(index, 4);
1581:   *index = -1;
1582:   PetscCall(DMLabelLookupStratum(label, value, &v));
1583:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1584:   PetscCall(DMLabelMakeValid_Private(label, v));
1585:   PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1586:   PetscCall(ISGetIndices(pointIS, &indices));
1587:   PetscCall(PetscFindInt(p, label->stratumSizes[v], indices, index));
1588:   PetscCall(ISRestoreIndices(pointIS, &indices));
1589:   PetscCall(ISDestroy(&pointIS));
1590:   PetscFunctionReturn(PETSC_SUCCESS);
1591: }

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

1596:   Not Collective

1598:   Input Parameters:
1599: + label - the `DMLabel`
1600: . start - the first point kept
1601: - end   - one more than the last point kept

1603:   Level: intermediate

1605: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1606: @*/
1607: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1608: {
1609:   PetscInt v;

1611:   PetscFunctionBegin;
1613:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1614:   PetscCall(DMLabelDestroyIndex(label));
1615:   PetscCall(DMLabelMakeAllValid_Private(label));
1616:   for (v = 0; v < label->numStrata; ++v) {
1617:     PetscCall(ISGeneralFilter(label->points[v], start, end));
1618:     PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1619:   }
1620:   PetscCall(DMLabelCreateIndex(label, start, end));
1621:   PetscFunctionReturn(PETSC_SUCCESS);
1622: }

1624: /*@
1625:   DMLabelPermute - Create a new label with permuted points

1627:   Not Collective

1629:   Input Parameters:
1630: + label       - the `DMLabel`
1631: - permutation - the point permutation

1633:   Output Parameter:
1634: . labelNew - the new label containing the permuted points

1636:   Level: intermediate

1638: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1639: @*/
1640: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1641: {
1642:   const PetscInt *perm;
1643:   PetscInt        numValues, numPoints, v, q;

1645:   PetscFunctionBegin;
1648:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1649:   PetscCall(DMLabelMakeAllValid_Private(label));
1650:   PetscCall(DMLabelDuplicate(label, labelNew));
1651:   PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1652:   PetscCall(ISGetLocalSize(permutation, &numPoints));
1653:   PetscCall(ISGetIndices(permutation, &perm));
1654:   for (v = 0; v < numValues; ++v) {
1655:     const PetscInt  size = (*labelNew)->stratumSizes[v];
1656:     const PetscInt *points;
1657:     PetscInt       *pointsNew;

1659:     PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1660:     PetscCall(PetscCalloc1(size, &pointsNew));
1661:     for (q = 0; q < size; ++q) {
1662:       const PetscInt point = points[q];

1664:       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);
1665:       pointsNew[q] = perm[point];
1666:     }
1667:     PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1668:     PetscCall(PetscSortInt(size, pointsNew));
1669:     PetscCall(ISDestroy(&((*labelNew)->points[v])));
1670:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1671:       PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1672:       PetscCall(PetscFree(pointsNew));
1673:     } else {
1674:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1675:     }
1676:     PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1677:   }
1678:   PetscCall(ISRestoreIndices(permutation, &perm));
1679:   if (label->bt) {
1680:     PetscCall(PetscBTDestroy(&label->bt));
1681:     PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1682:   }
1683:   PetscFunctionReturn(PETSC_SUCCESS);
1684: }

1686: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1687: {
1688:   MPI_Comm     comm;
1689:   PetscInt     s, l, nroots, nleaves, offset, size;
1690:   PetscInt    *remoteOffsets, *rootStrata, *rootIdx;
1691:   PetscSection rootSection;
1692:   PetscSF      labelSF;

1694:   PetscFunctionBegin;
1695:   if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1696:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1697:   /* Build a section of stratum values per point, generate the according SF
1698:      and distribute point-wise stratum values to leaves. */
1699:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1700:   PetscCall(PetscSectionCreate(comm, &rootSection));
1701:   PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1702:   if (label) {
1703:     for (s = 0; s < label->numStrata; ++s) {
1704:       const PetscInt *points;

1706:       PetscCall(ISGetIndices(label->points[s], &points));
1707:       for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1708:       PetscCall(ISRestoreIndices(label->points[s], &points));
1709:     }
1710:   }
1711:   PetscCall(PetscSectionSetUp(rootSection));
1712:   /* Create a point-wise array of stratum values */
1713:   PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1714:   PetscCall(PetscMalloc1(size, &rootStrata));
1715:   PetscCall(PetscCalloc1(nroots, &rootIdx));
1716:   if (label) {
1717:     for (s = 0; s < label->numStrata; ++s) {
1718:       const PetscInt *points;

1720:       PetscCall(ISGetIndices(label->points[s], &points));
1721:       for (l = 0; l < label->stratumSizes[s]; l++) {
1722:         const PetscInt p = points[l];
1723:         PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1724:         rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1725:       }
1726:       PetscCall(ISRestoreIndices(label->points[s], &points));
1727:     }
1728:   }
1729:   /* Build SF that maps label points to remote processes */
1730:   PetscCall(PetscSectionCreate(comm, leafSection));
1731:   PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1732:   PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1733:   PetscCall(PetscFree(remoteOffsets));
1734:   /* Send the strata for each point over the derived SF */
1735:   PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1736:   PetscCall(PetscMalloc1(size, leafStrata));
1737:   PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1738:   PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1739:   /* Clean up */
1740:   PetscCall(PetscFree(rootStrata));
1741:   PetscCall(PetscFree(rootIdx));
1742:   PetscCall(PetscSectionDestroy(&rootSection));
1743:   PetscCall(PetscSFDestroy(&labelSF));
1744:   PetscFunctionReturn(PETSC_SUCCESS);
1745: }

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

1750:   Collective

1752:   Input Parameters:
1753: + label - the `DMLabel`
1754: - sf    - the map from old to new distribution

1756:   Output Parameter:
1757: . labelNew - the new redistributed label

1759:   Level: intermediate

1761: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1762: @*/
1763: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1764: {
1765:   MPI_Comm     comm;
1766:   PetscSection leafSection;
1767:   PetscInt     p, pStart, pEnd, s, size, dof, offset, stratum;
1768:   PetscInt    *leafStrata, *strataIdx;
1769:   PetscInt   **points;
1770:   const char  *lname = NULL;
1771:   char        *name;
1772:   PetscInt     nameSize;
1773:   PetscHSetI   stratumHash;
1774:   size_t       len = 0;
1775:   PetscMPIInt  rank;

1777:   PetscFunctionBegin;
1779:   if (label) {
1781:     PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1782:     PetscCall(DMLabelMakeAllValid_Private(label));
1783:   }
1784:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1785:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1786:   /* Bcast name */
1787:   if (rank == 0) {
1788:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1789:     PetscCall(PetscStrlen(lname, &len));
1790:   }
1791:   nameSize = len;
1792:   PetscCallMPI(MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm));
1793:   PetscCall(PetscMalloc1(nameSize + 1, &name));
1794:   if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1795:   PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1796:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1797:   PetscCall(PetscFree(name));
1798:   /* Bcast defaultValue */
1799:   if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1800:   PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
1801:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1802:   PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
1803:   /* Determine received stratum values and initialise new label*/
1804:   PetscCall(PetscHSetICreate(&stratumHash));
1805:   PetscCall(PetscSectionGetStorageSize(leafSection, &size));
1806:   for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
1807:   PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
1808:   PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
1809:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1810:   PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
1811:   /* Turn leafStrata into indices rather than stratum values */
1812:   offset = 0;
1813:   PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
1814:   PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
1815:   for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
1816:   for (p = 0; p < size; ++p) {
1817:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1818:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1819:         leafStrata[p] = s;
1820:         break;
1821:       }
1822:     }
1823:   }
1824:   /* Rebuild the point strata on the receiver */
1825:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
1826:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1827:   for (p = pStart; p < pEnd; p++) {
1828:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1829:     PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1830:     for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1831:   }
1832:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
1833:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
1834:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
1835:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1836:     PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
1837:     PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s])));
1838:   }
1839:   /* Insert points into new strata */
1840:   PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
1841:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1842:   for (p = pStart; p < pEnd; p++) {
1843:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1844:     PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1845:     for (s = 0; s < dof; s++) {
1846:       stratum                               = leafStrata[offset + s];
1847:       points[stratum][strataIdx[stratum]++] = p;
1848:     }
1849:   }
1850:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1851:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &(points[s][0]), PETSC_OWN_POINTER, &((*labelNew)->points[s])));
1852:     PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
1853:   }
1854:   PetscCall(PetscFree(points));
1855:   PetscCall(PetscHSetIDestroy(&stratumHash));
1856:   PetscCall(PetscFree(leafStrata));
1857:   PetscCall(PetscFree(strataIdx));
1858:   PetscCall(PetscSectionDestroy(&leafSection));
1859:   PetscFunctionReturn(PETSC_SUCCESS);
1860: }

1862: /*@
1863:   DMLabelGather - Gather all label values from leafs into roots

1865:   Collective

1867:   Input Parameters:
1868: + label - the `DMLabel`
1869: - sf    - the `PetscSF` communication map

1871:   Output Parameter:
1872: . labelNew - the new `DMLabel` with localised leaf values

1874:   Level: developer

1876:   Note:
1877:   This is the inverse operation to `DMLabelDistribute()`.

1879: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
1880: @*/
1881: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1882: {
1883:   MPI_Comm        comm;
1884:   PetscSection    rootSection;
1885:   PetscSF         sfLabel;
1886:   PetscSFNode    *rootPoints, *leafPoints;
1887:   PetscInt        p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1888:   const PetscInt *rootDegree, *ilocal;
1889:   PetscInt       *rootStrata;
1890:   const char     *lname;
1891:   char           *name;
1892:   PetscInt        nameSize;
1893:   size_t          len = 0;
1894:   PetscMPIInt     rank, size;

1896:   PetscFunctionBegin;
1899:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1900:   PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1901:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1902:   PetscCallMPI(MPI_Comm_size(comm, &size));
1903:   /* Bcast name */
1904:   if (rank == 0) {
1905:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1906:     PetscCall(PetscStrlen(lname, &len));
1907:   }
1908:   nameSize = len;
1909:   PetscCallMPI(MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm));
1910:   PetscCall(PetscMalloc1(nameSize + 1, &name));
1911:   if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1912:   PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1913:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1914:   PetscCall(PetscFree(name));
1915:   /* Gather rank/index pairs of leaves into local roots to build
1916:      an inverse, multi-rooted SF. Note that this ignores local leaf
1917:      indexing due to the use of the multiSF in PetscSFGather. */
1918:   PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
1919:   PetscCall(PetscMalloc1(nroots, &leafPoints));
1920:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1921:   for (p = 0; p < nleaves; p++) {
1922:     PetscInt ilp = ilocal ? ilocal[p] : p;

1924:     leafPoints[ilp].index = ilp;
1925:     leafPoints[ilp].rank  = rank;
1926:   }
1927:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
1928:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
1929:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1930:   PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
1931:   PetscCall(PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints));
1932:   PetscCall(PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints));
1933:   PetscCall(PetscSFCreate(comm, &sfLabel));
1934:   PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
1935:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1936:   PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
1937:   /* Rebuild the point strata on the receiver */
1938:   for (p = 0, idx = 0; p < nroots; p++) {
1939:     for (d = 0; d < rootDegree[p]; d++) {
1940:       PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
1941:       PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
1942:       for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
1943:     }
1944:     idx += rootDegree[p];
1945:   }
1946:   PetscCall(PetscFree(leafPoints));
1947:   PetscCall(PetscFree(rootStrata));
1948:   PetscCall(PetscSectionDestroy(&rootSection));
1949:   PetscCall(PetscSFDestroy(&sfLabel));
1950:   PetscFunctionReturn(PETSC_SUCCESS);
1951: }

1953: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
1954: {
1955:   const PetscInt *degree;
1956:   const PetscInt *points;
1957:   PetscInt        Nr, r, Nl, l, val, defVal;

1959:   PetscFunctionBegin;
1960:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
1961:   /* Add in leaves */
1962:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
1963:   for (l = 0; l < Nl; ++l) {
1964:     PetscCall(DMLabelGetValue(label, points[l], &val));
1965:     if (val != defVal) valArray[points[l]] = val;
1966:   }
1967:   /* Add in shared roots */
1968:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &degree));
1969:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &degree));
1970:   for (r = 0; r < Nr; ++r) {
1971:     if (degree[r]) {
1972:       PetscCall(DMLabelGetValue(label, r, &val));
1973:       if (val != defVal) valArray[r] = val;
1974:     }
1975:   }
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }

1979: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1980: {
1981:   const PetscInt *degree;
1982:   const PetscInt *points;
1983:   PetscInt        Nr, r, Nl, l, val, defVal;

1985:   PetscFunctionBegin;
1986:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
1987:   /* Read out leaves */
1988:   PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
1989:   for (l = 0; l < Nl; ++l) {
1990:     const PetscInt p    = points[l];
1991:     const PetscInt cval = valArray[p];

1993:     if (cval != defVal) {
1994:       PetscCall(DMLabelGetValue(label, p, &val));
1995:       if (val == defVal) {
1996:         PetscCall(DMLabelSetValue(label, p, cval));
1997:         if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
1998:       }
1999:     }
2000:   }
2001:   /* Read out shared roots */
2002:   PetscCall(PetscSFComputeDegreeBegin(pointSF, &degree));
2003:   PetscCall(PetscSFComputeDegreeEnd(pointSF, &degree));
2004:   for (r = 0; r < Nr; ++r) {
2005:     if (degree[r]) {
2006:       const PetscInt cval = valArray[r];

2008:       if (cval != defVal) {
2009:         PetscCall(DMLabelGetValue(label, r, &val));
2010:         if (val == defVal) {
2011:           PetscCall(DMLabelSetValue(label, r, cval));
2012:           if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2013:         }
2014:       }
2015:     }
2016:   }
2017:   PetscFunctionReturn(PETSC_SUCCESS);
2018: }

2020: /*@
2021:   DMLabelPropagateBegin - Setup a cycle of label propagation

2023:   Collective

2025:   Input Parameters:
2026: + label - The `DMLabel` to propagate across processes
2027: - sf    - The `PetscSF` describing parallel layout of the label points

2029:   Level: intermediate

2031: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2032: @*/
2033: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2034: {
2035:   PetscInt    Nr, r, defVal;
2036:   PetscMPIInt size;

2038:   PetscFunctionBegin;
2039:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2040:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2041:   if (size > 1) {
2042:     PetscCall(DMLabelGetDefaultValue(label, &defVal));
2043:     PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2044:     if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2045:     for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2046:   }
2047:   PetscFunctionReturn(PETSC_SUCCESS);
2048: }

2050: /*@
2051:   DMLabelPropagateEnd - Tear down a cycle of label propagation

2053:   Collective

2055:   Input Parameters:
2056: + label   - The `DMLabel` to propagate across processes
2057: - pointSF - The `PetscSF` describing parallel layout of the label points

2059:   Level: intermediate

2061: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2062: @*/
2063: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2064: {
2065:   PetscFunctionBegin;
2066:   PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2067:   PetscCall(PetscFree(label->propArray));
2068:   label->propArray = NULL;
2069:   PetscFunctionReturn(PETSC_SUCCESS);
2070: }

2072: /*@C
2073:   DMLabelPropagatePush - Tear down a cycle of label propagation

2075:   Collective

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

2083:   Calling sequence of `markPoint`:
2084: + label - The `DMLabel`
2085: . p     - The point being marked
2086: . val   - The label value for `p`
2087: - ctx   - An optional user context

2089:   Level: intermediate

2091: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2092: @*/
2093: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx)
2094: {
2095:   PetscInt   *valArray = label->propArray, Nr;
2096:   PetscMPIInt size;

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

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

2111:        In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2112:        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
2113:        edge to the queue.
2114:     */
2115:     PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2116:     PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2117:     PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2118:     PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2119:     PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2120:     PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2121:   }
2122:   PetscFunctionReturn(PETSC_SUCCESS);
2123: }

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

2128:   Not Collective

2130:   Input Parameter:
2131: . label - the `DMLabel`

2133:   Output Parameters:
2134: + section - the section giving offsets for each stratum
2135: - is      - An `IS` containing all the label points

2137:   Level: developer

2139: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2140: @*/
2141: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2142: {
2143:   IS              vIS;
2144:   const PetscInt *values;
2145:   PetscInt       *points;
2146:   PetscInt        nV, vS = 0, vE = 0, v, N;

2148:   PetscFunctionBegin;
2150:   PetscCall(DMLabelGetNumValues(label, &nV));
2151:   PetscCall(DMLabelGetValueIS(label, &vIS));
2152:   PetscCall(ISGetIndices(vIS, &values));
2153:   if (nV) {
2154:     vS = values[0];
2155:     vE = values[0] + 1;
2156:   }
2157:   for (v = 1; v < nV; ++v) {
2158:     vS = PetscMin(vS, values[v]);
2159:     vE = PetscMax(vE, values[v] + 1);
2160:   }
2161:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2162:   PetscCall(PetscSectionSetChart(*section, vS, vE));
2163:   for (v = 0; v < nV; ++v) {
2164:     PetscInt n;

2166:     PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2167:     PetscCall(PetscSectionSetDof(*section, values[v], n));
2168:   }
2169:   PetscCall(PetscSectionSetUp(*section));
2170:   PetscCall(PetscSectionGetStorageSize(*section, &N));
2171:   PetscCall(PetscMalloc1(N, &points));
2172:   for (v = 0; v < nV; ++v) {
2173:     IS              is;
2174:     const PetscInt *spoints;
2175:     PetscInt        dof, off, p;

2177:     PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2178:     PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2179:     PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2180:     PetscCall(ISGetIndices(is, &spoints));
2181:     for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2182:     PetscCall(ISRestoreIndices(is, &spoints));
2183:     PetscCall(ISDestroy(&is));
2184:   }
2185:   PetscCall(ISRestoreIndices(vIS, &values));
2186:   PetscCall(ISDestroy(&vIS));
2187:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2188:   PetscFunctionReturn(PETSC_SUCCESS);
2189: }

2191: /*@C
2192:   DMLabelRegister - Adds a new label component implementation

2194:   Not Collective

2196:   Input Parameters:
2197: + name        - The name of a new user-defined creation routine
2198: - create_func - The creation routine itself

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

2203:   Example Usage:
2204: .vb
2205:   DMLabelRegister("my_label", MyLabelCreate);
2206: .ve

2208:   Then, your label type can be chosen with the procedural interface via
2209: .vb
2210:   DMLabelCreate(MPI_Comm, DMLabel *);
2211:   DMLabelSetType(DMLabel, "my_label");
2212: .ve
2213:   or at runtime via the option
2214: .vb
2215:   -dm_label_type my_label
2216: .ve

2218:   Level: advanced

2220: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2221: @*/
2222: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2223: {
2224:   PetscFunctionBegin;
2225:   PetscCall(DMInitializePackage());
2226:   PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2227:   PetscFunctionReturn(PETSC_SUCCESS);
2228: }

2230: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2231: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);

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

2236:   Not Collective

2238:   Level: advanced

2240: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2241: @*/
2242: PetscErrorCode DMLabelRegisterAll(void)
2243: {
2244:   PetscFunctionBegin;
2245:   if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2246:   DMLabelRegisterAllCalled = PETSC_TRUE;

2248:   PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2249:   PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2250:   PetscFunctionReturn(PETSC_SUCCESS);
2251: }

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

2256:   Level: developer

2258: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2259: @*/
2260: PetscErrorCode DMLabelRegisterDestroy(void)
2261: {
2262:   PetscFunctionBegin;
2263:   PetscCall(PetscFunctionListDestroy(&DMLabelList));
2264:   DMLabelRegisterAllCalled = PETSC_FALSE;
2265:   PetscFunctionReturn(PETSC_SUCCESS);
2266: }

2268: /*@C
2269:   DMLabelSetType - Sets the particular implementation for a label.

2271:   Collective

2273:   Input Parameters:
2274: + label  - The label
2275: - method - The name of the label type

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

2280:   Level: intermediate

2282: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2283: @*/
2284: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2285: {
2286:   PetscErrorCode (*r)(DMLabel);
2287:   PetscBool match;

2289:   PetscFunctionBegin;
2291:   PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2292:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

2294:   PetscCall(DMLabelRegisterAll());
2295:   PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2296:   PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);

2298:   PetscTryTypeMethod(label, destroy);
2299:   PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2300:   PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2301:   PetscCall((*r)(label));
2302:   PetscFunctionReturn(PETSC_SUCCESS);
2303: }

2305: /*@C
2306:   DMLabelGetType - Gets the type name (as a string) from the label.

2308:   Not Collective

2310:   Input Parameter:
2311: . label - The `DMLabel`

2313:   Output Parameter:
2314: . type - The `DMLabel` type name

2316:   Level: intermediate

2318: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2319: @*/
2320: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2321: {
2322:   PetscFunctionBegin;
2324:   PetscAssertPointer(type, 2);
2325:   PetscCall(DMLabelRegisterAll());
2326:   *type = ((PetscObject)label)->type_name;
2327:   PetscFunctionReturn(PETSC_SUCCESS);
2328: }

2330: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2331: {
2332:   PetscFunctionBegin;
2333:   label->ops->view         = DMLabelView_Concrete;
2334:   label->ops->setup        = NULL;
2335:   label->ops->duplicate    = DMLabelDuplicate_Concrete;
2336:   label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2337:   PetscFunctionReturn(PETSC_SUCCESS);
2338: }

2340: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2341: {
2342:   PetscFunctionBegin;
2344:   PetscCall(DMLabelInitialize_Concrete(label));
2345:   PetscFunctionReturn(PETSC_SUCCESS);
2346: }

2348: /*@
2349:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2350:   the local section and an `PetscSF` describing the section point overlap.

2352:   Collective

2354:   Input Parameters:
2355: + s                  - The `PetscSection` for the local field layout
2356: . sf                 - The `PetscSF` describing parallel layout of the section points
2357: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2358: . label              - The label specifying the points
2359: - labelValue         - The label stratum specifying the points

2361:   Output Parameter:
2362: . gsection - The `PetscSection` for the global field layout

2364:   Level: developer

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

2369: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2370: @*/
2371: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2372: {
2373:   PetscInt *neg = NULL, *tmpOff = NULL;
2374:   PetscInt  pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

2376:   PetscFunctionBegin;
2380:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2381:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2382:   PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2383:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2384:   if (nroots >= 0) {
2385:     PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2386:     PetscCall(PetscCalloc1(nroots, &neg));
2387:     if (nroots > pEnd - pStart) {
2388:       PetscCall(PetscCalloc1(nroots, &tmpOff));
2389:     } else {
2390:       tmpOff = &(*gsection)->atlasDof[-pStart];
2391:     }
2392:   }
2393:   /* Mark ghost points with negative dof */
2394:   for (p = pStart; p < pEnd; ++p) {
2395:     PetscInt value;

2397:     PetscCall(DMLabelGetValue(label, p, &value));
2398:     if (value != labelValue) continue;
2399:     PetscCall(PetscSectionGetDof(s, p, &dof));
2400:     PetscCall(PetscSectionSetDof(*gsection, p, dof));
2401:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2402:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2403:     if (neg) neg[p] = -(dof + 1);
2404:   }
2405:   PetscCall(PetscSectionSetUpBC(*gsection));
2406:   if (nroots >= 0) {
2407:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2408:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2409:     if (nroots > pEnd - pStart) {
2410:       for (p = pStart; p < pEnd; ++p) {
2411:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2412:       }
2413:     }
2414:   }
2415:   /* Calculate new sizes, get process offset, and calculate point offsets */
2416:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2417:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2418:     (*gsection)->atlasOff[p] = off;
2419:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2420:   }
2421:   PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2422:   globalOff -= off;
2423:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2424:     (*gsection)->atlasOff[p] += globalOff;
2425:     if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2426:   }
2427:   /* Put in negative offsets for ghost points */
2428:   if (nroots >= 0) {
2429:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2430:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2431:     if (nroots > pEnd - pStart) {
2432:       for (p = pStart; p < pEnd; ++p) {
2433:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2434:       }
2435:     }
2436:   }
2437:   if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2438:   PetscCall(PetscFree(neg));
2439:   PetscFunctionReturn(PETSC_SUCCESS);
2440: }

2442: typedef struct _n_PetscSectionSym_Label {
2443:   DMLabel              label;
2444:   PetscCopyMode       *modes;
2445:   PetscInt            *sizes;
2446:   const PetscInt    ***perms;
2447:   const PetscScalar ***rots;
2448:   PetscInt (*minMaxOrients)[2];
2449:   PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2450: } PetscSectionSym_Label;

2452: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2453: {
2454:   PetscInt               i, j;
2455:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;

2457:   PetscFunctionBegin;
2458:   for (i = 0; i <= sl->numStrata; i++) {
2459:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2460:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2461:         if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2462:         if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2463:       }
2464:       if (sl->perms[i]) {
2465:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

2467:         PetscCall(PetscFree(perms));
2468:       }
2469:       if (sl->rots[i]) {
2470:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

2472:         PetscCall(PetscFree(rots));
2473:       }
2474:     }
2475:   }
2476:   PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2477:   PetscCall(DMLabelDestroy(&sl->label));
2478:   sl->numStrata = 0;
2479:   PetscFunctionReturn(PETSC_SUCCESS);
2480: }

2482: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2483: {
2484:   PetscFunctionBegin;
2485:   PetscCall(PetscSectionSymLabelReset(sym));
2486:   PetscCall(PetscFree(sym->data));
2487:   PetscFunctionReturn(PETSC_SUCCESS);
2488: }

2490: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2491: {
2492:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2493:   PetscBool              isAscii;
2494:   DMLabel                label = sl->label;
2495:   const char            *name;

2497:   PetscFunctionBegin;
2498:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2499:   if (isAscii) {
2500:     PetscInt          i, j, k;
2501:     PetscViewerFormat format;

2503:     PetscCall(PetscViewerGetFormat(viewer, &format));
2504:     if (label) {
2505:       PetscCall(PetscViewerGetFormat(viewer, &format));
2506:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2507:         PetscCall(PetscViewerASCIIPushTab(viewer));
2508:         PetscCall(DMLabelView(label, viewer));
2509:         PetscCall(PetscViewerASCIIPopTab(viewer));
2510:       } else {
2511:         PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2512:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Label '%s'\n", name));
2513:       }
2514:     } else {
2515:       PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2516:     }
2517:     PetscCall(PetscViewerASCIIPushTab(viewer));
2518:     for (i = 0; i <= sl->numStrata; i++) {
2519:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

2521:       if (!(sl->perms[i] || sl->rots[i])) {
2522:         PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2523:       } else {
2524:         PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2525:         PetscCall(PetscViewerASCIIPushTab(viewer));
2526:         PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2527:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2528:           PetscCall(PetscViewerASCIIPushTab(viewer));
2529:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2530:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2531:               PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2532:             } else {
2533:               PetscInt tab;

2535:               PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2536:               PetscCall(PetscViewerASCIIPushTab(viewer));
2537:               PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2538:               if (sl->perms[i] && sl->perms[i][j]) {
2539:                 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2540:                 PetscCall(PetscViewerASCIISetTab(viewer, 0));
2541:                 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2542:                 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2543:                 PetscCall(PetscViewerASCIISetTab(viewer, tab));
2544:               }
2545:               if (sl->rots[i] && sl->rots[i][j]) {
2546:                 PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations:  "));
2547:                 PetscCall(PetscViewerASCIISetTab(viewer, 0));
2548: #if defined(PETSC_USE_COMPLEX)
2549:                 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])));
2550: #else
2551:                 for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2552: #endif
2553:                 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2554:                 PetscCall(PetscViewerASCIISetTab(viewer, tab));
2555:               }
2556:               PetscCall(PetscViewerASCIIPopTab(viewer));
2557:             }
2558:           }
2559:           PetscCall(PetscViewerASCIIPopTab(viewer));
2560:         }
2561:         PetscCall(PetscViewerASCIIPopTab(viewer));
2562:       }
2563:     }
2564:     PetscCall(PetscViewerASCIIPopTab(viewer));
2565:   }
2566:   PetscFunctionReturn(PETSC_SUCCESS);
2567: }

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

2572:   Logically

2574:   Input Parameters:
2575: + sym   - the section symmetries
2576: - label - the `DMLabel` describing the types of points

2578:   Level: developer:

2580: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2581: @*/
2582: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2583: {
2584:   PetscSectionSym_Label *sl;

2586:   PetscFunctionBegin;
2588:   sl = (PetscSectionSym_Label *)sym->data;
2589:   if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2590:   if (label) {
2591:     sl->label = label;
2592:     PetscCall(PetscObjectReference((PetscObject)label));
2593:     PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2594:     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));
2595:     PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2596:     PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2597:     PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2598:     PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2599:     PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2600:   }
2601:   PetscFunctionReturn(PETSC_SUCCESS);
2602: }

2604: /*@C
2605:   PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum

2607:   Logically Collective

2609:   Input Parameters:
2610: + sym     - the section symmetries
2611: - stratum - the stratum value in the label that we are assigning symmetries for

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

2620:   Level: developer

2622: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2623: @*/
2624: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2625: {
2626:   PetscSectionSym_Label *sl;
2627:   const char            *name;
2628:   PetscInt               i;

2630:   PetscFunctionBegin;
2632:   sl = (PetscSectionSym_Label *)sym->data;
2633:   PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2634:   for (i = 0; i <= sl->numStrata; i++) {
2635:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2637:     if (stratum == value) break;
2638:   }
2639:   PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2640:   PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2641:   if (size) {
2642:     PetscAssertPointer(size, 3);
2643:     *size = sl->sizes[i];
2644:   }
2645:   if (minOrient) {
2646:     PetscAssertPointer(minOrient, 4);
2647:     *minOrient = sl->minMaxOrients[i][0];
2648:   }
2649:   if (maxOrient) {
2650:     PetscAssertPointer(maxOrient, 5);
2651:     *maxOrient = sl->minMaxOrients[i][1];
2652:   }
2653:   if (perms) {
2654:     PetscAssertPointer(perms, 6);
2655:     *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;
2656:   }
2657:   if (rots) {
2658:     PetscAssertPointer(rots, 7);
2659:     *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;
2660:   }
2661:   PetscFunctionReturn(PETSC_SUCCESS);
2662: }

2664: /*@C
2665:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

2667:   Logically

2669:   Input Parameters:
2670: + sym       - the section symmetries
2671: . stratum   - the stratum value in the label that we are assigning symmetries for
2672: . size      - the number of dofs for points in the `stratum` of the label
2673: . minOrient - the smallest orientation for a point in this `stratum`
2674: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2675: . mode      - how `sym` should copy the `perms` and `rots` arrays
2676: . perms     - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation.  A `NULL` permutation is the identity
2677: - 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

2679:   Level: developer

2681: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2682: @*/
2683: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2684: {
2685:   PetscSectionSym_Label *sl;
2686:   const char            *name;
2687:   PetscInt               i, j, k;

2689:   PetscFunctionBegin;
2691:   sl = (PetscSectionSym_Label *)sym->data;
2692:   PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2693:   for (i = 0; i <= sl->numStrata; i++) {
2694:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2696:     if (stratum == value) break;
2697:   }
2698:   PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2699:   PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2700:   sl->sizes[i]            = size;
2701:   sl->modes[i]            = mode;
2702:   sl->minMaxOrients[i][0] = minOrient;
2703:   sl->minMaxOrients[i][1] = maxOrient;
2704:   if (mode == PETSC_COPY_VALUES) {
2705:     if (perms) {
2706:       PetscInt **ownPerms;

2708:       PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2709:       for (j = 0; j < maxOrient - minOrient; j++) {
2710:         if (perms[j]) {
2711:           PetscCall(PetscMalloc1(size, &ownPerms[j]));
2712:           for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2713:         }
2714:       }
2715:       sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2716:     }
2717:     if (rots) {
2718:       PetscScalar **ownRots;

2720:       PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2721:       for (j = 0; j < maxOrient - minOrient; j++) {
2722:         if (rots[j]) {
2723:           PetscCall(PetscMalloc1(size, &ownRots[j]));
2724:           for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2725:         }
2726:       }
2727:       sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2728:     }
2729:   } else {
2730:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2731:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2732:   }
2733:   PetscFunctionReturn(PETSC_SUCCESS);
2734: }

2736: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2737: {
2738:   PetscInt               i, j, numStrata;
2739:   PetscSectionSym_Label *sl;
2740:   DMLabel                label;

2742:   PetscFunctionBegin;
2743:   sl        = (PetscSectionSym_Label *)sym->data;
2744:   numStrata = sl->numStrata;
2745:   label     = sl->label;
2746:   for (i = 0; i < numPoints; i++) {
2747:     PetscInt point = points[2 * i];
2748:     PetscInt ornt  = points[2 * i + 1];

2750:     for (j = 0; j < numStrata; j++) {
2751:       if (label->validIS[j]) {
2752:         PetscInt k;

2754:         PetscCall(ISLocate(label->points[j], point, &k));
2755:         if (k >= 0) break;
2756:       } else {
2757:         PetscBool has;

2759:         PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2760:         if (has) break;
2761:       }
2762:     }
2763:     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],
2764:                j < numStrata ? label->stratumValues[j] : label->defaultValue);
2765:     if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2766:     if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2767:   }
2768:   PetscFunctionReturn(PETSC_SUCCESS);
2769: }

2771: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2772: {
2773:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2774:   IS                     valIS;
2775:   const PetscInt        *values;
2776:   PetscInt               Nv, v;

2778:   PetscFunctionBegin;
2779:   PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2780:   PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2781:   PetscCall(ISGetIndices(valIS, &values));
2782:   for (v = 0; v < Nv; ++v) {
2783:     const PetscInt      val = values[v];
2784:     PetscInt            size, minOrient, maxOrient;
2785:     const PetscInt    **perms;
2786:     const PetscScalar **rots;

2788:     PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2789:     PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
2790:   }
2791:   PetscCall(ISDestroy(&valIS));
2792:   PetscFunctionReturn(PETSC_SUCCESS);
2793: }

2795: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2796: {
2797:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2798:   DMLabel                dlabel;

2800:   PetscFunctionBegin;
2801:   PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
2802:   PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
2803:   PetscCall(DMLabelDestroy(&dlabel));
2804:   PetscCall(PetscSectionSymCopy(sym, *dsym));
2805:   PetscFunctionReturn(PETSC_SUCCESS);
2806: }

2808: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2809: {
2810:   PetscSectionSym_Label *sl;

2812:   PetscFunctionBegin;
2813:   PetscCall(PetscNew(&sl));
2814:   sym->ops->getpoints  = PetscSectionSymGetPoints_Label;
2815:   sym->ops->distribute = PetscSectionSymDistribute_Label;
2816:   sym->ops->copy       = PetscSectionSymCopy_Label;
2817:   sym->ops->view       = PetscSectionSymView_Label;
2818:   sym->ops->destroy    = PetscSectionSymDestroy_Label;
2819:   sym->data            = (void *)sl;
2820:   PetscFunctionReturn(PETSC_SUCCESS);
2821: }

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

2826:   Collective

2828:   Input Parameters:
2829: + comm  - the MPI communicator for the new symmetry
2830: - label - the label defining the strata

2832:   Output Parameter:
2833: . sym - the section symmetries

2835:   Level: developer

2837: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2838: @*/
2839: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2840: {
2841:   PetscFunctionBegin;
2842:   PetscCall(DMInitializePackage());
2843:   PetscCall(PetscSectionSymCreate(comm, sym));
2844:   PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
2845:   PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
2846:   PetscFunctionReturn(PETSC_SUCCESS);
2847: }