Actual source code: dmlabel.c

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

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

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

 13:   Collective

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

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

 22:   Level: beginner

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

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

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

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

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

 57:   Collective

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

 62:   Level: intermediate

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

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

 77:   Not collective

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

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

 86:   Level: developer

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

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

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

126:   Not Collective

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

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

134:   Level: developer

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

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

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

150:   Not Collective

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

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

159:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

313:   Level: beginner

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

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

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

331:   Not Collective

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

338:   Level: beginner

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

346:   PetscFunctionBegin;
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: PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
448: {
449:   PetscBool iascii;

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

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

460:   Collective

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

466:   Level: intermediate

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

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

484:   Not Collective

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

489:   Level: beginner

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

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

521: /*@
522:   DMLabelDestroy - Destroys a `DMLabel`

524:   Collective

526:   Input Parameter:
527: . label - The `DMLabel`

529:   Level: beginner

531: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
532: @*/
533: PetscErrorCode DMLabelDestroy(DMLabel *label)
534: {
535:   PetscFunctionBegin;
536:   if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
538:   if (--((PetscObject)(*label))->refct > 0) {
539:     *label = NULL;
540:     PetscFunctionReturn(PETSC_SUCCESS);
541:   }
542:   PetscCall(DMLabelReset(*label));
543:   PetscCall(PetscHMapIDestroy(&(*label)->hmap));
544:   PetscCall(PetscHeaderDestroy(label));
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

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

561: /*@
562:   DMLabelDuplicate - Duplicates a `DMLabel`

564:   Collective

566:   Input Parameter:
567: . label - The `DMLabel`

569:   Output Parameter:
570: . labelnew - new label

572:   Level: intermediate

574: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
575: @*/
576: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
577: {
578:   const char *name;

580:   PetscFunctionBegin;
582:   PetscCall(DMLabelMakeAllValid_Private(label));
583:   PetscCall(PetscObjectGetName((PetscObject)label, &name));
584:   PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));

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

605: /*@C
606:   DMLabelCompare - Compare two `DMLabel` objects

608:   Collective; No Fortran Support

610:   Input Parameters:
611: + comm - Comm over which to compare labels
612: . l0 - First `DMLabel`
613: - l1 - Second `DMLabel`

615:   Output Parameters
616: + equal   - (Optional) Flag whether the two labels are equal
617: - message - (Optional) Message describing the difference

619:   Level: intermediate

621:   Notes:
622:   The output flag equal is the same on all processes.
623:   If it is passed as `NULL` and difference is found, an error is thrown on all processes.
624:   Make sure to pass `NULL` on all processes.

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

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

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

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

644:   PetscFunctionBegin;
649:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
650:   PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
651:   PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
652:   {
653:     PetscInt v0, v1;

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

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

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

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

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

715:   Not Collective

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

720:   Level: intermediate

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

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

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

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

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

753:   Not Collective

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

760:   Level: intermediate

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

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

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

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

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

798:   Not Collective

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

803:   Level: intermediate

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

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

820:   Not Collective

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

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

829:   Level: intermediate

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

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

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

855:   Not Collective

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

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

864:   Level: developer

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

872:   PetscFunctionBegin;
875:   PetscCall(DMLabelLookupStratum(label, value, &v));
876:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
877:   PetscFunctionReturn(PETSC_SUCCESS);
878: }

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

883:   Not Collective

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

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

892:   Level: developer

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

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

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

916:   Not Collective

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

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

926:   Level: intermediate

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

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

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

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

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

967:   Not Collective

969:   Input parameter:
970: . label - a `DMLabel` object

972:   Output parameter:
973: . defaultValue - the default value

975:   Level: beginner

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

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

991:   Not Collective

993:   Input parameter:
994: . label - a `DMLabel` object

996:   Output parameter:
997: . defaultValue - the default value

999:   Level: beginner

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

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

1015:   Not Collective

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

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

1024:   Level: intermediate

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

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

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

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

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

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

1069:   Not Collective

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

1076:   Level: intermediate

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

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

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

1099:   Not Collective

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

1106:   Level: intermediate

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

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

1121:   if (label->bt) {
1122:     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);
1123:     PetscCall(PetscBTClear(label->bt, point - label->pStart));
1124:   }

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

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

1135:   Not Collective

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

1142:   Level: intermediate

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

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

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

1170:   Not Collective

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

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

1178:   Level: intermediate

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

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

1194:   Not Collective

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

1199:   Output Parameter:
1200: . is    - the value `IS`

1202:   Level: intermediate

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

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

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

1223:   Not Collective

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

1228:   Output Parameter:
1229: . is    - the value `IS`

1231:   Level: intermediate

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

1237: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1238: @*/
1239: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1240: {
1241:   PetscInt  i, j;
1242:   PetscInt *valuesArr;

1244:   PetscFunctionBegin;
1247:   PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1248:   for (i = 0, j = 0; i < label->numStrata; i++) {
1249:     PetscInt n;

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

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

1266:   Not Collective

1268:   Input Parameters:
1269: + label - the `DMLabel`
1270: - value - the value

1272:   Output Parameter:
1273: . index - the index of value in the list of values

1275:   Level: intermediate

1277: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1278: @*/
1279: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1280: {
1281:   PetscInt v;

1283:   PetscFunctionBegin;
1286:   /* Do not assume they are sorted */
1287:   for (v = 0; v < label->numStrata; ++v)
1288:     if (label->stratumValues[v] == value) break;
1289:   if (v >= label->numStrata) *index = -1;
1290:   else *index = v;
1291:   PetscFunctionReturn(PETSC_SUCCESS);
1292: }

1294: /*@
1295:   DMLabelHasStratum - Determine whether points exist with the given value

1297:   Not Collective

1299:   Input Parameters:
1300: + label - the `DMLabel`
1301: - value - the stratum value

1303:   Output Parameter:
1304: . exists - Flag saying whether points exist

1306:   Level: intermediate

1308: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1309: @*/
1310: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1311: {
1312:   PetscInt v;

1314:   PetscFunctionBegin;
1317:   PetscCall(DMLabelLookupStratum(label, value, &v));
1318:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1319:   PetscFunctionReturn(PETSC_SUCCESS);
1320: }

1322: /*@
1323:   DMLabelGetStratumSize - Get the size of a stratum

1325:   Not Collective

1327:   Input Parameters:
1328: + label - the `DMLabel`
1329: - value - the stratum value

1331:   Output Parameter:
1332: . size - The number of points in the stratum

1334:   Level: intermediate

1336: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1337: @*/
1338: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1339: {
1340:   PetscInt v;

1342:   PetscFunctionBegin;
1345:   PetscCall(DMLabelLookupStratum(label, value, &v));
1346:   PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1347:   PetscFunctionReturn(PETSC_SUCCESS);
1348: }

1350: /*@
1351:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1353:   Not Collective

1355:   Input Parameters:
1356: + label - the `DMLabel`
1357: - value - the stratum value

1359:   Output Parameters:
1360: + start - the smallest point in the stratum
1361: - end - the largest point in the stratum

1363:   Level: intermediate

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

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

1394: PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1395: {
1396:   PetscFunctionBegin;
1397:   PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1398:   *pointIS = label->points[v];
1399:   PetscFunctionReturn(PETSC_SUCCESS);
1400: }

1402: /*@
1403:   DMLabelGetStratumIS - Get an `IS` with the stratum points

1405:   Not Collective

1407:   Input Parameters:
1408: + label - the `DMLabel`
1409: - value - the stratum value

1411:   Output Parameter:
1412: . points - The stratum points

1414:   Level: intermediate

1416:   Notes:
1417:   The output `IS` should be destroyed when no longer needed.
1418:   Returns `NULL` if the stratum is empty.

1420: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1421: @*/
1422: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1423: {
1424:   PetscInt v;

1426:   PetscFunctionBegin;
1429:   *points = NULL;
1430:   PetscCall(DMLabelLookupStratum(label, value, &v));
1431:   if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1432:   PetscCall(DMLabelMakeValid_Private(label, v));
1433:   PetscUseTypeMethod(label, getstratumis, v, points);
1434:   PetscFunctionReturn(PETSC_SUCCESS);
1435: }

1437: /*@
1438:   DMLabelSetStratumIS - Set the stratum points using an `IS`

1440:   Not Collective

1442:   Input Parameters:
1443: + label - the `DMLabel`
1444: . value - the stratum value
1445: - points - The stratum points

1447:   Level: intermediate

1449: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1450: @*/
1451: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1452: {
1453:   PetscInt v;

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

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

1476:       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);
1477:       PetscCall(PetscBTSet(label->bt, point - label->pStart));
1478:     }
1479:   }
1480:   PetscFunctionReturn(PETSC_SUCCESS);
1481: }

1483: /*@
1484:   DMLabelClearStratum - Remove a stratum

1486:   Not Collective

1488:   Input Parameters:
1489: + label - the `DMLabel`
1490: - value - the stratum value

1492:   Level: intermediate

1494: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1495: @*/
1496: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1497: {
1498:   PetscInt v;

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

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

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

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

1533:   Not Collective

1535:   Input Parameters:
1536: + label  - The `DMLabel`
1537: . value  - The label value for all points
1538: . pStart - The first point
1539: - pEnd   - A point beyond all marked points

1541:   Level: intermediate

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

1546: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1547: @*/
1548: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1549: {
1550:   IS pIS;

1552:   PetscFunctionBegin;
1553:   PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1554:   PetscCall(DMLabelSetStratumIS(label, value, pIS));
1555:   PetscCall(ISDestroy(&pIS));
1556:   PetscFunctionReturn(PETSC_SUCCESS);
1557: }

1559: /*@
1560:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

1562:   Not Collective

1564:   Input Parameters:
1565: + label  - The `DMLabel`
1566: . value  - The label value
1567: - p      - A point with this value

1569:   Output Parameter:
1570: . 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

1572:   Level: intermediate

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

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

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

1600:   Not Collective

1602:   Input Parameters:
1603: + label - the `DMLabel`
1604: . start - the first point kept
1605: - end - one more than the last point kept

1607:   Level: intermediate

1609: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1610: @*/
1611: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1612: {
1613:   PetscInt v;

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

1628: /*@
1629:   DMLabelPermute - Create a new label with permuted points

1631:   Not Collective

1633:   Input Parameters:
1634: + label - the `DMLabel`
1635: - permutation - the point permutation

1637:   Output Parameter:
1638: . labelnew - the new label containing the permuted points

1640:   Level: intermediate

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

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

1663:     PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1664:     PetscCall(PetscCalloc1(size, &pointsNew));
1665:     for (q = 0; q < size; ++q) {
1666:       const PetscInt point = points[q];

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

1690: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1691: {
1692:   MPI_Comm     comm;
1693:   PetscInt     s, l, nroots, nleaves, offset, size;
1694:   PetscInt    *remoteOffsets, *rootStrata, *rootIdx;
1695:   PetscSection rootSection;
1696:   PetscSF      labelSF;

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

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

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

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

1754:   Collective

1756:   Input Parameters:
1757: + label - the `DMLabel`
1758: - sf    - the map from old to new distribution

1760:   Output Parameter:
1761: . labelnew - the new redistributed label

1763:   Level: intermediate

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

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

1866: /*@
1867:   DMLabelGather - Gather all label values from leafs into roots

1869:   Collective

1871:   Input Parameters:
1872: + label - the `DMLabel`
1873: - sf - the `PetscSF` communication map

1875:   Output Parameter:
1876: . labelNew - the new `DMLabel` with localised leaf values

1878:   Level: developer

1880:   Note:
1881:   This is the inverse operation to `DMLabelDistribute()`.

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

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

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

1957: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
1958: {
1959:   const PetscInt *degree;
1960:   const PetscInt *points;
1961:   PetscInt        Nr, r, Nl, l, val, defVal;

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

1983: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1984: {
1985:   const PetscInt *degree;
1986:   const PetscInt *points;
1987:   PetscInt        Nr, r, Nl, l, val, defVal;

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

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

2012:       if (cval != defVal) {
2013:         PetscCall(DMLabelGetValue(label, r, &val));
2014:         if (val == defVal) {
2015:           PetscCall(DMLabelSetValue(label, r, cval));
2016:           if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2017:         }
2018:       }
2019:     }
2020:   }
2021:   PetscFunctionReturn(PETSC_SUCCESS);
2022: }

2024: /*@
2025:   DMLabelPropagateBegin - Setup a cycle of label propagation

2027:   Collective

2029:   Input Parameters:
2030: + label - The `DMLabel` to propagate across processes
2031: - sf    - The `PetscSF` describing parallel layout of the label points

2033:   Level: intermediate

2035: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2036: @*/
2037: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2038: {
2039:   PetscInt    Nr, r, defVal;
2040:   PetscMPIInt size;

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

2054: /*@
2055:   DMLabelPropagateEnd - Tear down a cycle of label propagation

2057:   Collective

2059:   Input Parameters:
2060: + label - The `DMLabel` to propagate across processes
2061: - sf    - The `PetscSF` describing parallel layout of the label points

2063:   Level: intermediate

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

2076: /*@C
2077:   DMLabelPropagatePush - Tear down a cycle of label propagation

2079:   Collective

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

2087:   Calling sequence of `markPoint`:
2088: $ PetscErrorCode markPoint(DMLabel label, PetscInt p, PetscInt val, void *ctx);
2089: + label - The `DMLabel`
2090: . p     - The point being marked
2091: . val   - The label value for `p`
2092: - ctx   - An optional user context

2094:   Level: intermediate

2096: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2097: @*/
2098: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
2099: {
2100:   PetscInt   *valArray = label->propArray, Nr;
2101:   PetscMPIInt size;

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

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

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

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

2133:   Not Collective

2135:   Input Parameter:
2136: . label - the `DMLabel`

2138:   Output Parameters:
2139: + section - the section giving offsets for each stratum
2140: - is - An `IS` containing all the label points

2142:   Level: developer

2144: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2145: @*/
2146: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2147: {
2148:   IS              vIS;
2149:   const PetscInt *values;
2150:   PetscInt       *points;
2151:   PetscInt        nV, vS = 0, vE = 0, v, N;

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

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

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

2196: /*@C
2197:   DMLabelRegister - Adds a new label component implementation

2199:   Not Collective

2201:   Input Parameters:
2202: + name        - The name of a new user-defined creation routine
2203: - create_func - The creation routine itself

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

2208:   Sample usage:
2209: .vb
2210:   DMLabelRegister("my_label", MyLabelCreate);
2211: .ve

2213:   Then, your label type can be chosen with the procedural interface via
2214: .vb
2215:   DMLabelCreate(MPI_Comm, DMLabel *);
2216:   DMLabelSetType(DMLabel, "my_label");
2217: .ve
2218:   or at runtime via the option
2219: .vb
2220:   -dm_label_type my_label
2221: .ve

2223:   Level: advanced

2225: .seealso: `DMLabel`, `DM`, `DMLabel`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2226: @*/
2227: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2228: {
2229:   PetscFunctionBegin;
2230:   PetscCall(DMInitializePackage());
2231:   PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2232:   PetscFunctionReturn(PETSC_SUCCESS);
2233: }

2235: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2236: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);

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

2241:   Not Collective

2243:   Level: advanced

2245: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2246: @*/
2247: PetscErrorCode DMLabelRegisterAll(void)
2248: {
2249:   PetscFunctionBegin;
2250:   if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2251:   DMLabelRegisterAllCalled = PETSC_TRUE;

2253:   PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2254:   PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2255:   PetscFunctionReturn(PETSC_SUCCESS);
2256: }

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

2261:   Level: developer

2263: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2264: @*/
2265: PetscErrorCode DMLabelRegisterDestroy(void)
2266: {
2267:   PetscFunctionBegin;
2268:   PetscCall(PetscFunctionListDestroy(&DMLabelList));
2269:   DMLabelRegisterAllCalled = PETSC_FALSE;
2270:   PetscFunctionReturn(PETSC_SUCCESS);
2271: }

2273: /*@C
2274:   DMLabelSetType - Sets the particular implementation for a label.

2276:   Collective

2278:   Input Parameters:
2279: + label  - The label
2280: - method - The name of the label type

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

2285:   Level: intermediate

2287: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2288: @*/
2289: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2290: {
2291:   PetscErrorCode (*r)(DMLabel);
2292:   PetscBool match;

2294:   PetscFunctionBegin;
2296:   PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2297:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

2299:   PetscCall(DMLabelRegisterAll());
2300:   PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2301:   PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);

2303:   PetscTryTypeMethod(label, destroy);
2304:   PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2305:   PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2306:   PetscCall((*r)(label));
2307:   PetscFunctionReturn(PETSC_SUCCESS);
2308: }

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

2313:   Not Collective

2315:   Input Parameter:
2316: . label  - The `DMLabel`

2318:   Output Parameter:
2319: . type - The `DMLabel` type name

2321:   Level: intermediate

2323: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2324: @*/
2325: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2326: {
2327:   PetscFunctionBegin;
2330:   PetscCall(DMLabelRegisterAll());
2331:   *type = ((PetscObject)label)->type_name;
2332:   PetscFunctionReturn(PETSC_SUCCESS);
2333: }

2335: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2336: {
2337:   PetscFunctionBegin;
2338:   label->ops->view         = DMLabelView_Concrete;
2339:   label->ops->setup        = NULL;
2340:   label->ops->duplicate    = DMLabelDuplicate_Concrete;
2341:   label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2342:   PetscFunctionReturn(PETSC_SUCCESS);
2343: }

2345: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2346: {
2347:   PetscFunctionBegin;
2349:   PetscCall(DMLabelInitialize_Concrete(label));
2350:   PetscFunctionReturn(PETSC_SUCCESS);
2351: }

2353: /*@
2354:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2355:   the local section and an `PetscSF` describing the section point overlap.

2357:   Collective

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

2366:   Output Parameter:
2367: . gsection - The `PetscSection` for the global field layout

2369:   Level: developer

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

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

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

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

2447: typedef struct _n_PetscSectionSym_Label {
2448:   DMLabel              label;
2449:   PetscCopyMode       *modes;
2450:   PetscInt            *sizes;
2451:   const PetscInt    ***perms;
2452:   const PetscScalar ***rots;
2453:   PetscInt (*minMaxOrients)[2];
2454:   PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2455: } PetscSectionSym_Label;

2457: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2458: {
2459:   PetscInt               i, j;
2460:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;

2462:   PetscFunctionBegin;
2463:   for (i = 0; i <= sl->numStrata; i++) {
2464:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2465:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2466:         if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2467:         if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2468:       }
2469:       if (sl->perms[i]) {
2470:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

2472:         PetscCall(PetscFree(perms));
2473:       }
2474:       if (sl->rots[i]) {
2475:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

2477:         PetscCall(PetscFree(rots));
2478:       }
2479:     }
2480:   }
2481:   PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2482:   PetscCall(DMLabelDestroy(&sl->label));
2483:   sl->numStrata = 0;
2484:   PetscFunctionReturn(PETSC_SUCCESS);
2485: }

2487: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2488: {
2489:   PetscFunctionBegin;
2490:   PetscCall(PetscSectionSymLabelReset(sym));
2491:   PetscCall(PetscFree(sym->data));
2492:   PetscFunctionReturn(PETSC_SUCCESS);
2493: }

2495: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2496: {
2497:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2498:   PetscBool              isAscii;
2499:   DMLabel                label = sl->label;
2500:   const char            *name;

2502:   PetscFunctionBegin;
2503:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2504:   if (isAscii) {
2505:     PetscInt          i, j, k;
2506:     PetscViewerFormat format;

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

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

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

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

2577:   Logically

2579:   Input parameters:
2580: + sym - the section symmetries
2581: - label - the `DMLabel` describing the types of points

2583:   Level: developer:

2585: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2586: @*/
2587: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2588: {
2589:   PetscSectionSym_Label *sl;

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

2609: /*@C
2610:   PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum

2612:   Logically Collective

2614:   Input Parameters:
2615: + sym       - the section symmetries
2616: - stratum   - the stratum value in the label that we are assigning symmetries for

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

2625:   Level: developer

2627: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2628: @*/
2629: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2630: {
2631:   PetscSectionSym_Label *sl;
2632:   const char            *name;
2633:   PetscInt               i;

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

2642:     if (stratum == value) break;
2643:   }
2644:   PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2645:   PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2646:   if (size) {
2648:     *size = sl->sizes[i];
2649:   }
2650:   if (minOrient) {
2652:     *minOrient = sl->minMaxOrients[i][0];
2653:   }
2654:   if (maxOrient) {
2656:     *maxOrient = sl->minMaxOrients[i][1];
2657:   }
2658:   if (perms) {
2660:     *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;
2661:   }
2662:   if (rots) {
2664:     *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;
2665:   }
2666:   PetscFunctionReturn(PETSC_SUCCESS);
2667: }

2669: /*@C
2670:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

2672:   Logically

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

2684:   Level: developer

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

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

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

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

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

2741: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2742: {
2743:   PetscInt               i, j, numStrata;
2744:   PetscSectionSym_Label *sl;
2745:   DMLabel                label;

2747:   PetscFunctionBegin;
2748:   sl        = (PetscSectionSym_Label *)sym->data;
2749:   numStrata = sl->numStrata;
2750:   label     = sl->label;
2751:   for (i = 0; i < numPoints; i++) {
2752:     PetscInt point = points[2 * i];
2753:     PetscInt ornt  = points[2 * i + 1];

2755:     for (j = 0; j < numStrata; j++) {
2756:       if (label->validIS[j]) {
2757:         PetscInt k;

2759:         PetscCall(ISLocate(label->points[j], point, &k));
2760:         if (k >= 0) break;
2761:       } else {
2762:         PetscBool has;

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

2776: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2777: {
2778:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2779:   IS                     valIS;
2780:   const PetscInt        *values;
2781:   PetscInt               Nv, v;

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

2793:     PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2794:     PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
2795:   }
2796:   PetscCall(ISDestroy(&valIS));
2797:   PetscFunctionReturn(PETSC_SUCCESS);
2798: }

2800: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2801: {
2802:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2803:   DMLabel                dlabel;

2805:   PetscFunctionBegin;
2806:   PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
2807:   PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
2808:   PetscCall(DMLabelDestroy(&dlabel));
2809:   PetscCall(PetscSectionSymCopy(sym, *dsym));
2810:   PetscFunctionReturn(PETSC_SUCCESS);
2811: }

2813: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2814: {
2815:   PetscSectionSym_Label *sl;

2817:   PetscFunctionBegin;
2818:   PetscCall(PetscNew(&sl));
2819:   sym->ops->getpoints  = PetscSectionSymGetPoints_Label;
2820:   sym->ops->distribute = PetscSectionSymDistribute_Label;
2821:   sym->ops->copy       = PetscSectionSymCopy_Label;
2822:   sym->ops->view       = PetscSectionSymView_Label;
2823:   sym->ops->destroy    = PetscSectionSymDestroy_Label;
2824:   sym->data            = (void *)sl;
2825:   PetscFunctionReturn(PETSC_SUCCESS);
2826: }

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

2831:   Collective

2833:   Input Parameters:
2834: + comm - the MPI communicator for the new symmetry
2835: - label - the label defining the strata

2837:   Output Parameter:
2838: . sym - the section symmetries

2840:   Level: developer

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