Actual source code: dmlabel.c
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/sectionimpl.h>
4: #include <petscsf.h>
5: #include <petscsection.h>
7: PetscFunctionList DMLabelList = NULL;
8: PetscBool DMLabelRegisterAllCalled = PETSC_FALSE;
10: /*@
11: DMLabelCreate - Create a `DMLabel` object, which is a multimap
13: Collective
15: Input Parameters:
16: + comm - The communicator, usually `PETSC_COMM_SELF`
17: - name - The label name
19: Output Parameter:
20: . label - The `DMLabel`
22: Level: beginner
24: Notes:
25: The label name is actually usually the `PetscObject` name.
26: One can get/set it with `PetscObjectGetName()`/`PetscObjectSetName()`.
28: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`
29: @*/
30: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
31: {
32: PetscFunctionBegin;
33: PetscAssertPointer(label, 3);
34: PetscCall(DMInitializePackage());
36: PetscCall(PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView));
37: (*label)->numStrata = 0;
38: (*label)->defaultValue = -1;
39: (*label)->stratumValues = NULL;
40: (*label)->validIS = NULL;
41: (*label)->stratumSizes = NULL;
42: (*label)->points = NULL;
43: (*label)->ht = NULL;
44: (*label)->pStart = -1;
45: (*label)->pEnd = -1;
46: (*label)->bt = NULL;
47: PetscCall(PetscHMapICreate(&(*label)->hmap));
48: PetscCall(PetscObjectSetName((PetscObject)*label, name));
49: PetscCall(DMLabelSetType(*label, DMLABELCONCRETE));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@
54: DMLabelSetUp - SetUp a `DMLabel` object
56: Collective
58: Input Parameters:
59: . label - The `DMLabel`
61: Level: intermediate
63: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
64: @*/
65: PetscErrorCode DMLabelSetUp(DMLabel label)
66: {
67: PetscFunctionBegin;
69: PetscTryTypeMethod(label, setup);
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*
74: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
76: Not collective
78: Input parameter:
79: + label - The `DMLabel`
80: - v - The stratum value
82: Output parameter:
83: . label - The `DMLabel` with stratum in sorted list format
85: Level: developer
87: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
88: */
89: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
90: {
91: IS is;
92: PetscInt off = 0, *pointArray, p;
94: PetscFunctionBegin;
95: if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
96: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
97: PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
98: PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
99: PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
100: PetscCall(PetscHSetIClear(label->ht[v]));
101: PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
102: if (label->bt) {
103: for (p = 0; p < label->stratumSizes[v]; ++p) {
104: const PetscInt point = pointArray[p];
105: PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
106: PetscCall(PetscBTSet(label->bt, point - label->pStart));
107: }
108: }
109: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
110: PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
111: PetscCall(PetscFree(pointArray));
112: } else {
113: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
114: }
115: PetscCall(PetscObjectSetName((PetscObject)is, "indices"));
116: label->points[v] = is;
117: label->validIS[v] = PETSC_TRUE;
118: PetscCall(PetscObjectStateIncrease((PetscObject)label));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*
123: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
125: Not Collective
127: Input parameter:
128: . label - The `DMLabel`
130: Output parameter:
131: . label - The `DMLabel` with all strata in sorted list format
133: Level: developer
135: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
136: */
137: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
138: {
139: PetscInt v;
141: PetscFunctionBegin;
142: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeValid_Private(label, v));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*
147: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
149: Not Collective
151: Input parameter:
152: + label - The `DMLabel`
153: - v - The stratum value
155: Output parameter:
156: . label - The `DMLabel` with stratum in hash format
158: Level: developer
160: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
161: */
162: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
163: {
164: PetscInt p;
165: const PetscInt *points;
167: PetscFunctionBegin;
168: if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
169: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeInvalid_Private", v);
170: if (label->points[v]) {
171: PetscCall(ISGetIndices(label->points[v], &points));
172: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
173: PetscCall(ISRestoreIndices(label->points[v], &points));
174: PetscCall(ISDestroy(&label->points[v]));
175: }
176: label->validIS[v] = PETSC_FALSE;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
181: {
182: PetscInt v;
184: PetscFunctionBegin;
185: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeInvalid_Private(label, v));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
190: #define DMLABEL_LOOKUP_THRESHOLD 16
191: #endif
193: PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
194: {
195: PetscInt v;
197: PetscFunctionBegin;
198: *index = -1;
199: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD || label->readonly) {
200: for (v = 0; v < label->numStrata; ++v)
201: if (label->stratumValues[v] == value) {
202: *index = v;
203: break;
204: }
205: } else {
206: PetscCall(PetscHMapIGet(label->hmap, value, index));
207: }
208: if (PetscDefined(USE_DEBUG) && !label->readonly) { /* Check strata hash map consistency */
209: PetscInt len, loc = -1;
210: PetscCall(PetscHMapIGetSize(label->hmap, &len));
211: PetscCheck(len == label->numStrata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
212: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
213: PetscCall(PetscHMapIGet(label->hmap, value, &loc));
214: } else {
215: for (v = 0; v < label->numStrata; ++v)
216: if (label->stratumValues[v] == value) {
217: loc = v;
218: break;
219: }
220: }
221: PetscCheck(loc == *index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
222: }
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
227: {
228: PetscInt v;
229: PetscInt *tmpV;
230: PetscInt *tmpS;
231: PetscHSetI *tmpH, ht;
232: IS *tmpP, is;
233: PetscBool *tmpB;
234: PetscHMapI hmap = label->hmap;
236: PetscFunctionBegin;
237: v = label->numStrata;
238: tmpV = label->stratumValues;
239: tmpS = label->stratumSizes;
240: tmpH = label->ht;
241: tmpP = label->points;
242: tmpB = label->validIS;
243: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
244: PetscInt *oldV = tmpV;
245: PetscInt *oldS = tmpS;
246: PetscHSetI *oldH = tmpH;
247: IS *oldP = tmpP;
248: PetscBool *oldB = tmpB;
249: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV));
250: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS));
251: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpH), &tmpH));
252: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpP), &tmpP));
253: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB));
254: PetscCall(PetscArraycpy(tmpV, oldV, v));
255: PetscCall(PetscArraycpy(tmpS, oldS, v));
256: PetscCall(PetscArraycpy(tmpH, oldH, v));
257: PetscCall(PetscArraycpy(tmpP, oldP, v));
258: PetscCall(PetscArraycpy(tmpB, oldB, v));
259: PetscCall(PetscFree(oldV));
260: PetscCall(PetscFree(oldS));
261: PetscCall(PetscFree(oldH));
262: PetscCall(PetscFree(oldP));
263: PetscCall(PetscFree(oldB));
264: }
265: label->numStrata = v + 1;
266: label->stratumValues = tmpV;
267: label->stratumSizes = tmpS;
268: label->ht = tmpH;
269: label->points = tmpP;
270: label->validIS = tmpB;
271: PetscCall(PetscHSetICreate(&ht));
272: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
273: PetscCall(PetscHMapISet(hmap, value, v));
274: tmpV[v] = value;
275: tmpS[v] = 0;
276: tmpH[v] = ht;
277: tmpP[v] = is;
278: tmpB[v] = PETSC_TRUE;
279: PetscCall(PetscObjectStateIncrease((PetscObject)label));
280: *index = v;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
285: {
286: PetscFunctionBegin;
287: PetscCall(DMLabelLookupStratum(label, value, index));
288: if (*index < 0) PetscCall(DMLabelNewStratum(label, value, index));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
293: {
294: PetscFunctionBegin;
295: *size = 0;
296: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
297: if (label->readonly || label->validIS[v]) {
298: *size = label->stratumSizes[v];
299: } else {
300: PetscCall(PetscHSetIGetSize(label->ht[v], size));
301: }
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: DMLabelAddStratum - Adds a new stratum value in a `DMLabel`
308: Input Parameters:
309: + label - The `DMLabel`
310: - value - The stratum value
312: Level: beginner
314: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
315: @*/
316: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
317: {
318: PetscInt v;
320: PetscFunctionBegin;
322: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
323: PetscCall(DMLabelLookupAddStratum(label, value, &v));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: DMLabelAddStrata - Adds new stratum values in a `DMLabel`
330: Not Collective
332: Input Parameters:
333: + label - The `DMLabel`
334: . numStrata - The number of stratum values
335: - stratumValues - The stratum values
337: Level: beginner
339: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
340: @*/
341: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
342: {
343: PetscInt *values, v;
345: PetscFunctionBegin;
347: if (numStrata) PetscAssertPointer(stratumValues, 3);
348: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
349: PetscCall(PetscMalloc1(numStrata, &values));
350: PetscCall(PetscArraycpy(values, stratumValues, numStrata));
351: PetscCall(PetscSortRemoveDupsInt(&numStrata, values));
352: if (!label->numStrata) { /* Fast preallocation */
353: PetscInt *tmpV;
354: PetscInt *tmpS;
355: PetscHSetI *tmpH, ht;
356: IS *tmpP, is;
357: PetscBool *tmpB;
358: PetscHMapI hmap = label->hmap;
360: PetscCall(PetscMalloc1(numStrata, &tmpV));
361: PetscCall(PetscMalloc1(numStrata, &tmpS));
362: PetscCall(PetscCalloc1(numStrata, &tmpH));
363: PetscCall(PetscCalloc1(numStrata, &tmpP));
364: PetscCall(PetscMalloc1(numStrata, &tmpB));
365: label->numStrata = numStrata;
366: label->stratumValues = tmpV;
367: label->stratumSizes = tmpS;
368: label->ht = tmpH;
369: label->points = tmpP;
370: label->validIS = tmpB;
371: for (v = 0; v < numStrata; ++v) {
372: PetscCall(PetscHSetICreate(&ht));
373: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
374: PetscCall(PetscHMapISet(hmap, values[v], v));
375: tmpV[v] = values[v];
376: tmpS[v] = 0;
377: tmpH[v] = ht;
378: tmpP[v] = is;
379: tmpB[v] = PETSC_TRUE;
380: }
381: PetscCall(PetscObjectStateIncrease((PetscObject)label));
382: } else {
383: for (v = 0; v < numStrata; ++v) PetscCall(DMLabelAddStratum(label, values[v]));
384: }
385: PetscCall(PetscFree(values));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*@
390: DMLabelAddStrataIS - Adds new stratum values in a `DMLabel`
392: Not Collective
394: Input Parameters:
395: + label - The `DMLabel`
396: - valueIS - Index set with stratum values
398: Level: beginner
400: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
401: @*/
402: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
403: {
404: PetscInt numStrata;
405: const PetscInt *stratumValues;
407: PetscFunctionBegin;
410: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
411: PetscCall(ISGetLocalSize(valueIS, &numStrata));
412: PetscCall(ISGetIndices(valueIS, &stratumValues));
413: PetscCall(DMLabelAddStrata(label, numStrata, stratumValues));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode DMLabelView_Concrete_Ascii(DMLabel label, PetscViewer viewer)
418: {
419: PetscInt v;
420: PetscMPIInt rank;
422: PetscFunctionBegin;
423: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
424: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
425: if (label) {
426: const char *name;
428: PetscCall(PetscObjectGetName((PetscObject)label, &name));
429: PetscCall(PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name));
430: if (label->bt) PetscCall(PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd));
431: for (v = 0; v < label->numStrata; ++v) {
432: const PetscInt value = label->stratumValues[v];
433: const PetscInt *points;
434: PetscInt p;
436: PetscCall(ISGetIndices(label->points[v], &points));
437: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
438: PetscCall(ISRestoreIndices(label->points[v], &points));
439: }
440: }
441: PetscCall(PetscViewerFlush(viewer));
442: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
447: {
448: PetscBool iascii;
450: PetscFunctionBegin;
451: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
452: if (iascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /*@
457: DMLabelView - View the label
459: Collective
461: Input Parameters:
462: + label - The `DMLabel`
463: - viewer - The `PetscViewer`
465: Level: intermediate
467: .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
468: @*/
469: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
470: {
471: PetscFunctionBegin;
473: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer));
475: PetscCall(DMLabelMakeAllValid_Private(label));
476: PetscUseTypeMethod(label, view, viewer);
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*@
481: DMLabelReset - Destroys internal data structures in a `DMLabel`
483: Not Collective
485: Input Parameter:
486: . label - The `DMLabel`
488: Level: beginner
490: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
491: @*/
492: PetscErrorCode DMLabelReset(DMLabel label)
493: {
494: PetscInt v;
496: PetscFunctionBegin;
498: for (v = 0; v < label->numStrata; ++v) {
499: if (label->ht[v]) PetscCall(PetscHSetIDestroy(&label->ht[v]));
500: PetscCall(ISDestroy(&label->points[v]));
501: }
502: label->numStrata = 0;
503: PetscCall(PetscFree(label->stratumValues));
504: PetscCall(PetscFree(label->stratumSizes));
505: PetscCall(PetscFree(label->ht));
506: PetscCall(PetscFree(label->points));
507: PetscCall(PetscFree(label->validIS));
508: PetscCall(PetscHMapIReset(label->hmap));
509: label->pStart = -1;
510: label->pEnd = -1;
511: PetscCall(PetscBTDestroy(&label->bt));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: /*@
516: DMLabelDestroy - Destroys a `DMLabel`
518: Collective
520: Input Parameter:
521: . label - The `DMLabel`
523: Level: beginner
525: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
526: @*/
527: PetscErrorCode DMLabelDestroy(DMLabel *label)
528: {
529: PetscFunctionBegin;
530: if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
532: if (--((PetscObject)*label)->refct > 0) {
533: *label = NULL;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
536: PetscCall(DMLabelReset(*label));
537: PetscCall(PetscHMapIDestroy(&(*label)->hmap));
538: PetscCall(PetscHeaderDestroy(label));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
543: {
544: PetscFunctionBegin;
545: for (PetscInt v = 0; v < label->numStrata; ++v) {
546: PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
547: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
548: (*labelnew)->points[v] = label->points[v];
549: }
550: PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
551: PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@
556: DMLabelDuplicate - Duplicates a `DMLabel`
558: Collective
560: Input Parameter:
561: . label - The `DMLabel`
563: Output Parameter:
564: . labelnew - new label
566: Level: intermediate
568: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
569: @*/
570: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
571: {
572: const char *name;
574: PetscFunctionBegin;
576: PetscCall(DMLabelMakeAllValid_Private(label));
577: PetscCall(PetscObjectGetName((PetscObject)label, &name));
578: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));
580: (*labelnew)->numStrata = label->numStrata;
581: (*labelnew)->defaultValue = label->defaultValue;
582: (*labelnew)->readonly = label->readonly;
583: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
584: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
585: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
586: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
587: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
588: for (PetscInt v = 0; v < label->numStrata; ++v) {
589: (*labelnew)->stratumValues[v] = label->stratumValues[v];
590: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
591: (*labelnew)->validIS[v] = PETSC_TRUE;
592: }
593: (*labelnew)->pStart = -1;
594: (*labelnew)->pEnd = -1;
595: (*labelnew)->bt = NULL;
596: PetscUseTypeMethod(label, duplicate, labelnew);
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*@C
601: DMLabelCompare - Compare two `DMLabel` objects
603: Collective; No Fortran Support
605: Input Parameters:
606: + comm - Comm over which to compare labels
607: . l0 - First `DMLabel`
608: - l1 - Second `DMLabel`
610: Output Parameters:
611: + equal - (Optional) Flag whether the two labels are equal
612: - message - (Optional) Message describing the difference
614: Level: intermediate
616: Notes:
617: The output flag equal is the same on all processes.
618: If it is passed as `NULL` and difference is found, an error is thrown on all processes.
619: Make sure to pass `NULL` on all processes.
621: The output message is set independently on each rank.
622: It is set to `NULL` if no difference was found on the current rank. It must be freed by user.
623: If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner.
624: Make sure to pass `NULL` on all processes.
626: For the comparison, we ignore the order of stratum values, and strata with no points.
628: The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel.
630: Developer Note:
631: Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()`
633: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
634: @*/
635: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
636: {
637: const char *name0, *name1;
638: char msg[PETSC_MAX_PATH_LEN] = "";
639: PetscBool eq;
640: PetscMPIInt rank;
642: PetscFunctionBegin;
645: if (equal) PetscAssertPointer(equal, 4);
646: if (message) PetscAssertPointer(message, 5);
647: PetscCallMPI(MPI_Comm_rank(comm, &rank));
648: PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
649: PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
650: {
651: PetscInt v0, v1;
653: PetscCall(DMLabelGetDefaultValue(l0, &v0));
654: PetscCall(DMLabelGetDefaultValue(l1, &v1));
655: eq = (PetscBool)(v0 == v1);
656: 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));
657: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
658: if (!eq) goto finish;
659: }
660: {
661: IS is0, is1;
663: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
664: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
665: PetscCall(ISEqual(is0, is1, &eq));
666: PetscCall(ISDestroy(&is0));
667: PetscCall(ISDestroy(&is1));
668: if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
669: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
670: if (!eq) goto finish;
671: }
672: {
673: PetscInt i, nValues;
675: PetscCall(DMLabelGetNumValues(l0, &nValues));
676: for (i = 0; i < nValues; i++) {
677: const PetscInt v = l0->stratumValues[i];
678: PetscInt n;
679: IS is0, is1;
681: PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
682: if (!n) continue;
683: PetscCall(DMLabelGetStratumIS(l0, v, &is0));
684: PetscCall(DMLabelGetStratumIS(l1, v, &is1));
685: PetscCall(ISEqualUnsorted(is0, is1, &eq));
686: PetscCall(ISDestroy(&is0));
687: PetscCall(ISDestroy(&is1));
688: if (!eq) {
689: 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));
690: break;
691: }
692: }
693: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
694: }
695: finish:
696: /* If message output arg not set, print to stderr */
697: if (message) {
698: *message = NULL;
699: if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
700: } else {
701: if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
702: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
703: }
704: /* If same output arg not ser and labels are not equal, throw error */
705: if (equal) *equal = eq;
706: else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
713: Not Collective
715: Input Parameter:
716: . label - The `DMLabel`
718: Level: intermediate
720: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
721: @*/
722: PetscErrorCode DMLabelComputeIndex(DMLabel label)
723: {
724: PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v;
726: PetscFunctionBegin;
728: PetscCall(DMLabelMakeAllValid_Private(label));
729: for (v = 0; v < label->numStrata; ++v) {
730: const PetscInt *points;
731: PetscInt i;
733: PetscCall(ISGetIndices(label->points[v], &points));
734: for (i = 0; i < label->stratumSizes[v]; ++i) {
735: const PetscInt point = points[i];
737: pStart = PetscMin(point, pStart);
738: pEnd = PetscMax(point + 1, pEnd);
739: }
740: PetscCall(ISRestoreIndices(label->points[v], &points));
741: }
742: label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
743: label->pEnd = pEnd;
744: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: /*@
749: DMLabelCreateIndex - Create an index structure for membership determination
751: Not Collective
753: Input Parameters:
754: + label - The `DMLabel`
755: . pStart - The smallest point
756: - pEnd - The largest point + 1
758: Level: intermediate
760: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
761: @*/
762: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
763: {
764: PetscInt v;
766: PetscFunctionBegin;
768: PetscCall(DMLabelDestroyIndex(label));
769: PetscCall(DMLabelMakeAllValid_Private(label));
770: label->pStart = pStart;
771: label->pEnd = pEnd;
772: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
773: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
774: for (v = 0; v < label->numStrata; ++v) {
775: IS pointIS;
776: const PetscInt *points;
777: PetscInt i;
779: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
780: PetscCall(ISGetIndices(pointIS, &points));
781: for (i = 0; i < label->stratumSizes[v]; ++i) {
782: const PetscInt point = points[i];
784: 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);
785: PetscCall(PetscBTSet(label->bt, point - pStart));
786: }
787: PetscCall(ISRestoreIndices(label->points[v], &points));
788: PetscCall(ISDestroy(&pointIS));
789: }
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: DMLabelDestroyIndex - Destroy the index structure
796: Not Collective
798: Input Parameter:
799: . label - the `DMLabel`
801: Level: intermediate
803: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
804: @*/
805: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
806: {
807: PetscFunctionBegin;
809: label->pStart = -1;
810: label->pEnd = -1;
811: PetscCall(PetscBTDestroy(&label->bt));
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /*@
816: DMLabelGetBounds - Return the smallest and largest point in the label
818: Not Collective
820: Input Parameter:
821: . label - the `DMLabel`
823: Output Parameters:
824: + pStart - The smallest point
825: - pEnd - The largest point + 1
827: Level: intermediate
829: Note:
830: This will compute an index for the label if one does not exist.
832: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
833: @*/
834: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
835: {
836: PetscFunctionBegin;
838: if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
839: if (pStart) {
840: PetscAssertPointer(pStart, 2);
841: *pStart = label->pStart;
842: }
843: if (pEnd) {
844: PetscAssertPointer(pEnd, 3);
845: *pEnd = label->pEnd;
846: }
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@
851: DMLabelHasValue - Determine whether a label assigns the value to any point
853: Not Collective
855: Input Parameters:
856: + label - the `DMLabel`
857: - value - the value
859: Output Parameter:
860: . contains - Flag indicating whether the label maps this value to any point
862: Level: developer
864: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
865: @*/
866: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
867: {
868: PetscInt v;
870: PetscFunctionBegin;
872: PetscAssertPointer(contains, 3);
873: PetscCall(DMLabelLookupStratum(label, value, &v));
874: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: /*@
879: DMLabelHasPoint - Determine whether a label assigns a value to a point
881: Not Collective
883: Input Parameters:
884: + label - the `DMLabel`
885: - point - the point
887: Output Parameter:
888: . contains - Flag indicating whether the label maps this point to a value
890: Level: developer
892: Note:
893: The user must call `DMLabelCreateIndex()` before this function.
895: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
896: @*/
897: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
898: {
899: PetscFunctionBeginHot;
901: PetscAssertPointer(contains, 3);
902: PetscCall(DMLabelMakeAllValid_Private(label));
903: if (PetscDefined(USE_DEBUG)) {
904: PetscCheck(label->bt, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
905: 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);
906: }
907: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: /*@
912: DMLabelStratumHasPoint - Return true if the stratum contains a point
914: Not Collective
916: Input Parameters:
917: + label - the `DMLabel`
918: . value - the stratum value
919: - point - the point
921: Output Parameter:
922: . contains - true if the stratum contains the point
924: Level: intermediate
926: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
927: @*/
928: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
929: {
930: PetscFunctionBeginHot;
932: PetscAssertPointer(contains, 4);
933: if (value == label->defaultValue) {
934: PetscInt pointVal;
936: PetscCall(DMLabelGetValue(label, point, &pointVal));
937: *contains = (PetscBool)(pointVal == value);
938: } else {
939: PetscInt v;
941: PetscCall(DMLabelLookupStratum(label, value, &v));
942: if (v >= 0) {
943: if (label->validIS[v] || label->readonly) {
944: IS is;
945: PetscInt i;
947: PetscUseTypeMethod(label, getstratumis, v, &is);
948: PetscCall(ISLocate(is, point, &i));
949: PetscCall(ISDestroy(&is));
950: *contains = (PetscBool)(i >= 0);
951: } else {
952: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
953: }
954: } else { // value is not present
955: *contains = PETSC_FALSE;
956: }
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: /*@
962: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
963: When a label is created, it is initialized to -1.
965: Not Collective
967: Input Parameter:
968: . label - a `DMLabel` object
970: Output Parameter:
971: . defaultValue - the default value
973: Level: beginner
975: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
976: @*/
977: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
978: {
979: PetscFunctionBegin;
981: *defaultValue = label->defaultValue;
982: PetscFunctionReturn(PETSC_SUCCESS);
983: }
985: /*@
986: DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
987: When a label is created, it is initialized to -1.
989: Not Collective
991: Input Parameter:
992: . label - a `DMLabel` object
994: Output Parameter:
995: . defaultValue - the default value
997: Level: beginner
999: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1000: @*/
1001: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1002: {
1003: PetscFunctionBegin;
1005: label->defaultValue = defaultValue;
1006: PetscFunctionReturn(PETSC_SUCCESS);
1007: }
1009: /*@
1010: 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
1011: `DMLabelSetDefaultValue()`)
1013: Not Collective
1015: Input Parameters:
1016: + label - the `DMLabel`
1017: - point - the point
1019: Output Parameter:
1020: . value - The point value, or the default value (-1 by default)
1022: Level: intermediate
1024: Note:
1025: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1026: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1028: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1029: @*/
1030: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1031: {
1032: PetscInt v;
1034: PetscFunctionBeginHot;
1036: PetscAssertPointer(value, 3);
1037: *value = label->defaultValue;
1038: for (v = 0; v < label->numStrata; ++v) {
1039: if (label->validIS[v] || label->readonly) {
1040: IS is;
1041: PetscInt i;
1043: PetscUseTypeMethod(label, getstratumis, v, &is);
1044: PetscCall(ISLocate(label->points[v], point, &i));
1045: PetscCall(ISDestroy(&is));
1046: if (i >= 0) {
1047: *value = label->stratumValues[v];
1048: break;
1049: }
1050: } else {
1051: PetscBool has;
1053: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1054: if (has) {
1055: *value = label->stratumValues[v];
1056: break;
1057: }
1058: }
1059: }
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }
1063: /*@
1064: 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
1065: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1067: Not Collective
1069: Input Parameters:
1070: + label - the `DMLabel`
1071: . point - the point
1072: - value - The point value
1074: Level: intermediate
1076: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1077: @*/
1078: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1079: {
1080: PetscInt v;
1082: PetscFunctionBegin;
1084: /* Find label value, add new entry if needed */
1085: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1086: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1087: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1088: /* Set key */
1089: PetscCall(DMLabelMakeInvalid_Private(label, v));
1090: PetscCall(PetscHSetIAdd(label->ht[v], point));
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: /*@
1095: DMLabelClearValue - Clear the value a label assigns to a point
1097: Not Collective
1099: Input Parameters:
1100: + label - the `DMLabel`
1101: . point - the point
1102: - value - The point value
1104: Level: intermediate
1106: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1107: @*/
1108: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1109: {
1110: PetscInt v;
1112: PetscFunctionBegin;
1114: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1115: /* Find label value */
1116: PetscCall(DMLabelLookupStratum(label, value, &v));
1117: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1119: if (label->bt) {
1120: 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);
1121: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1122: }
1124: /* Delete key */
1125: PetscCall(DMLabelMakeInvalid_Private(label, v));
1126: PetscCall(PetscHSetIDel(label->ht[v], point));
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1130: /*@
1131: DMLabelInsertIS - Set all points in the `IS` to a value
1133: Not Collective
1135: Input Parameters:
1136: + label - the `DMLabel`
1137: . is - the point `IS`
1138: - value - The point value
1140: Level: intermediate
1142: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1143: @*/
1144: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1145: {
1146: PetscInt v, n, p;
1147: const PetscInt *points;
1149: PetscFunctionBegin;
1152: /* Find label value, add new entry if needed */
1153: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1154: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1155: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1156: /* Set keys */
1157: PetscCall(DMLabelMakeInvalid_Private(label, v));
1158: PetscCall(ISGetLocalSize(is, &n));
1159: PetscCall(ISGetIndices(is, &points));
1160: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1161: PetscCall(ISRestoreIndices(is, &points));
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: /*@
1166: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1168: Not Collective
1170: Input Parameter:
1171: . label - the `DMLabel`
1173: Output Parameter:
1174: . numValues - the number of values
1176: Level: intermediate
1178: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1179: @*/
1180: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1181: {
1182: PetscFunctionBegin;
1184: PetscAssertPointer(numValues, 2);
1185: *numValues = label->numStrata;
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: /*@
1190: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1192: Not Collective
1194: Input Parameter:
1195: . label - the `DMLabel`
1197: Output Parameter:
1198: . values - the value `IS`
1200: Level: intermediate
1202: Notes:
1203: The output `IS` should be destroyed when no longer needed.
1204: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1205: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1207: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1208: @*/
1209: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1210: {
1211: PetscFunctionBegin;
1213: PetscAssertPointer(values, 2);
1214: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1215: PetscFunctionReturn(PETSC_SUCCESS);
1216: }
1218: /*@
1219: DMLabelGetValueBounds - Return the smallest and largest value in the label
1221: Not Collective
1223: Input Parameter:
1224: . label - the `DMLabel`
1226: Output Parameters:
1227: + minValue - The smallest value
1228: - maxValue - The largest value
1230: Level: intermediate
1232: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1233: @*/
1234: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
1235: {
1236: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;
1238: PetscFunctionBegin;
1240: for (PetscInt v = 0; v < label->numStrata; ++v) {
1241: min = PetscMin(min, label->stratumValues[v]);
1242: max = PetscMax(max, label->stratumValues[v]);
1243: }
1244: if (minValue) {
1245: PetscAssertPointer(minValue, 2);
1246: *minValue = min;
1247: }
1248: if (maxValue) {
1249: PetscAssertPointer(maxValue, 3);
1250: *maxValue = max;
1251: }
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1258: Not Collective
1260: Input Parameter:
1261: . label - the `DMLabel`
1263: Output Parameter:
1264: . values - the value `IS`
1266: Level: intermediate
1268: Notes:
1269: The output `IS` should be destroyed when no longer needed.
1270: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1272: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1273: @*/
1274: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1275: {
1276: PetscInt i, j;
1277: PetscInt *valuesArr;
1279: PetscFunctionBegin;
1281: PetscAssertPointer(values, 2);
1282: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1283: for (i = 0, j = 0; i < label->numStrata; i++) {
1284: PetscInt n;
1286: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1287: if (n) valuesArr[j++] = label->stratumValues[i];
1288: }
1289: if (j == label->numStrata) {
1290: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1291: } else {
1292: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1293: }
1294: PetscCall(PetscFree(valuesArr));
1295: PetscFunctionReturn(PETSC_SUCCESS);
1296: }
1298: /*@
1299: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1301: Not Collective
1303: Input Parameters:
1304: + label - the `DMLabel`
1305: - value - the value
1307: Output Parameter:
1308: . index - the index of value in the list of values
1310: Level: intermediate
1312: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1313: @*/
1314: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1315: {
1316: PetscInt v;
1318: PetscFunctionBegin;
1320: PetscAssertPointer(index, 3);
1321: /* Do not assume they are sorted */
1322: for (v = 0; v < label->numStrata; ++v)
1323: if (label->stratumValues[v] == value) break;
1324: if (v >= label->numStrata) *index = -1;
1325: else *index = v;
1326: PetscFunctionReturn(PETSC_SUCCESS);
1327: }
1329: /*@
1330: DMLabelHasStratum - Determine whether points exist with the given value
1332: Not Collective
1334: Input Parameters:
1335: + label - the `DMLabel`
1336: - value - the stratum value
1338: Output Parameter:
1339: . exists - Flag saying whether points exist
1341: Level: intermediate
1343: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1344: @*/
1345: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1346: {
1347: PetscInt v;
1349: PetscFunctionBegin;
1351: PetscAssertPointer(exists, 3);
1352: PetscCall(DMLabelLookupStratum(label, value, &v));
1353: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /*@
1358: DMLabelGetStratumSize - Get the size of a stratum
1360: Not Collective
1362: Input Parameters:
1363: + label - the `DMLabel`
1364: - value - the stratum value
1366: Output Parameter:
1367: . size - The number of points in the stratum
1369: Level: intermediate
1371: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1372: @*/
1373: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1374: {
1375: PetscInt v;
1377: PetscFunctionBegin;
1379: PetscAssertPointer(size, 3);
1380: PetscCall(DMLabelLookupStratum(label, value, &v));
1381: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1382: PetscFunctionReturn(PETSC_SUCCESS);
1383: }
1385: /*@
1386: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1388: Not Collective
1390: Input Parameters:
1391: + label - the `DMLabel`
1392: - value - the stratum value
1394: Output Parameters:
1395: + start - the smallest point in the stratum
1396: - end - the largest point in the stratum
1398: Level: intermediate
1400: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1401: @*/
1402: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1403: {
1404: IS is;
1405: PetscInt v, min, max;
1407: PetscFunctionBegin;
1409: if (start) {
1410: PetscAssertPointer(start, 3);
1411: *start = -1;
1412: }
1413: if (end) {
1414: PetscAssertPointer(end, 4);
1415: *end = -1;
1416: }
1417: PetscCall(DMLabelLookupStratum(label, value, &v));
1418: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1419: PetscCall(DMLabelMakeValid_Private(label, v));
1420: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1421: PetscUseTypeMethod(label, getstratumis, v, &is);
1422: PetscCall(ISGetMinMax(is, &min, &max));
1423: PetscCall(ISDestroy(&is));
1424: if (start) *start = min;
1425: if (end) *end = max + 1;
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1429: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1430: {
1431: PetscFunctionBegin;
1432: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1433: *pointIS = label->points[v];
1434: PetscFunctionReturn(PETSC_SUCCESS);
1435: }
1437: /*@
1438: DMLabelGetStratumIS - Get an `IS` with the stratum points
1440: Not Collective
1442: Input Parameters:
1443: + label - the `DMLabel`
1444: - value - the stratum value
1446: Output Parameter:
1447: . points - The stratum points
1449: Level: intermediate
1451: Notes:
1452: The output `IS` should be destroyed when no longer needed.
1453: Returns `NULL` if the stratum is empty.
1455: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1456: @*/
1457: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1458: {
1459: PetscInt v;
1461: PetscFunctionBegin;
1463: PetscAssertPointer(points, 3);
1464: *points = NULL;
1465: PetscCall(DMLabelLookupStratum(label, value, &v));
1466: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1467: PetscCall(DMLabelMakeValid_Private(label, v));
1468: PetscUseTypeMethod(label, getstratumis, v, points);
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1472: /*@
1473: DMLabelSetStratumIS - Set the stratum points using an `IS`
1475: Not Collective
1477: Input Parameters:
1478: + label - the `DMLabel`
1479: . value - the stratum value
1480: - is - The stratum points
1482: Level: intermediate
1484: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1485: @*/
1486: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1487: {
1488: PetscInt v;
1490: PetscFunctionBegin;
1493: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1494: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1495: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1496: PetscCall(DMLabelClearStratum(label, value));
1497: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1498: PetscCall(PetscObjectReference((PetscObject)is));
1499: PetscCall(ISDestroy(&label->points[v]));
1500: label->points[v] = is;
1501: label->validIS[v] = PETSC_TRUE;
1502: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1503: if (label->bt) {
1504: const PetscInt *points;
1505: PetscInt p;
1507: PetscCall(ISGetIndices(is, &points));
1508: for (p = 0; p < label->stratumSizes[v]; ++p) {
1509: const PetscInt point = points[p];
1511: 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);
1512: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1513: }
1514: }
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: /*@
1519: DMLabelClearStratum - Remove a stratum
1521: Not Collective
1523: Input Parameters:
1524: + label - the `DMLabel`
1525: - value - the stratum value
1527: Level: intermediate
1529: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1530: @*/
1531: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1532: {
1533: PetscInt v;
1535: PetscFunctionBegin;
1537: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1538: PetscCall(DMLabelLookupStratum(label, value, &v));
1539: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1540: if (label->validIS[v]) {
1541: if (label->bt) {
1542: PetscInt i;
1543: const PetscInt *points;
1545: PetscCall(ISGetIndices(label->points[v], &points));
1546: for (i = 0; i < label->stratumSizes[v]; ++i) {
1547: const PetscInt point = points[i];
1549: 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);
1550: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1551: }
1552: PetscCall(ISRestoreIndices(label->points[v], &points));
1553: }
1554: label->stratumSizes[v] = 0;
1555: PetscCall(ISDestroy(&label->points[v]));
1556: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1557: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1558: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1559: } else {
1560: PetscCall(PetscHSetIClear(label->ht[v]));
1561: }
1562: PetscFunctionReturn(PETSC_SUCCESS);
1563: }
1565: /*@
1566: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1568: Not Collective
1570: Input Parameters:
1571: + label - The `DMLabel`
1572: . value - The label value for all points
1573: . pStart - The first point
1574: - pEnd - A point beyond all marked points
1576: Level: intermediate
1578: Note:
1579: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1581: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1582: @*/
1583: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1584: {
1585: IS pIS;
1587: PetscFunctionBegin;
1588: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1589: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1590: PetscCall(ISDestroy(&pIS));
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: /*@
1595: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1597: Not Collective
1599: Input Parameters:
1600: + label - The `DMLabel`
1601: . value - The label value
1602: - p - A point with this value
1604: Output Parameter:
1605: . 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
1607: Level: intermediate
1609: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1610: @*/
1611: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1612: {
1613: IS pointIS;
1614: const PetscInt *indices;
1615: PetscInt v;
1617: PetscFunctionBegin;
1619: PetscAssertPointer(index, 4);
1620: *index = -1;
1621: PetscCall(DMLabelLookupStratum(label, value, &v));
1622: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1623: PetscCall(DMLabelMakeValid_Private(label, v));
1624: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1625: PetscCall(ISGetIndices(pointIS, &indices));
1626: PetscCall(PetscFindInt(p, label->stratumSizes[v], indices, index));
1627: PetscCall(ISRestoreIndices(pointIS, &indices));
1628: PetscCall(ISDestroy(&pointIS));
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: /*@
1633: DMLabelFilter - Remove all points outside of [`start`, `end`)
1635: Not Collective
1637: Input Parameters:
1638: + label - the `DMLabel`
1639: . start - the first point kept
1640: - end - one more than the last point kept
1642: Level: intermediate
1644: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1645: @*/
1646: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1647: {
1648: PetscInt v;
1650: PetscFunctionBegin;
1652: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1653: PetscCall(DMLabelDestroyIndex(label));
1654: PetscCall(DMLabelMakeAllValid_Private(label));
1655: for (v = 0; v < label->numStrata; ++v) {
1656: PetscCall(ISGeneralFilter(label->points[v], start, end));
1657: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1658: }
1659: PetscCall(DMLabelCreateIndex(label, start, end));
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: /*@
1664: DMLabelPermute - Create a new label with permuted points
1666: Not Collective
1668: Input Parameters:
1669: + label - the `DMLabel`
1670: - permutation - the point permutation
1672: Output Parameter:
1673: . labelNew - the new label containing the permuted points
1675: Level: intermediate
1677: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1678: @*/
1679: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1680: {
1681: const PetscInt *perm;
1682: PetscInt numValues, numPoints, v, q;
1684: PetscFunctionBegin;
1687: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1688: PetscCall(DMLabelMakeAllValid_Private(label));
1689: PetscCall(DMLabelDuplicate(label, labelNew));
1690: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1691: PetscCall(ISGetLocalSize(permutation, &numPoints));
1692: PetscCall(ISGetIndices(permutation, &perm));
1693: for (v = 0; v < numValues; ++v) {
1694: const PetscInt size = (*labelNew)->stratumSizes[v];
1695: const PetscInt *points;
1696: PetscInt *pointsNew;
1698: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1699: PetscCall(PetscCalloc1(size, &pointsNew));
1700: for (q = 0; q < size; ++q) {
1701: const PetscInt point = points[q];
1703: 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);
1704: pointsNew[q] = perm[point];
1705: }
1706: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1707: PetscCall(PetscSortInt(size, pointsNew));
1708: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1709: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1710: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1711: PetscCall(PetscFree(pointsNew));
1712: } else {
1713: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1714: }
1715: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1716: }
1717: PetscCall(ISRestoreIndices(permutation, &perm));
1718: if (label->bt) {
1719: PetscCall(PetscBTDestroy(&label->bt));
1720: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1721: }
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }
1725: /*@
1726: DMLabelPermuteValues - Permute the values in a label
1728: Not collective
1730: Input Parameters:
1731: + label - the `DMLabel`
1732: - permutation - the value permutation, permutation[old value] = new value
1734: Output Parameter:
1735: . label - the `DMLabel` now with permuted values
1737: Note:
1738: The modification is done in-place
1740: Level: intermediate
1742: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1743: @*/
1744: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1745: {
1746: PetscInt Nv, Np;
1748: PetscFunctionBegin;
1751: PetscCall(DMLabelGetNumValues(label, &Nv));
1752: PetscCall(ISGetLocalSize(permutation, &Np));
1753: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1754: if (PetscDefined(USE_DEBUG)) {
1755: PetscBool flg;
1756: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1757: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1758: }
1759: PetscCall(DMLabelRewriteValues(label, permutation));
1760: PetscFunctionReturn(PETSC_SUCCESS);
1761: }
1763: /*@
1764: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1766: Not collective
1768: Input Parameters:
1769: + label - the `DMLabel`
1770: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1772: Output Parameter:
1773: . label - the `DMLabel` now with permuted values
1775: Note:
1776: The modification is done in-place
1778: Level: intermediate
1780: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1781: @*/
1782: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1783: {
1784: const PetscInt *perm;
1785: PetscInt Nv, Np;
1787: PetscFunctionBegin;
1790: PetscCall(DMLabelMakeAllValid_Private(label));
1791: PetscCall(DMLabelGetNumValues(label, &Nv));
1792: PetscCall(ISGetLocalSize(permutation, &Np));
1793: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1794: PetscCall(ISGetIndices(permutation, &perm));
1795: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1796: PetscCall(ISRestoreIndices(permutation, &perm));
1797: PetscFunctionReturn(PETSC_SUCCESS);
1798: }
1800: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1801: {
1802: MPI_Comm comm;
1803: PetscInt s, l, nroots, nleaves, offset, size;
1804: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1805: PetscSection rootSection;
1806: PetscSF labelSF;
1808: PetscFunctionBegin;
1809: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1810: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1811: /* Build a section of stratum values per point, generate the according SF
1812: and distribute point-wise stratum values to leaves. */
1813: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1814: PetscCall(PetscSectionCreate(comm, &rootSection));
1815: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1816: if (label) {
1817: for (s = 0; s < label->numStrata; ++s) {
1818: const PetscInt *points;
1820: PetscCall(ISGetIndices(label->points[s], &points));
1821: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1822: PetscCall(ISRestoreIndices(label->points[s], &points));
1823: }
1824: }
1825: PetscCall(PetscSectionSetUp(rootSection));
1826: /* Create a point-wise array of stratum values */
1827: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1828: PetscCall(PetscMalloc1(size, &rootStrata));
1829: PetscCall(PetscCalloc1(nroots, &rootIdx));
1830: if (label) {
1831: for (s = 0; s < label->numStrata; ++s) {
1832: const PetscInt *points;
1834: PetscCall(ISGetIndices(label->points[s], &points));
1835: for (l = 0; l < label->stratumSizes[s]; l++) {
1836: const PetscInt p = points[l];
1837: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1838: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1839: }
1840: PetscCall(ISRestoreIndices(label->points[s], &points));
1841: }
1842: }
1843: /* Build SF that maps label points to remote processes */
1844: PetscCall(PetscSectionCreate(comm, leafSection));
1845: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1846: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1847: PetscCall(PetscFree(remoteOffsets));
1848: /* Send the strata for each point over the derived SF */
1849: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1850: PetscCall(PetscMalloc1(size, leafStrata));
1851: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1852: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1853: /* Clean up */
1854: PetscCall(PetscFree(rootStrata));
1855: PetscCall(PetscFree(rootIdx));
1856: PetscCall(PetscSectionDestroy(&rootSection));
1857: PetscCall(PetscSFDestroy(&labelSF));
1858: PetscFunctionReturn(PETSC_SUCCESS);
1859: }
1861: /*@
1862: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1864: Collective
1866: Input Parameters:
1867: + label - the `DMLabel`
1868: - sf - the map from old to new distribution
1870: Output Parameter:
1871: . labelNew - the new redistributed label
1873: Level: intermediate
1875: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1876: @*/
1877: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1878: {
1879: MPI_Comm comm;
1880: PetscSection leafSection;
1881: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1882: PetscInt *leafStrata, *strataIdx;
1883: PetscInt **points;
1884: const char *lname = NULL;
1885: char *name;
1886: PetscMPIInt nameSize;
1887: PetscHSetI stratumHash;
1888: size_t len = 0;
1889: PetscMPIInt rank;
1891: PetscFunctionBegin;
1893: if (label) {
1895: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1896: PetscCall(DMLabelMakeAllValid_Private(label));
1897: }
1898: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1899: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1900: /* Bcast name */
1901: if (rank == 0) {
1902: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1903: PetscCall(PetscStrlen(lname, &len));
1904: }
1905: PetscCall(PetscMPIIntCast(len, &nameSize));
1906: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
1907: PetscCall(PetscMalloc1(nameSize + 1, &name));
1908: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1909: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1910: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1911: PetscCall(PetscFree(name));
1912: /* Bcast defaultValue */
1913: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1914: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
1915: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1916: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
1917: /* Determine received stratum values and initialise new label*/
1918: PetscCall(PetscHSetICreate(&stratumHash));
1919: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
1920: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
1921: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
1922: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
1923: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1924: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
1925: /* Turn leafStrata into indices rather than stratum values */
1926: offset = 0;
1927: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
1928: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
1929: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
1930: for (p = 0; p < size; ++p) {
1931: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1932: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1933: leafStrata[p] = s;
1934: break;
1935: }
1936: }
1937: }
1938: /* Rebuild the point strata on the receiver */
1939: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
1940: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1941: for (p = pStart; p < pEnd; p++) {
1942: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1943: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1944: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1945: }
1946: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
1947: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
1948: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
1949: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1950: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
1951: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
1952: }
1953: /* Insert points into new strata */
1954: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
1955: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1956: for (p = pStart; p < pEnd; p++) {
1957: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1958: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1959: for (s = 0; s < dof; s++) {
1960: stratum = leafStrata[offset + s];
1961: points[stratum][strataIdx[stratum]++] = p;
1962: }
1963: }
1964: for (s = 0; s < (*labelNew)->numStrata; s++) {
1965: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
1966: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
1967: }
1968: PetscCall(PetscFree(points));
1969: PetscCall(PetscHSetIDestroy(&stratumHash));
1970: PetscCall(PetscFree(leafStrata));
1971: PetscCall(PetscFree(strataIdx));
1972: PetscCall(PetscSectionDestroy(&leafSection));
1973: PetscFunctionReturn(PETSC_SUCCESS);
1974: }
1976: /*@
1977: DMLabelGather - Gather all label values from leafs into roots
1979: Collective
1981: Input Parameters:
1982: + label - the `DMLabel`
1983: - sf - the `PetscSF` communication map
1985: Output Parameter:
1986: . labelNew - the new `DMLabel` with localised leaf values
1988: Level: developer
1990: Note:
1991: This is the inverse operation to `DMLabelDistribute()`.
1993: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
1994: @*/
1995: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1996: {
1997: MPI_Comm comm;
1998: PetscSection rootSection;
1999: PetscSF sfLabel;
2000: PetscSFNode *rootPoints, *leafPoints;
2001: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2002: const PetscInt *rootDegree, *ilocal;
2003: PetscInt *rootStrata;
2004: const char *lname;
2005: char *name;
2006: PetscMPIInt nameSize;
2007: size_t len = 0;
2008: PetscMPIInt rank, size;
2010: PetscFunctionBegin;
2013: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2014: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2015: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2016: PetscCallMPI(MPI_Comm_size(comm, &size));
2017: /* Bcast name */
2018: if (rank == 0) {
2019: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2020: PetscCall(PetscStrlen(lname, &len));
2021: }
2022: PetscCall(PetscMPIIntCast(len, &nameSize));
2023: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2024: PetscCall(PetscMalloc1(nameSize + 1, &name));
2025: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2026: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2027: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2028: PetscCall(PetscFree(name));
2029: /* Gather rank/index pairs of leaves into local roots to build
2030: an inverse, multi-rooted SF. Note that this ignores local leaf
2031: indexing due to the use of the multiSF in PetscSFGather. */
2032: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2033: PetscCall(PetscMalloc1(nroots, &leafPoints));
2034: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2035: for (p = 0; p < nleaves; p++) {
2036: PetscInt ilp = ilocal ? ilocal[p] : p;
2038: leafPoints[ilp].index = ilp;
2039: leafPoints[ilp].rank = rank;
2040: }
2041: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2042: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2043: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2044: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2045: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2046: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2047: PetscCall(PetscSFCreate(comm, &sfLabel));
2048: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2049: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2050: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2051: /* Rebuild the point strata on the receiver */
2052: for (p = 0, idx = 0; p < nroots; p++) {
2053: for (d = 0; d < rootDegree[p]; d++) {
2054: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2055: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2056: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2057: }
2058: idx += rootDegree[p];
2059: }
2060: PetscCall(PetscFree(leafPoints));
2061: PetscCall(PetscFree(rootStrata));
2062: PetscCall(PetscSectionDestroy(&rootSection));
2063: PetscCall(PetscSFDestroy(&sfLabel));
2064: PetscFunctionReturn(PETSC_SUCCESS);
2065: }
2067: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2068: {
2069: const PetscInt *degree;
2070: const PetscInt *points;
2071: PetscInt Nr, r, Nl, l, val, defVal;
2073: PetscFunctionBegin;
2074: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2075: /* Add in leaves */
2076: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2077: for (l = 0; l < Nl; ++l) {
2078: PetscCall(DMLabelGetValue(label, points[l], &val));
2079: if (val != defVal) valArray[points[l]] = val;
2080: }
2081: /* Add in shared roots */
2082: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2083: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2084: for (r = 0; r < Nr; ++r) {
2085: if (degree[r]) {
2086: PetscCall(DMLabelGetValue(label, r, &val));
2087: if (val != defVal) valArray[r] = val;
2088: }
2089: }
2090: PetscFunctionReturn(PETSC_SUCCESS);
2091: }
2093: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
2094: {
2095: const PetscInt *degree;
2096: const PetscInt *points;
2097: PetscInt Nr, r, Nl, l, val, defVal;
2099: PetscFunctionBegin;
2100: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2101: /* Read out leaves */
2102: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2103: for (l = 0; l < Nl; ++l) {
2104: const PetscInt p = points[l];
2105: const PetscInt cval = valArray[p];
2107: if (cval != defVal) {
2108: PetscCall(DMLabelGetValue(label, p, &val));
2109: if (val == defVal) {
2110: PetscCall(DMLabelSetValue(label, p, cval));
2111: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2112: }
2113: }
2114: }
2115: /* Read out shared roots */
2116: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2117: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2118: for (r = 0; r < Nr; ++r) {
2119: if (degree[r]) {
2120: const PetscInt cval = valArray[r];
2122: if (cval != defVal) {
2123: PetscCall(DMLabelGetValue(label, r, &val));
2124: if (val == defVal) {
2125: PetscCall(DMLabelSetValue(label, r, cval));
2126: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2127: }
2128: }
2129: }
2130: }
2131: PetscFunctionReturn(PETSC_SUCCESS);
2132: }
2134: /*@
2135: DMLabelPropagateBegin - Setup a cycle of label propagation
2137: Collective
2139: Input Parameters:
2140: + label - The `DMLabel` to propagate across processes
2141: - sf - The `PetscSF` describing parallel layout of the label points
2143: Level: intermediate
2145: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2146: @*/
2147: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2148: {
2149: PetscInt Nr, r, defVal;
2150: PetscMPIInt size;
2152: PetscFunctionBegin;
2153: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2154: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2155: if (size > 1) {
2156: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2157: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2158: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2159: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2160: }
2161: PetscFunctionReturn(PETSC_SUCCESS);
2162: }
2164: /*@
2165: DMLabelPropagateEnd - Tear down a cycle of label propagation
2167: Collective
2169: Input Parameters:
2170: + label - The `DMLabel` to propagate across processes
2171: - pointSF - The `PetscSF` describing parallel layout of the label points
2173: Level: intermediate
2175: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2176: @*/
2177: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2178: {
2179: PetscFunctionBegin;
2180: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2181: PetscCall(PetscFree(label->propArray));
2182: label->propArray = NULL;
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: }
2186: /*@C
2187: DMLabelPropagatePush - Tear down a cycle of label propagation
2189: Collective
2191: Input Parameters:
2192: + label - The `DMLabel` to propagate across processes
2193: . pointSF - The `PetscSF` describing parallel layout of the label points
2194: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2195: - ctx - An optional user context for the callback, or `NULL`
2197: Calling sequence of `markPoint`:
2198: + label - The `DMLabel`
2199: . p - The point being marked
2200: . val - The label value for `p`
2201: - ctx - An optional user context
2203: Level: intermediate
2205: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2206: @*/
2207: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx)
2208: {
2209: PetscInt *valArray = label->propArray, Nr;
2210: PetscMPIInt size;
2212: PetscFunctionBegin;
2213: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2214: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2215: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2216: if (size > 1 && Nr >= 0) {
2217: /* Communicate marked edges
2218: The current implementation allocates an array the size of the number of root. We put the label values into the
2219: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2221: TODO: We could use in-place communication with a different SF
2222: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2223: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2225: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2226: 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
2227: edge to the queue.
2228: */
2229: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2230: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2231: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2232: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2233: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2234: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2235: }
2236: PetscFunctionReturn(PETSC_SUCCESS);
2237: }
2239: /*@
2240: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2242: Not Collective
2244: Input Parameter:
2245: . label - the `DMLabel`
2247: Output Parameters:
2248: + section - the section giving offsets for each stratum
2249: - is - An `IS` containing all the label points
2251: Level: developer
2253: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2254: @*/
2255: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2256: {
2257: IS vIS;
2258: const PetscInt *values;
2259: PetscInt *points;
2260: PetscInt nV, vS = 0, vE = 0, v, N;
2262: PetscFunctionBegin;
2264: PetscCall(DMLabelGetNumValues(label, &nV));
2265: PetscCall(DMLabelGetValueIS(label, &vIS));
2266: PetscCall(ISGetIndices(vIS, &values));
2267: if (nV) {
2268: vS = values[0];
2269: vE = values[0] + 1;
2270: }
2271: for (v = 1; v < nV; ++v) {
2272: vS = PetscMin(vS, values[v]);
2273: vE = PetscMax(vE, values[v] + 1);
2274: }
2275: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2276: PetscCall(PetscSectionSetChart(*section, vS, vE));
2277: for (v = 0; v < nV; ++v) {
2278: PetscInt n;
2280: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2281: PetscCall(PetscSectionSetDof(*section, values[v], n));
2282: }
2283: PetscCall(PetscSectionSetUp(*section));
2284: PetscCall(PetscSectionGetStorageSize(*section, &N));
2285: PetscCall(PetscMalloc1(N, &points));
2286: for (v = 0; v < nV; ++v) {
2287: IS is;
2288: const PetscInt *spoints;
2289: PetscInt dof, off, p;
2291: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2292: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2293: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2294: PetscCall(ISGetIndices(is, &spoints));
2295: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2296: PetscCall(ISRestoreIndices(is, &spoints));
2297: PetscCall(ISDestroy(&is));
2298: }
2299: PetscCall(ISRestoreIndices(vIS, &values));
2300: PetscCall(ISDestroy(&vIS));
2301: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2302: PetscFunctionReturn(PETSC_SUCCESS);
2303: }
2305: /*@C
2306: DMLabelRegister - Adds a new label component implementation
2308: Not Collective
2310: Input Parameters:
2311: + name - The name of a new user-defined creation routine
2312: - create_func - The creation routine itself
2314: Notes:
2315: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2317: Example Usage:
2318: .vb
2319: DMLabelRegister("my_label", MyLabelCreate);
2320: .ve
2322: Then, your label type can be chosen with the procedural interface via
2323: .vb
2324: DMLabelCreate(MPI_Comm, DMLabel *);
2325: DMLabelSetType(DMLabel, "my_label");
2326: .ve
2327: or at runtime via the option
2328: .vb
2329: -dm_label_type my_label
2330: .ve
2332: Level: advanced
2334: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2335: @*/
2336: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2337: {
2338: PetscFunctionBegin;
2339: PetscCall(DMInitializePackage());
2340: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2341: PetscFunctionReturn(PETSC_SUCCESS);
2342: }
2344: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2345: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2347: /*@C
2348: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2350: Not Collective
2352: Level: advanced
2354: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2355: @*/
2356: PetscErrorCode DMLabelRegisterAll(void)
2357: {
2358: PetscFunctionBegin;
2359: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2360: DMLabelRegisterAllCalled = PETSC_TRUE;
2362: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2363: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2364: PetscFunctionReturn(PETSC_SUCCESS);
2365: }
2367: /*@C
2368: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2370: Level: developer
2372: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2373: @*/
2374: PetscErrorCode DMLabelRegisterDestroy(void)
2375: {
2376: PetscFunctionBegin;
2377: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2378: DMLabelRegisterAllCalled = PETSC_FALSE;
2379: PetscFunctionReturn(PETSC_SUCCESS);
2380: }
2382: /*@
2383: DMLabelSetType - Sets the particular implementation for a label.
2385: Collective
2387: Input Parameters:
2388: + label - The label
2389: - method - The name of the label type
2391: Options Database Key:
2392: . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType`
2394: Level: intermediate
2396: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2397: @*/
2398: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2399: {
2400: PetscErrorCode (*r)(DMLabel);
2401: PetscBool match;
2403: PetscFunctionBegin;
2405: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2406: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2408: PetscCall(DMLabelRegisterAll());
2409: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2410: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2412: PetscTryTypeMethod(label, destroy);
2413: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2414: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2415: PetscCall((*r)(label));
2416: PetscFunctionReturn(PETSC_SUCCESS);
2417: }
2419: /*@
2420: DMLabelGetType - Gets the type name (as a string) from the label.
2422: Not Collective
2424: Input Parameter:
2425: . label - The `DMLabel`
2427: Output Parameter:
2428: . type - The `DMLabel` type name
2430: Level: intermediate
2432: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2433: @*/
2434: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2435: {
2436: PetscFunctionBegin;
2438: PetscAssertPointer(type, 2);
2439: PetscCall(DMLabelRegisterAll());
2440: *type = ((PetscObject)label)->type_name;
2441: PetscFunctionReturn(PETSC_SUCCESS);
2442: }
2444: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2445: {
2446: PetscFunctionBegin;
2447: label->ops->view = DMLabelView_Concrete;
2448: label->ops->setup = NULL;
2449: label->ops->duplicate = DMLabelDuplicate_Concrete;
2450: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2451: PetscFunctionReturn(PETSC_SUCCESS);
2452: }
2454: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2455: {
2456: PetscFunctionBegin;
2458: PetscCall(DMLabelInitialize_Concrete(label));
2459: PetscFunctionReturn(PETSC_SUCCESS);
2460: }
2462: /*@
2463: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2464: the local section and an `PetscSF` describing the section point overlap.
2466: Collective
2468: Input Parameters:
2469: + s - The `PetscSection` for the local field layout
2470: . sf - The `PetscSF` describing parallel layout of the section points
2471: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2472: . label - The label specifying the points
2473: - labelValue - The label stratum specifying the points
2475: Output Parameter:
2476: . gsection - The `PetscSection` for the global field layout
2478: Level: developer
2480: Note:
2481: This gives negative sizes and offsets to points not owned by this process
2483: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2484: @*/
2485: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2486: {
2487: PetscInt *neg = NULL, *tmpOff = NULL;
2488: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2490: PetscFunctionBegin;
2494: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2495: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2496: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2497: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2498: if (nroots >= 0) {
2499: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2500: PetscCall(PetscCalloc1(nroots, &neg));
2501: if (nroots > pEnd - pStart) {
2502: PetscCall(PetscCalloc1(nroots, &tmpOff));
2503: } else {
2504: tmpOff = &(*gsection)->atlasDof[-pStart];
2505: }
2506: }
2507: /* Mark ghost points with negative dof */
2508: for (p = pStart; p < pEnd; ++p) {
2509: PetscInt value;
2511: PetscCall(DMLabelGetValue(label, p, &value));
2512: if (value != labelValue) continue;
2513: PetscCall(PetscSectionGetDof(s, p, &dof));
2514: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2515: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2516: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2517: if (neg) neg[p] = -(dof + 1);
2518: }
2519: PetscCall(PetscSectionSetUpBC(*gsection));
2520: if (nroots >= 0) {
2521: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2522: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2523: if (nroots > pEnd - pStart) {
2524: for (p = pStart; p < pEnd; ++p) {
2525: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2526: }
2527: }
2528: }
2529: /* Calculate new sizes, get process offset, and calculate point offsets */
2530: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2531: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2532: (*gsection)->atlasOff[p] = off;
2533: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2534: }
2535: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2536: globalOff -= off;
2537: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2538: (*gsection)->atlasOff[p] += globalOff;
2539: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2540: }
2541: /* Put in negative offsets for ghost points */
2542: if (nroots >= 0) {
2543: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2544: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2545: if (nroots > pEnd - pStart) {
2546: for (p = pStart; p < pEnd; ++p) {
2547: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2548: }
2549: }
2550: }
2551: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2552: PetscCall(PetscFree(neg));
2553: PetscFunctionReturn(PETSC_SUCCESS);
2554: }
2556: typedef struct _n_PetscSectionSym_Label {
2557: DMLabel label;
2558: PetscCopyMode *modes;
2559: PetscInt *sizes;
2560: const PetscInt ***perms;
2561: const PetscScalar ***rots;
2562: PetscInt (*minMaxOrients)[2];
2563: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2564: } PetscSectionSym_Label;
2566: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2567: {
2568: PetscInt i, j;
2569: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2571: PetscFunctionBegin;
2572: for (i = 0; i <= sl->numStrata; i++) {
2573: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2574: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2575: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2576: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2577: }
2578: if (sl->perms[i]) {
2579: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2581: PetscCall(PetscFree(perms));
2582: }
2583: if (sl->rots[i]) {
2584: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2586: PetscCall(PetscFree(rots));
2587: }
2588: }
2589: }
2590: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2591: PetscCall(DMLabelDestroy(&sl->label));
2592: sl->numStrata = 0;
2593: PetscFunctionReturn(PETSC_SUCCESS);
2594: }
2596: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2597: {
2598: PetscFunctionBegin;
2599: PetscCall(PetscSectionSymLabelReset(sym));
2600: PetscCall(PetscFree(sym->data));
2601: PetscFunctionReturn(PETSC_SUCCESS);
2602: }
2604: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2605: {
2606: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2607: PetscBool isAscii;
2608: DMLabel label = sl->label;
2609: const char *name;
2611: PetscFunctionBegin;
2612: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2613: if (isAscii) {
2614: PetscInt i, j, k;
2615: PetscViewerFormat format;
2617: PetscCall(PetscViewerGetFormat(viewer, &format));
2618: if (label) {
2619: PetscCall(PetscViewerGetFormat(viewer, &format));
2620: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2621: PetscCall(PetscViewerASCIIPushTab(viewer));
2622: PetscCall(DMLabelView(label, viewer));
2623: PetscCall(PetscViewerASCIIPopTab(viewer));
2624: } else {
2625: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2626: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2627: }
2628: } else {
2629: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2630: }
2631: PetscCall(PetscViewerASCIIPushTab(viewer));
2632: for (i = 0; i <= sl->numStrata; i++) {
2633: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2635: if (!(sl->perms[i] || sl->rots[i])) {
2636: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2637: } else {
2638: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2639: PetscCall(PetscViewerASCIIPushTab(viewer));
2640: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2641: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2642: PetscCall(PetscViewerASCIIPushTab(viewer));
2643: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2644: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2645: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2646: } else {
2647: PetscInt tab;
2649: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2650: PetscCall(PetscViewerASCIIPushTab(viewer));
2651: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2652: if (sl->perms[i] && sl->perms[i][j]) {
2653: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2654: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2655: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2656: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2657: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2658: }
2659: if (sl->rots[i] && sl->rots[i][j]) {
2660: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2661: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2662: #if defined(PETSC_USE_COMPLEX)
2663: 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])));
2664: #else
2665: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2666: #endif
2667: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2668: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2669: }
2670: PetscCall(PetscViewerASCIIPopTab(viewer));
2671: }
2672: }
2673: PetscCall(PetscViewerASCIIPopTab(viewer));
2674: }
2675: PetscCall(PetscViewerASCIIPopTab(viewer));
2676: }
2677: }
2678: PetscCall(PetscViewerASCIIPopTab(viewer));
2679: }
2680: PetscFunctionReturn(PETSC_SUCCESS);
2681: }
2683: /*@
2684: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2686: Logically
2688: Input Parameters:
2689: + sym - the section symmetries
2690: - label - the `DMLabel` describing the types of points
2692: Level: developer:
2694: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2695: @*/
2696: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2697: {
2698: PetscSectionSym_Label *sl;
2700: PetscFunctionBegin;
2702: sl = (PetscSectionSym_Label *)sym->data;
2703: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2704: if (label) {
2705: sl->label = label;
2706: PetscCall(PetscObjectReference((PetscObject)label));
2707: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2708: 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));
2709: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2710: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2711: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2712: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2713: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2714: }
2715: PetscFunctionReturn(PETSC_SUCCESS);
2716: }
2718: /*@C
2719: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2721: Logically Collective
2723: Input Parameters:
2724: + sym - the section symmetries
2725: - stratum - the stratum value in the label that we are assigning symmetries for
2727: Output Parameters:
2728: + size - the number of dofs for points in the `stratum` of the label
2729: . minOrient - the smallest orientation for a point in this `stratum`
2730: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2731: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2732: - 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
2734: Level: developer
2736: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2737: @*/
2738: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2739: {
2740: PetscSectionSym_Label *sl;
2741: const char *name;
2742: PetscInt i;
2744: PetscFunctionBegin;
2746: sl = (PetscSectionSym_Label *)sym->data;
2747: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2748: for (i = 0; i <= sl->numStrata; i++) {
2749: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2751: if (stratum == value) break;
2752: }
2753: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2754: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2755: if (size) {
2756: PetscAssertPointer(size, 3);
2757: *size = sl->sizes[i];
2758: }
2759: if (minOrient) {
2760: PetscAssertPointer(minOrient, 4);
2761: *minOrient = sl->minMaxOrients[i][0];
2762: }
2763: if (maxOrient) {
2764: PetscAssertPointer(maxOrient, 5);
2765: *maxOrient = sl->minMaxOrients[i][1];
2766: }
2767: if (perms) {
2768: PetscAssertPointer(perms, 6);
2769: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2770: }
2771: if (rots) {
2772: PetscAssertPointer(rots, 7);
2773: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2774: }
2775: PetscFunctionReturn(PETSC_SUCCESS);
2776: }
2778: /*@C
2779: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2781: Logically
2783: Input Parameters:
2784: + sym - the section symmetries
2785: . stratum - the stratum value in the label that we are assigning symmetries for
2786: . size - the number of dofs for points in the `stratum` of the label
2787: . minOrient - the smallest orientation for a point in this `stratum`
2788: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2789: . mode - how `sym` should copy the `perms` and `rots` arrays
2790: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2791: - 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
2793: Level: developer
2795: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2796: @*/
2797: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2798: {
2799: PetscSectionSym_Label *sl;
2800: const char *name;
2801: PetscInt i, j, k;
2803: PetscFunctionBegin;
2805: sl = (PetscSectionSym_Label *)sym->data;
2806: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2807: for (i = 0; i <= sl->numStrata; i++) {
2808: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2810: if (stratum == value) break;
2811: }
2812: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2813: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2814: sl->sizes[i] = size;
2815: sl->modes[i] = mode;
2816: sl->minMaxOrients[i][0] = minOrient;
2817: sl->minMaxOrients[i][1] = maxOrient;
2818: if (mode == PETSC_COPY_VALUES) {
2819: if (perms) {
2820: PetscInt **ownPerms;
2822: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2823: for (j = 0; j < maxOrient - minOrient; j++) {
2824: if (perms[j]) {
2825: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2826: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2827: }
2828: }
2829: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2830: }
2831: if (rots) {
2832: PetscScalar **ownRots;
2834: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2835: for (j = 0; j < maxOrient - minOrient; j++) {
2836: if (rots[j]) {
2837: PetscCall(PetscMalloc1(size, &ownRots[j]));
2838: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2839: }
2840: }
2841: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2842: }
2843: } else {
2844: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2845: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2846: }
2847: PetscFunctionReturn(PETSC_SUCCESS);
2848: }
2850: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2851: {
2852: PetscInt i, j, numStrata;
2853: PetscSectionSym_Label *sl;
2854: DMLabel label;
2856: PetscFunctionBegin;
2857: sl = (PetscSectionSym_Label *)sym->data;
2858: numStrata = sl->numStrata;
2859: label = sl->label;
2860: for (i = 0; i < numPoints; i++) {
2861: PetscInt point = points[2 * i];
2862: PetscInt ornt = points[2 * i + 1];
2864: for (j = 0; j < numStrata; j++) {
2865: if (label->validIS[j]) {
2866: PetscInt k;
2868: PetscCall(ISLocate(label->points[j], point, &k));
2869: if (k >= 0) break;
2870: } else {
2871: PetscBool has;
2873: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2874: if (has) break;
2875: }
2876: }
2877: 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],
2878: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2879: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2880: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2881: }
2882: PetscFunctionReturn(PETSC_SUCCESS);
2883: }
2885: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2886: {
2887: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2888: IS valIS;
2889: const PetscInt *values;
2890: PetscInt Nv, v;
2892: PetscFunctionBegin;
2893: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2894: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2895: PetscCall(ISGetIndices(valIS, &values));
2896: for (v = 0; v < Nv; ++v) {
2897: const PetscInt val = values[v];
2898: PetscInt size, minOrient, maxOrient;
2899: const PetscInt **perms;
2900: const PetscScalar **rots;
2902: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2903: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
2904: }
2905: PetscCall(ISDestroy(&valIS));
2906: PetscFunctionReturn(PETSC_SUCCESS);
2907: }
2909: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2910: {
2911: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2912: DMLabel dlabel;
2914: PetscFunctionBegin;
2915: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
2916: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
2917: PetscCall(DMLabelDestroy(&dlabel));
2918: PetscCall(PetscSectionSymCopy(sym, *dsym));
2919: PetscFunctionReturn(PETSC_SUCCESS);
2920: }
2922: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2923: {
2924: PetscSectionSym_Label *sl;
2926: PetscFunctionBegin;
2927: PetscCall(PetscNew(&sl));
2928: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2929: sym->ops->distribute = PetscSectionSymDistribute_Label;
2930: sym->ops->copy = PetscSectionSymCopy_Label;
2931: sym->ops->view = PetscSectionSymView_Label;
2932: sym->ops->destroy = PetscSectionSymDestroy_Label;
2933: sym->data = (void *)sl;
2934: PetscFunctionReturn(PETSC_SUCCESS);
2935: }
2937: /*@
2938: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2940: Collective
2942: Input Parameters:
2943: + comm - the MPI communicator for the new symmetry
2944: - label - the label defining the strata
2946: Output Parameter:
2947: . sym - the section symmetries
2949: Level: developer
2951: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2952: @*/
2953: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2954: {
2955: PetscFunctionBegin;
2956: PetscCall(DMInitializePackage());
2957: PetscCall(PetscSectionSymCreate(comm, sym));
2958: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
2959: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
2960: PetscFunctionReturn(PETSC_SUCCESS);
2961: }