Actual source code: dmlabel.c

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

  7: PetscFunctionList DMLabelList              = NULL;
  8: PetscBool         DMLabelRegisterAllCalled = PETSC_FALSE;

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

 13:   Collective

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

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

 22:   Level: beginner

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

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

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

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

 56:   Collective

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

 61:   Level: intermediate

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

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

 76:   Not collective

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

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

 85:   Level: developer

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

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

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

126:   Not Collective

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

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

134:   Level: developer

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

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

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

150:   Not Collective

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

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

159:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

313:   Level: beginner

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

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

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

331:   Not Collective

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

338:   Level: beginner

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

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

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

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

393:   Not Collective

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

399:   Level: beginner

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

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

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

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

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

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

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

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

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

460:   Collective

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

466:   Level: intermediate

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

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

484:   Not Collective

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

489:   Level: beginner

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

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

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

519:   Collective

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

524:   Level: beginner

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

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

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

559:   Collective

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

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

567:   Level: intermediate

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

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

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

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

604:   Collective; No Fortran Support

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

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

615:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

714:   Not Collective

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

719:   Level: intermediate

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

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

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

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

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

752:   Not Collective

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

759:   Level: intermediate

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

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

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

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

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

797:   Not Collective

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

802:   Level: intermediate

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

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

819:   Not Collective

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

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

828:   Level: intermediate

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

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

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

854:   Not Collective

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

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

863:   Level: developer

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

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

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

882:   Not Collective

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

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

891:   Level: developer

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

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

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

915:   Not Collective

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

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

925:   Level: intermediate

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

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

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

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

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

966:   Not Collective

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

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

974:   Level: beginner

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

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

990:   Not Collective

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

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

998:   Level: beginner

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

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

1014:   Not Collective

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

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

1023:   Level: intermediate

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

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

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

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

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

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

1068:   Not Collective

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

1075:   Level: intermediate

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

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

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

1098:   Not Collective

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

1105:   Level: intermediate

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

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

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

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

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

1134:   Not Collective

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

1141:   Level: intermediate

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

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

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

1169:   Not Collective

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

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

1177:   Level: intermediate

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

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

1193:   Not Collective

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

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

1201:   Level: intermediate

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

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

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

1222:   Not Collective

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

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

1231:   Level: intermediate

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

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

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

1259:   Not Collective

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

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

1267:   Level: intermediate

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

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

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

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

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

1302:   Not Collective

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

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

1311:   Level: intermediate

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

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

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

1333:   Not Collective

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

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

1342:   Level: intermediate

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

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

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

1361:   Not Collective

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

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

1370:   Level: intermediate

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

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

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

1389:   Not Collective

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

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

1399:   Level: intermediate

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

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

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

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

1441:   Not Collective

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

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

1450:   Level: intermediate

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

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

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

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

1476:   Not Collective

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

1483:   Level: intermediate

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

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

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

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

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

1522:   Not Collective

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

1528:   Level: intermediate

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

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

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

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

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

1569:   Not Collective

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

1577:   Level: intermediate

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

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

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

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

1598:   Not Collective

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

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

1608:   Level: intermediate

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

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

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

1633:   Not Collective

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

1640:   Level: intermediate

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

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

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

1664:   Not Collective

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

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

1673:   Level: intermediate

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

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

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

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

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

1726:   Not collective

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

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

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

1738:   Level: intermediate

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

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

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

1764:   Not collective

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

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

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

1776:   Level: intermediate

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

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

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

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

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

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

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

1862:   Collective

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

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

1871:   Level: intermediate

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

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

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

1977:   Collective

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

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

1986:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

2135:   Collective

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

2141:   Level: intermediate

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

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

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

2165:   Collective

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

2171:   Level: intermediate

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

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

2187:   Collective

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

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

2201:   Level: intermediate

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

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

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

2223:        In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2224:        values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new
2225:        edge to the queue.
2226:     */
2227:     PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2228:     PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2229:     PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2230:     PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2231:     PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2232:     PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2233:   }
2234:   PetscFunctionReturn(PETSC_SUCCESS);
2235: }

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

2240:   Not Collective

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

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

2249:   Level: developer

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

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

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

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

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

2306:   Not Collective

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

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

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

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

2330:   Level: advanced

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

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

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

2348:   Not Collective

2350:   Level: advanced

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

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

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

2368:   Level: developer

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

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

2383:   Collective

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

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

2392:   Level: intermediate

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

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

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

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

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

2420:   Not Collective

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

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

2428:   Level: intermediate

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

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

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

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

2464:   Collective

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

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

2476:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2684:   Logically

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

2690:   Level: developer:

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

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

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

2719:   Logically Collective

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

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

2732:   Level: developer

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

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

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

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

2779:   Logically

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

2791:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2938:   Collective

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

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

2947:   Level: developer

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