Actual source code: dmlabel.c
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/sectionimpl.h>
4: #include <petscsf.h>
5: #include <petscsection.h>
7: PetscFunctionList DMLabelList = NULL;
8: PetscBool DMLabelRegisterAllCalled = PETSC_FALSE;
10: /*@
11: DMLabelCreate - Create a `DMLabel` object, which is a multimap
13: Collective
15: Input Parameters:
16: + comm - The communicator, usually `PETSC_COMM_SELF`
17: - name - The label name
19: Output Parameter:
20: . label - The `DMLabel`
22: Level: beginner
24: Notes:
25: The label name is actually usually the `PetscObject` name.
26: One can get/set it with `PetscObjectGetName()`/`PetscObjectSetName()`.
28: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`
29: @*/
30: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
31: {
32: PetscFunctionBegin;
33: PetscAssertPointer(label, 3);
34: PetscCall(DMInitializePackage());
36: PetscCall(PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView));
37: (*label)->numStrata = 0;
38: (*label)->defaultValue = -1;
39: (*label)->stratumValues = NULL;
40: (*label)->validIS = NULL;
41: (*label)->stratumSizes = NULL;
42: (*label)->points = NULL;
43: (*label)->ht = NULL;
44: (*label)->pStart = -1;
45: (*label)->pEnd = -1;
46: (*label)->bt = NULL;
47: PetscCall(PetscHMapICreate(&(*label)->hmap));
48: PetscCall(PetscObjectSetName((PetscObject)*label, name));
49: PetscCall(DMLabelSetType(*label, DMLABELCONCRETE));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@
54: DMLabelSetUp - SetUp a `DMLabel` object
56: Collective
58: Input Parameters:
59: . label - The `DMLabel`
61: Level: intermediate
63: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
64: @*/
65: PetscErrorCode DMLabelSetUp(DMLabel label)
66: {
67: PetscFunctionBegin;
69: PetscTryTypeMethod(label, setup);
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*
74: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
76: Not collective
78: Input parameter:
79: + label - The `DMLabel`
80: - v - The stratum value
82: Output parameter:
83: . label - The `DMLabel` with stratum in sorted list format
85: Level: developer
87: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
88: */
89: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
90: {
91: IS is;
92: PetscInt off = 0, *pointArray, p;
94: PetscFunctionBegin;
95: if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
96: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
97: PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
98: PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
99: PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
100: PetscCall(PetscHSetIClear(label->ht[v]));
101: PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
102: if (label->bt) {
103: for (p = 0; p < label->stratumSizes[v]; ++p) {
104: const PetscInt point = pointArray[p];
105: PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
106: PetscCall(PetscBTSet(label->bt, point - label->pStart));
107: }
108: }
109: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
110: PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
111: PetscCall(PetscFree(pointArray));
112: } else {
113: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
114: }
115: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, PETSC_TRUE));
116: PetscCall(PetscObjectSetName((PetscObject)is, "indices"));
117: label->points[v] = is;
118: label->validIS[v] = PETSC_TRUE;
119: PetscCall(PetscObjectStateIncrease((PetscObject)label));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*
124: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
126: Not Collective
128: Input parameter:
129: . label - The `DMLabel`
131: Output parameter:
132: . label - The `DMLabel` with all strata in sorted list format
134: Level: developer
136: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
137: */
138: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
139: {
140: PetscInt v;
142: PetscFunctionBegin;
143: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeValid_Private(label, v));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*
148: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
150: Not Collective
152: Input parameter:
153: + label - The `DMLabel`
154: - v - The stratum value
156: Output parameter:
157: . label - The `DMLabel` with stratum in hash format
159: Level: developer
161: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
162: */
163: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
164: {
165: 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 (PetscInt 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;
435: PetscCall(ISGetIndices(label->points[v], &points));
436: for (PetscInt p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
437: PetscCall(ISRestoreIndices(label->points[v], &points));
438: }
439: }
440: PetscCall(PetscViewerFlush(viewer));
441: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
446: {
447: PetscBool isascii;
449: PetscFunctionBegin;
450: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
451: if (isascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: /*@
456: DMLabelView - View the label
458: Collective
460: Input Parameters:
461: + label - The `DMLabel`
462: - viewer - The `PetscViewer`
464: Level: intermediate
466: .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
467: @*/
468: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
469: {
470: PetscFunctionBegin;
472: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer));
474: PetscCall(DMLabelMakeAllValid_Private(label));
475: PetscUseTypeMethod(label, view, viewer);
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: DMLabelViewFromOptions - View a `DMLabel` in a particular way based on a request in the options database
482: Collective
484: Input Parameters:
485: + label - the `DMLabel` object
486: . obj - optional object that provides the prefix for the options database (if `NULL` then the prefix in `obj` is used)
487: - name - option string that is used to activate viewing
489: Options Database Key:
490: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
492: Level: intermediate
494: .seealso: [](ch_dmbase), `DMLabel`, `DMLabelView()`, `PetscObjectViewFromOptions()`, `DMLabelCreate()`
495: @*/
496: PetscErrorCode DMLabelViewFromOptions(DMLabel label, PeOp PetscObject obj, const char name[])
497: {
498: PetscFunctionBegin;
500: PetscCall(PetscObjectViewFromOptions((PetscObject)label, obj, name));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: /*@
505: DMLabelReset - Destroys internal data structures in a `DMLabel`
507: Not Collective
509: Input Parameter:
510: . label - The `DMLabel`
512: Level: beginner
514: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
515: @*/
516: PetscErrorCode DMLabelReset(DMLabel label)
517: {
518: PetscFunctionBegin;
520: for (PetscInt v = 0; v < label->numStrata; ++v) {
521: PetscCall(PetscHSetIDestroy(&label->ht[v]));
522: PetscCall(ISDestroy(&label->points[v]));
523: }
524: label->numStrata = 0;
525: PetscCall(PetscFree(label->stratumValues));
526: PetscCall(PetscFree(label->stratumSizes));
527: PetscCall(PetscFree(label->ht));
528: PetscCall(PetscFree(label->points));
529: PetscCall(PetscFree(label->validIS));
530: PetscCall(PetscHMapIReset(label->hmap));
531: label->pStart = -1;
532: label->pEnd = -1;
533: PetscCall(PetscBTDestroy(&label->bt));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: DMLabelDestroy - Destroys a `DMLabel`
540: Collective
542: Input Parameter:
543: . label - The `DMLabel`
545: Level: beginner
547: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
548: @*/
549: PetscErrorCode DMLabelDestroy(DMLabel *label)
550: {
551: PetscFunctionBegin;
552: if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
554: if (--((PetscObject)*label)->refct > 0) {
555: *label = NULL;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
558: PetscCall(DMLabelReset(*label));
559: PetscCall(PetscHMapIDestroy(&(*label)->hmap));
560: PetscCall(PetscHeaderDestroy(label));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
565: {
566: PetscFunctionBegin;
567: for (PetscInt v = 0; v < label->numStrata; ++v) {
568: PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
569: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
570: (*labelnew)->points[v] = label->points[v];
571: }
572: PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
573: PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: DMLabelDuplicate - Duplicates a `DMLabel`
580: Collective
582: Input Parameter:
583: . label - The `DMLabel`
585: Output Parameter:
586: . labelnew - new label
588: Level: intermediate
590: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
591: @*/
592: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
593: {
594: const char *name;
596: PetscFunctionBegin;
598: PetscCall(DMLabelMakeAllValid_Private(label));
599: PetscCall(PetscObjectGetName((PetscObject)label, &name));
600: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));
602: (*labelnew)->numStrata = label->numStrata;
603: (*labelnew)->defaultValue = label->defaultValue;
604: (*labelnew)->readonly = label->readonly;
605: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
606: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
607: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
608: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
609: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
610: for (PetscInt v = 0; v < label->numStrata; ++v) {
611: (*labelnew)->stratumValues[v] = label->stratumValues[v];
612: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
613: (*labelnew)->validIS[v] = PETSC_TRUE;
614: }
615: (*labelnew)->pStart = -1;
616: (*labelnew)->pEnd = -1;
617: (*labelnew)->bt = NULL;
618: PetscUseTypeMethod(label, duplicate, labelnew);
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /*@
623: DMLabelCompare - Compare two `DMLabel` objects
625: Collective; No Fortran Support
627: Input Parameters:
628: + comm - Comm over which to compare labels
629: . l0 - First `DMLabel`
630: - l1 - Second `DMLabel`
632: Output Parameters:
633: + equal - (Optional) Flag whether the two labels are equal
634: - message - (Optional) Message describing the difference
636: Level: intermediate
638: Notes:
639: The output flag equal is the same on all processes.
640: If it is passed as `NULL` and difference is found, an error is thrown on all processes.
641: Make sure to pass `NULL` on all processes.
643: The output message is set independently on each rank.
644: It is set to `NULL` if no difference was found on the current rank. It must be freed by user.
645: If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner.
646: Make sure to pass `NULL` on all processes.
648: For the comparison, we ignore the order of stratum values, and strata with no points.
650: The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel.
652: Developer Note:
653: Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()`
655: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
656: @*/
657: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
658: {
659: const char *name0, *name1;
660: char msg[PETSC_MAX_PATH_LEN] = "";
661: PetscBool eq;
662: PetscMPIInt rank;
664: PetscFunctionBegin;
667: if (equal) PetscAssertPointer(equal, 4);
668: if (message) PetscAssertPointer(message, 5);
669: PetscCallMPI(MPI_Comm_rank(comm, &rank));
670: PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
671: PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
672: {
673: PetscInt v0, v1;
675: PetscCall(DMLabelGetDefaultValue(l0, &v0));
676: PetscCall(DMLabelGetDefaultValue(l1, &v1));
677: eq = (PetscBool)(v0 == v1);
678: 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));
679: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
680: if (!eq) goto finish;
681: }
682: {
683: IS is0, is1;
685: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
686: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
687: PetscCall(ISEqual(is0, is1, &eq));
688: PetscCall(ISDestroy(&is0));
689: PetscCall(ISDestroy(&is1));
690: if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
691: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
692: if (!eq) goto finish;
693: }
694: {
695: PetscInt nValues;
697: PetscCall(DMLabelGetNumValues(l0, &nValues));
698: for (PetscInt i = 0; i < nValues; i++) {
699: const PetscInt v = l0->stratumValues[i];
700: PetscInt n;
701: IS is0, is1;
703: PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
704: if (!n) continue;
705: PetscCall(DMLabelGetStratumIS(l0, v, &is0));
706: PetscCall(DMLabelGetStratumIS(l1, v, &is1));
707: PetscCall(ISEqualUnsorted(is0, is1, &eq));
708: PetscCall(ISDestroy(&is0));
709: PetscCall(ISDestroy(&is1));
710: if (!eq) {
711: 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));
712: break;
713: }
714: }
715: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
716: }
717: finish:
718: /* If message output arg not set, print to stderr */
719: if (message) {
720: *message = NULL;
721: if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
722: } else {
723: if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
724: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
725: }
726: /* If same output arg not ser and labels are not equal, throw error */
727: if (equal) *equal = eq;
728: else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: /*@
733: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
735: Not Collective
737: Input Parameter:
738: . label - The `DMLabel`
740: Level: intermediate
742: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
743: @*/
744: PetscErrorCode DMLabelComputeIndex(DMLabel label)
745: {
746: PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v;
748: PetscFunctionBegin;
750: PetscCall(DMLabelMakeAllValid_Private(label));
751: for (v = 0; v < label->numStrata; ++v) {
752: const PetscInt *points;
754: PetscCall(ISGetIndices(label->points[v], &points));
755: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
756: const PetscInt point = points[i];
758: pStart = PetscMin(point, pStart);
759: pEnd = PetscMax(point + 1, pEnd);
760: }
761: PetscCall(ISRestoreIndices(label->points[v], &points));
762: }
763: label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
764: label->pEnd = pEnd;
765: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /*@
770: DMLabelCreateIndex - Create an index structure for membership determination
772: Not Collective
774: Input Parameters:
775: + label - The `DMLabel`
776: . pStart - The smallest point
777: - pEnd - The largest point + 1
779: Level: intermediate
781: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
782: @*/
783: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
784: {
785: PetscFunctionBegin;
787: PetscCall(DMLabelDestroyIndex(label));
788: PetscCall(DMLabelMakeAllValid_Private(label));
789: label->pStart = pStart;
790: label->pEnd = pEnd;
791: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
792: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
793: for (PetscInt v = 0; v < label->numStrata; ++v) {
794: IS pointIS;
795: const PetscInt *points;
797: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
798: PetscCall(ISGetIndices(pointIS, &points));
799: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
800: const PetscInt point = points[i];
802: 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);
803: PetscCall(PetscBTSet(label->bt, point - pStart));
804: }
805: PetscCall(ISRestoreIndices(label->points[v], &points));
806: PetscCall(ISDestroy(&pointIS));
807: }
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: /*@
812: DMLabelDestroyIndex - Destroy the index structure
814: Not Collective
816: Input Parameter:
817: . label - the `DMLabel`
819: Level: intermediate
821: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
822: @*/
823: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
824: {
825: PetscFunctionBegin;
827: label->pStart = -1;
828: label->pEnd = -1;
829: PetscCall(PetscBTDestroy(&label->bt));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@
834: DMLabelGetBounds - Return the smallest and largest point in the label
836: Not Collective
838: Input Parameter:
839: . label - the `DMLabel`
841: Output Parameters:
842: + pStart - The smallest point
843: - pEnd - The largest point + 1
845: Level: intermediate
847: Note:
848: This will compute an index for the label if one does not exist.
850: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
851: @*/
852: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
853: {
854: PetscFunctionBegin;
856: if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
857: if (pStart) {
858: PetscAssertPointer(pStart, 2);
859: *pStart = label->pStart;
860: }
861: if (pEnd) {
862: PetscAssertPointer(pEnd, 3);
863: *pEnd = label->pEnd;
864: }
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@
869: DMLabelHasValue - Determine whether a label assigns the value to any point
871: Not Collective
873: Input Parameters:
874: + label - the `DMLabel`
875: - value - the value
877: Output Parameter:
878: . contains - Flag indicating whether the label maps this value to any point
880: Level: developer
882: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
883: @*/
884: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
885: {
886: PetscInt v;
888: PetscFunctionBegin;
890: PetscAssertPointer(contains, 3);
891: PetscCall(DMLabelLookupStratum(label, value, &v));
892: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: /*@
897: DMLabelHasPoint - Determine whether a label assigns a value to a point
899: Not Collective
901: Input Parameters:
902: + label - the `DMLabel`
903: - point - the point
905: Output Parameter:
906: . contains - Flag indicating whether the label maps this point to a value
908: Level: developer
910: Note:
911: The user must call `DMLabelCreateIndex()` before this function.
913: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
914: @*/
915: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
916: {
917: PetscInt pStart, pEnd;
919: PetscFunctionBeginHot;
921: PetscAssertPointer(contains, 3);
922: /* DMLabelGetBounds() calls DMLabelCreateIndex() only if needed */
923: PetscCall(DMLabelGetBounds(label, &pStart, &pEnd));
924: PetscCall(DMLabelMakeAllValid_Private(label));
925: *contains = point >= pStart && point < pEnd && (PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE);
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@
930: DMLabelStratumHasPoint - Return true if the stratum contains a point
932: Not Collective
934: Input Parameters:
935: + label - the `DMLabel`
936: . value - the stratum value
937: - point - the point
939: Output Parameter:
940: . contains - true if the stratum contains the point
942: Level: intermediate
944: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
945: @*/
946: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
947: {
948: PetscFunctionBeginHot;
950: PetscAssertPointer(contains, 4);
951: if (value == label->defaultValue) {
952: PetscInt pointVal;
954: PetscCall(DMLabelGetValue(label, point, &pointVal));
955: *contains = (PetscBool)(pointVal == value);
956: } else {
957: PetscInt v;
959: PetscCall(DMLabelLookupStratum(label, value, &v));
960: if (v >= 0) {
961: if (label->validIS[v] || label->readonly) {
962: IS is;
963: PetscInt i;
965: PetscUseTypeMethod(label, getstratumis, v, &is);
966: PetscCall(ISLocate(is, point, &i));
967: PetscCall(ISDestroy(&is));
968: *contains = (PetscBool)(i >= 0);
969: } else {
970: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
971: }
972: } else { // value is not present
973: *contains = PETSC_FALSE;
974: }
975: }
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*@
980: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
981: When a label is created, it is initialized to -1.
983: Not Collective
985: Input Parameter:
986: . label - a `DMLabel` object
988: Output Parameter:
989: . defaultValue - the default value
991: Level: beginner
993: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
994: @*/
995: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
996: {
997: PetscFunctionBegin;
999: *defaultValue = label->defaultValue;
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
1005: When a label is created, it is initialized to -1.
1007: Not Collective
1009: Input Parameter:
1010: . label - a `DMLabel` object
1012: Output Parameter:
1013: . defaultValue - the default value
1015: Level: beginner
1017: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1018: @*/
1019: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1020: {
1021: PetscFunctionBegin;
1023: label->defaultValue = defaultValue;
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: /*@
1028: 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
1029: `DMLabelSetDefaultValue()`)
1031: Not Collective
1033: Input Parameters:
1034: + label - the `DMLabel`
1035: - point - the point
1037: Output Parameter:
1038: . value - The point value, or the default value (-1 by default)
1040: Level: intermediate
1042: Note:
1043: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1044: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1046: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1047: @*/
1048: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1049: {
1050: PetscInt v;
1052: PetscFunctionBeginHot;
1054: PetscAssertPointer(value, 3);
1055: *value = label->defaultValue;
1056: for (v = 0; v < label->numStrata; ++v) {
1057: if (label->validIS[v] || label->readonly) {
1058: IS is;
1059: PetscInt i;
1061: PetscUseTypeMethod(label, getstratumis, v, &is);
1062: PetscCall(ISLocate(label->points[v], point, &i));
1063: PetscCall(ISDestroy(&is));
1064: if (i >= 0) {
1065: *value = label->stratumValues[v];
1066: break;
1067: }
1068: } else {
1069: PetscBool has;
1071: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1072: if (has) {
1073: *value = label->stratumValues[v];
1074: break;
1075: }
1076: }
1077: }
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: /*@
1082: 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
1083: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1085: Not Collective
1087: Input Parameters:
1088: + label - the `DMLabel`
1089: . point - the point
1090: - value - The point value
1092: Level: intermediate
1094: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1095: @*/
1096: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1097: {
1098: PetscInt v;
1100: PetscFunctionBegin;
1102: /* Find label value, add new entry if needed */
1103: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1104: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1105: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1106: /* Set key */
1107: PetscCall(DMLabelMakeInvalid_Private(label, v));
1108: PetscCall(PetscHSetIAdd(label->ht[v], point));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: /*@
1113: DMLabelClearValue - Clear the value a label assigns to a point
1115: Not Collective
1117: Input Parameters:
1118: + label - the `DMLabel`
1119: . point - the point
1120: - value - The point value
1122: Level: intermediate
1124: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1125: @*/
1126: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1127: {
1128: PetscInt v;
1130: PetscFunctionBegin;
1132: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1133: /* Find label value */
1134: PetscCall(DMLabelLookupStratum(label, value, &v));
1135: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1137: if (label->bt && point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1139: /* Delete key */
1140: PetscCall(DMLabelMakeInvalid_Private(label, v));
1141: PetscCall(PetscHSetIDel(label->ht[v], point));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: /*@
1146: DMLabelInsertIS - Set all points in the `IS` to a value
1148: Not Collective
1150: Input Parameters:
1151: + label - the `DMLabel`
1152: . is - the point `IS`
1153: - value - The point value
1155: Level: intermediate
1157: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1158: @*/
1159: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1160: {
1161: PetscInt v, n, p;
1162: const PetscInt *points;
1164: PetscFunctionBegin;
1167: /* Find label value, add new entry if needed */
1168: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1169: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1170: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1171: /* Set keys */
1172: PetscCall(DMLabelMakeInvalid_Private(label, v));
1173: PetscCall(ISGetLocalSize(is, &n));
1174: PetscCall(ISGetIndices(is, &points));
1175: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1176: PetscCall(ISRestoreIndices(is, &points));
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: /*@
1181: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1183: Not Collective
1185: Input Parameter:
1186: . label - the `DMLabel`
1188: Output Parameter:
1189: . numValues - the number of values
1191: Level: intermediate
1193: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1194: @*/
1195: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1196: {
1197: PetscFunctionBegin;
1199: PetscAssertPointer(numValues, 2);
1200: *numValues = label->numStrata;
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1207: Not Collective
1209: Input Parameter:
1210: . label - the `DMLabel`
1212: Output Parameter:
1213: . values - the value `IS`
1215: Level: intermediate
1217: Notes:
1218: The `values` should be destroyed when no longer needed.
1220: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1222: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1224: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1225: @*/
1226: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1227: {
1228: PetscFunctionBegin;
1230: PetscAssertPointer(values, 2);
1231: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1232: PetscFunctionReturn(PETSC_SUCCESS);
1233: }
1235: /*@
1236: DMLabelGetValueBounds - Return the smallest and largest value in the label
1238: Not Collective
1240: Input Parameter:
1241: . label - the `DMLabel`
1243: Output Parameters:
1244: + minValue - The smallest value
1245: - maxValue - The largest value
1247: Level: intermediate
1249: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1250: @*/
1251: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
1252: {
1253: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;
1255: PetscFunctionBegin;
1257: for (PetscInt v = 0; v < label->numStrata; ++v) {
1258: min = PetscMin(min, label->stratumValues[v]);
1259: max = PetscMax(max, label->stratumValues[v]);
1260: }
1261: if (minValue) {
1262: PetscAssertPointer(minValue, 2);
1263: *minValue = min;
1264: }
1265: if (maxValue) {
1266: PetscAssertPointer(maxValue, 3);
1267: *maxValue = max;
1268: }
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*@
1273: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1275: Not Collective
1277: Input Parameter:
1278: . label - the `DMLabel`
1280: Output Parameter:
1281: . values - the value `IS`
1283: Level: intermediate
1285: Notes:
1286: The `values` should be destroyed when no longer needed.
1288: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1290: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1291: @*/
1292: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1293: {
1294: PetscInt i, j;
1295: PetscInt *valuesArr;
1297: PetscFunctionBegin;
1299: PetscAssertPointer(values, 2);
1300: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1301: for (i = 0, j = 0; i < label->numStrata; i++) {
1302: PetscInt n;
1304: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1305: if (n) valuesArr[j++] = label->stratumValues[i];
1306: }
1307: if (j == label->numStrata) {
1308: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1309: } else {
1310: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1311: }
1312: PetscCall(PetscFree(valuesArr));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: /*@
1317: DMLabelGetValueISGlobal - Get an `IS` of all values that the `DMlabel` takes across all ranks
1319: Collective
1321: Input Parameter:
1322: + comm - MPI communicator to collect values
1323: . label - the `DMLabel`, may be `NULL` for ranks in `comm` which do not have the corresponding `DMLabel`
1324: - get_nonempty - whether to get nonempty stratum values (akin to `DMLabelGetNonEmptyStratumValuesIS()`)
1326: Output Parameter:
1327: . values - the value `IS`
1329: Level: intermediate
1331: Notes:
1332: The `values` should be destroyed when no longer needed.
1334: This is similar to `DMLabelGetValueIS()` and `DMLabelGetNonEmptyStratumValuesIS()`, but gets the (nonempty) values across all ranks in `comm`.
1336: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1337: @*/
1338: PetscErrorCode DMLabelGetValueISGlobal(MPI_Comm comm, DMLabel label, PetscBool get_nonempty, IS *values)
1339: {
1340: PetscInt num_values_local = 0, num_values_global, minmax_values[2], minmax_values_loc[2] = {PETSC_INT_MAX, PETSC_INT_MIN};
1341: IS is_values = NULL;
1342: const PetscInt *values_local = NULL;
1343: PetscInt *values_global;
1345: PetscFunctionBegin;
1347: if (PetscDefined(USE_DEBUG)) {
1349: IS dummy;
1350: PetscCall(ISCreate(comm, &dummy));
1352: PetscCall(ISDestroy(&dummy));
1353: }
1354: PetscAssertPointer(values, 4);
1356: if (label) {
1357: if (get_nonempty) PetscCall(DMLabelGetNonEmptyStratumValuesIS(label, &is_values));
1358: else PetscCall(DMLabelGetValueIS(label, &is_values));
1359: PetscCall(ISGetIndices(is_values, &values_local));
1360: PetscCall(ISGetLocalSize(is_values, &num_values_local));
1361: }
1362: for (PetscInt i = 0; i < num_values_local; i++) {
1363: minmax_values_loc[0] = PetscMin(minmax_values_loc[0], values_local[i]);
1364: minmax_values_loc[1] = PetscMax(minmax_values_loc[1], values_local[i]);
1365: }
1367: PetscCall(PetscGlobalMinMaxInt(comm, minmax_values_loc, minmax_values));
1368: PetscInt value_range = minmax_values[1] - minmax_values[0] + 1;
1369: PetscBT global_values_bt;
1371: // Create a "ballot" where each rank marks which values they have into the PetscBT.
1372: // An Allreduce using bitwise-OR over the ranks then communicates which values are owned by a rank in comm
1373: PetscCall(PetscBTCreate(value_range, &global_values_bt));
1374: for (PetscInt i = 0; i < num_values_local; i++) PetscCall(PetscBTSet(global_values_bt, values_local[i] - minmax_values[0]));
1375: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, global_values_bt, PetscBTLength(value_range), MPI_CHAR, MPI_BOR, comm));
1376: {
1377: PetscCount num_values_global_count;
1378: num_values_global_count = PetscBTCountSet(global_values_bt, value_range);
1379: PetscCall(PetscIntCast(num_values_global_count, &num_values_global));
1380: }
1382: PetscCall(PetscMalloc1(num_values_global, &values_global));
1383: for (PetscInt i = 0, a = 0; i < value_range; i++) {
1384: if (PetscBTLookup(global_values_bt, i)) {
1385: values_global[a] = i + minmax_values[0];
1386: a++;
1387: }
1388: }
1389: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_values_global, values_global, PETSC_OWN_POINTER, values));
1391: PetscCall(PetscBTDestroy(&global_values_bt));
1392: if (is_values) {
1393: PetscCall(ISRestoreIndices(is_values, &values_local));
1394: PetscCall(ISDestroy(&is_values));
1395: }
1396: PetscFunctionReturn(PETSC_SUCCESS);
1397: }
1399: /*@
1400: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1402: Not Collective
1404: Input Parameters:
1405: + label - the `DMLabel`
1406: - value - the value
1408: Output Parameter:
1409: . index - the index of value in the list of values
1411: Level: intermediate
1413: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1414: @*/
1415: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1416: {
1417: PetscInt v;
1419: PetscFunctionBegin;
1421: PetscAssertPointer(index, 3);
1422: /* Do not assume they are sorted */
1423: for (v = 0; v < label->numStrata; ++v)
1424: if (label->stratumValues[v] == value) break;
1425: if (v >= label->numStrata) *index = -1;
1426: else *index = v;
1427: PetscFunctionReturn(PETSC_SUCCESS);
1428: }
1430: /*@
1431: DMLabelHasStratum - Determine whether points exist with the given value
1433: Not Collective
1435: Input Parameters:
1436: + label - the `DMLabel`
1437: - value - the stratum value
1439: Output Parameter:
1440: . exists - Flag saying whether points exist
1442: Level: intermediate
1444: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1445: @*/
1446: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1447: {
1448: PetscInt v;
1450: PetscFunctionBegin;
1452: PetscAssertPointer(exists, 3);
1453: PetscCall(DMLabelLookupStratum(label, value, &v));
1454: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: /*@
1459: DMLabelGetStratumSize - Get the size of a stratum
1461: Not Collective
1463: Input Parameters:
1464: + label - the `DMLabel`
1465: - value - the stratum value
1467: Output Parameter:
1468: . size - The number of points in the stratum
1470: Level: intermediate
1472: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1473: @*/
1474: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1475: {
1476: PetscInt v;
1478: PetscFunctionBegin;
1480: PetscAssertPointer(size, 3);
1481: PetscCall(DMLabelLookupStratum(label, value, &v));
1482: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /*@
1487: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1489: Not Collective
1491: Input Parameters:
1492: + label - the `DMLabel`
1493: - value - the stratum value
1495: Output Parameters:
1496: + start - the smallest point in the stratum
1497: - end - the largest point in the stratum
1499: Level: intermediate
1501: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1502: @*/
1503: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1504: {
1505: IS is;
1506: PetscInt v, min, max;
1508: PetscFunctionBegin;
1510: if (start) {
1511: PetscAssertPointer(start, 3);
1512: *start = -1;
1513: }
1514: if (end) {
1515: PetscAssertPointer(end, 4);
1516: *end = -1;
1517: }
1518: PetscCall(DMLabelLookupStratum(label, value, &v));
1519: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1520: PetscCall(DMLabelMakeValid_Private(label, v));
1521: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1522: PetscUseTypeMethod(label, getstratumis, v, &is);
1523: PetscCall(ISGetMinMax(is, &min, &max));
1524: PetscCall(ISDestroy(&is));
1525: if (start) *start = min;
1526: if (end) *end = max + 1;
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1531: {
1532: PetscFunctionBegin;
1533: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1534: *pointIS = label->points[v];
1535: PetscFunctionReturn(PETSC_SUCCESS);
1536: }
1538: /*@
1539: DMLabelGetStratumIS - Get an `IS` with the stratum points
1541: Not Collective
1543: Input Parameters:
1544: + label - the `DMLabel`
1545: - value - the stratum value
1547: Output Parameter:
1548: . points - The stratum points
1550: Level: intermediate
1552: Notes:
1553: The output `IS` should be destroyed when no longer needed.
1554: Returns `NULL` if the stratum is empty.
1556: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1557: @*/
1558: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1559: {
1560: PetscInt v;
1562: PetscFunctionBegin;
1564: PetscAssertPointer(points, 3);
1565: *points = NULL;
1566: PetscCall(DMLabelLookupStratum(label, value, &v));
1567: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1568: PetscCall(DMLabelMakeValid_Private(label, v));
1569: PetscUseTypeMethod(label, getstratumis, v, points);
1570: PetscFunctionReturn(PETSC_SUCCESS);
1571: }
1573: /*@
1574: DMLabelSetStratumIS - Set the stratum points using an `IS`
1576: Not Collective
1578: Input Parameters:
1579: + label - the `DMLabel`
1580: . value - the stratum value
1581: - is - The stratum points
1583: Level: intermediate
1585: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1586: @*/
1587: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1588: {
1589: PetscInt v;
1591: PetscFunctionBegin;
1594: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1595: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1596: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1597: PetscCall(DMLabelClearStratum(label, value));
1598: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1599: PetscCall(PetscObjectReference((PetscObject)is));
1600: PetscCall(ISDestroy(&label->points[v]));
1601: label->points[v] = is;
1602: label->validIS[v] = PETSC_TRUE;
1603: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1604: if (label->bt) {
1605: const PetscInt *points;
1607: PetscCall(ISGetIndices(is, &points));
1608: for (PetscInt p = 0; p < label->stratumSizes[v]; ++p) {
1609: const PetscInt point = points[p];
1611: 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);
1612: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1613: }
1614: }
1615: PetscFunctionReturn(PETSC_SUCCESS);
1616: }
1618: /*@
1619: DMLabelClearStratum - Remove a stratum
1621: Not Collective
1623: Input Parameters:
1624: + label - the `DMLabel`
1625: - value - the stratum value
1627: Level: intermediate
1629: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1630: @*/
1631: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1632: {
1633: PetscInt v;
1635: PetscFunctionBegin;
1637: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1638: PetscCall(DMLabelLookupStratum(label, value, &v));
1639: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1640: if (label->validIS[v]) {
1641: if (label->bt) {
1642: const PetscInt *points;
1644: PetscCall(ISGetIndices(label->points[v], &points));
1645: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
1646: const PetscInt point = points[i];
1648: if (point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1649: }
1650: PetscCall(ISRestoreIndices(label->points[v], &points));
1651: }
1652: label->stratumSizes[v] = 0;
1653: PetscCall(ISDestroy(&label->points[v]));
1654: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1655: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1656: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1657: } else {
1658: PetscCall(PetscHSetIClear(label->ht[v]));
1659: }
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: /*@
1664: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1666: Not Collective
1668: Input Parameters:
1669: + label - The `DMLabel`
1670: . value - The label value for all points
1671: . pStart - The first point
1672: - pEnd - A point beyond all marked points
1674: Level: intermediate
1676: Note:
1677: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1679: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1680: @*/
1681: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1682: {
1683: IS pIS;
1685: PetscFunctionBegin;
1686: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1687: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1688: PetscCall(ISDestroy(&pIS));
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1695: Not Collective
1697: Input Parameters:
1698: + label - The `DMLabel`
1699: . value - The label value
1700: - p - A point with this value
1702: Output Parameter:
1703: . 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
1705: Level: intermediate
1707: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1708: @*/
1709: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1710: {
1711: IS pointIS;
1712: PetscInt v;
1714: PetscFunctionBegin;
1716: PetscAssertPointer(index, 4);
1717: *index = -1;
1718: PetscCall(DMLabelLookupStratum(label, value, &v));
1719: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1720: PetscCall(DMLabelMakeValid_Private(label, v));
1721: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1722: PetscCall(ISLocate(pointIS, p, index));
1723: PetscCall(ISDestroy(&pointIS));
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }
1727: /*@
1728: DMLabelFilter - Remove all points outside of [`start`, `end`)
1730: Not Collective
1732: Input Parameters:
1733: + label - the `DMLabel`
1734: . start - the first point kept
1735: - end - one more than the last point kept
1737: Level: intermediate
1739: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1740: @*/
1741: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1742: {
1743: PetscInt v;
1745: PetscFunctionBegin;
1747: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1748: PetscCall(DMLabelDestroyIndex(label));
1749: PetscCall(DMLabelMakeAllValid_Private(label));
1750: for (v = 0; v < label->numStrata; ++v) {
1751: PetscCall(ISGeneralFilter(label->points[v], start, end));
1752: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1753: }
1754: PetscCall(DMLabelCreateIndex(label, start, end));
1755: PetscFunctionReturn(PETSC_SUCCESS);
1756: }
1758: /*@
1759: DMLabelPermute - Create a new label with permuted points
1761: Not Collective
1763: Input Parameters:
1764: + label - the `DMLabel`
1765: - permutation - the point permutation
1767: Output Parameter:
1768: . labelNew - the new label containing the permuted points
1770: Level: intermediate
1772: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1773: @*/
1774: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1775: {
1776: const PetscInt *perm;
1777: PetscInt numValues, numPoints, v, q;
1779: PetscFunctionBegin;
1782: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1783: PetscCall(DMLabelMakeAllValid_Private(label));
1784: PetscCall(DMLabelDuplicate(label, labelNew));
1785: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1786: PetscCall(ISGetLocalSize(permutation, &numPoints));
1787: PetscCall(ISGetIndices(permutation, &perm));
1788: for (v = 0; v < numValues; ++v) {
1789: const PetscInt size = (*labelNew)->stratumSizes[v];
1790: const PetscInt *points;
1791: PetscInt *pointsNew;
1793: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1794: PetscCall(PetscCalloc1(size, &pointsNew));
1795: for (q = 0; q < size; ++q) {
1796: const PetscInt point = points[q];
1798: 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);
1799: pointsNew[q] = perm[point];
1800: }
1801: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1802: PetscCall(PetscSortInt(size, pointsNew));
1803: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1804: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1805: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1806: PetscCall(PetscFree(pointsNew));
1807: } else {
1808: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1809: }
1810: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1811: }
1812: PetscCall(ISRestoreIndices(permutation, &perm));
1813: if (label->bt) {
1814: PetscCall(PetscBTDestroy(&label->bt));
1815: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1816: }
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*@
1821: DMLabelPermuteValues - Permute the values in a label
1823: Not collective
1825: Input Parameters:
1826: + label - the `DMLabel`
1827: - permutation - the value permutation, permutation[old value] = new value
1829: Output Parameter:
1830: . label - the `DMLabel` now with permuted values
1832: Note:
1833: The modification is done in-place
1835: Level: intermediate
1837: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1838: @*/
1839: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1840: {
1841: PetscInt Nv, Np;
1843: PetscFunctionBegin;
1846: PetscCall(DMLabelGetNumValues(label, &Nv));
1847: PetscCall(ISGetLocalSize(permutation, &Np));
1848: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1849: if (PetscDefined(USE_DEBUG)) {
1850: PetscBool flg;
1851: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1852: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1853: }
1854: PetscCall(DMLabelRewriteValues(label, permutation));
1855: PetscFunctionReturn(PETSC_SUCCESS);
1856: }
1858: /*@
1859: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1861: Not collective
1863: Input Parameters:
1864: + label - the `DMLabel`
1865: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1867: Output Parameter:
1868: . label - the `DMLabel` now with permuted values
1870: Note:
1871: The modification is done in-place
1873: Level: intermediate
1875: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1876: @*/
1877: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1878: {
1879: const PetscInt *perm;
1880: PetscInt Nv, Np;
1882: PetscFunctionBegin;
1885: PetscCall(DMLabelMakeAllValid_Private(label));
1886: PetscCall(DMLabelGetNumValues(label, &Nv));
1887: PetscCall(ISGetLocalSize(permutation, &Np));
1888: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1889: PetscCall(ISGetIndices(permutation, &perm));
1890: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1891: PetscCall(ISRestoreIndices(permutation, &perm));
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1896: {
1897: MPI_Comm comm;
1898: PetscInt s, l, nroots, nleaves, offset, size;
1899: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1900: PetscSection rootSection;
1901: PetscSF labelSF;
1903: PetscFunctionBegin;
1904: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1905: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1906: /* Build a section of stratum values per point, generate the according SF
1907: and distribute point-wise stratum values to leaves. */
1908: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1909: PetscCall(PetscSectionCreate(comm, &rootSection));
1910: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1911: if (label) {
1912: for (s = 0; s < label->numStrata; ++s) {
1913: const PetscInt *points;
1915: PetscCall(ISGetIndices(label->points[s], &points));
1916: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1917: PetscCall(ISRestoreIndices(label->points[s], &points));
1918: }
1919: }
1920: PetscCall(PetscSectionSetUp(rootSection));
1921: /* Create a point-wise array of stratum values */
1922: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1923: PetscCall(PetscMalloc1(size, &rootStrata));
1924: PetscCall(PetscCalloc1(nroots, &rootIdx));
1925: if (label) {
1926: for (s = 0; s < label->numStrata; ++s) {
1927: const PetscInt *points;
1929: PetscCall(ISGetIndices(label->points[s], &points));
1930: for (l = 0; l < label->stratumSizes[s]; l++) {
1931: const PetscInt p = points[l];
1932: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1933: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1934: }
1935: PetscCall(ISRestoreIndices(label->points[s], &points));
1936: }
1937: }
1938: /* Build SF that maps label points to remote processes */
1939: PetscCall(PetscSectionCreate(comm, leafSection));
1940: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1941: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1942: PetscCall(PetscFree(remoteOffsets));
1943: /* Send the strata for each point over the derived SF */
1944: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1945: PetscCall(PetscMalloc1(size, leafStrata));
1946: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1947: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1948: /* Clean up */
1949: PetscCall(PetscFree(rootStrata));
1950: PetscCall(PetscFree(rootIdx));
1951: PetscCall(PetscSectionDestroy(&rootSection));
1952: PetscCall(PetscSFDestroy(&labelSF));
1953: PetscFunctionReturn(PETSC_SUCCESS);
1954: }
1956: /*@
1957: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1959: Collective
1961: Input Parameters:
1962: + label - the `DMLabel`
1963: - sf - the map from old to new distribution
1965: Output Parameter:
1966: . labelNew - the new redistributed label
1968: Level: intermediate
1970: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1971: @*/
1972: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1973: {
1974: MPI_Comm comm;
1975: PetscSection leafSection;
1976: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1977: PetscInt *leafStrata, *strataIdx;
1978: PetscInt **points;
1979: const char *lname = NULL;
1980: char *name;
1981: PetscMPIInt nameSize;
1982: PetscHSetI stratumHash;
1983: size_t len = 0;
1984: PetscMPIInt rank;
1986: PetscFunctionBegin;
1988: if (label) {
1990: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1991: PetscCall(DMLabelMakeAllValid_Private(label));
1992: }
1993: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1994: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1995: /* Bcast name */
1996: if (rank == 0) {
1997: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1998: PetscCall(PetscStrlen(lname, &len));
1999: }
2000: PetscCall(PetscMPIIntCast(len, &nameSize));
2001: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2002: PetscCall(PetscMalloc1(nameSize + 1, &name));
2003: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2004: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2005: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2006: PetscCall(PetscFree(name));
2007: /* Bcast defaultValue */
2008: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
2009: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
2010: /* Distribute stratum values over the SF and get the point mapping on the receiver */
2011: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
2012: /* Determine received stratum values and initialise new label*/
2013: PetscCall(PetscHSetICreate(&stratumHash));
2014: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
2015: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
2016: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
2017: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
2018: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
2019: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
2020: /* Turn leafStrata into indices rather than stratum values */
2021: offset = 0;
2022: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
2023: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
2024: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
2025: for (p = 0; p < size; ++p) {
2026: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2027: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
2028: leafStrata[p] = s;
2029: break;
2030: }
2031: }
2032: }
2033: /* Rebuild the point strata on the receiver */
2034: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
2035: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2036: for (p = pStart; p < pEnd; p++) {
2037: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2038: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2039: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
2040: }
2041: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
2042: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
2043: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
2044: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2045: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
2046: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
2047: }
2048: /* Insert points into new strata */
2049: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
2050: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2051: for (p = pStart; p < pEnd; p++) {
2052: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2053: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2054: for (s = 0; s < dof; s++) {
2055: stratum = leafStrata[offset + s];
2056: points[stratum][strataIdx[stratum]++] = p;
2057: }
2058: }
2059: for (s = 0; s < (*labelNew)->numStrata; s++) {
2060: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
2061: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
2062: }
2063: PetscCall(PetscFree(points));
2064: PetscCall(PetscHSetIDestroy(&stratumHash));
2065: PetscCall(PetscFree(leafStrata));
2066: PetscCall(PetscFree(strataIdx));
2067: PetscCall(PetscSectionDestroy(&leafSection));
2068: PetscFunctionReturn(PETSC_SUCCESS);
2069: }
2071: /*@
2072: DMLabelGather - Gather all label values from leafs into roots
2074: Collective
2076: Input Parameters:
2077: + label - the `DMLabel`
2078: - sf - the `PetscSF` communication map
2080: Output Parameter:
2081: . labelNew - the new `DMLabel` with localised leaf values
2083: Level: developer
2085: Note:
2086: This is the inverse operation to `DMLabelDistribute()`.
2088: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2089: @*/
2090: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
2091: {
2092: MPI_Comm comm;
2093: PetscSection rootSection;
2094: PetscSF sfLabel;
2095: PetscSFNode *rootPoints, *leafPoints;
2096: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2097: const PetscInt *rootDegree, *ilocal;
2098: PetscInt *rootStrata;
2099: const char *lname;
2100: char *name;
2101: PetscMPIInt nameSize;
2102: size_t len = 0;
2103: PetscMPIInt rank, size;
2105: PetscFunctionBegin;
2108: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2109: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2110: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2111: PetscCallMPI(MPI_Comm_size(comm, &size));
2112: /* Bcast name */
2113: if (rank == 0) {
2114: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2115: PetscCall(PetscStrlen(lname, &len));
2116: }
2117: PetscCall(PetscMPIIntCast(len, &nameSize));
2118: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2119: PetscCall(PetscMalloc1(nameSize + 1, &name));
2120: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2121: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2122: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2123: PetscCall(PetscFree(name));
2124: /* Gather rank/index pairs of leaves into local roots to build
2125: an inverse, multi-rooted SF. Note that this ignores local leaf
2126: indexing due to the use of the multiSF in PetscSFGather. */
2127: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2128: PetscCall(PetscMalloc1(nroots, &leafPoints));
2129: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2130: for (p = 0; p < nleaves; p++) {
2131: PetscInt ilp = ilocal ? ilocal[p] : p;
2133: leafPoints[ilp].index = ilp;
2134: leafPoints[ilp].rank = rank;
2135: }
2136: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2137: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2138: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2139: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2140: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2141: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2142: PetscCall(PetscSFCreate(comm, &sfLabel));
2143: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2144: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2145: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2146: /* Rebuild the point strata on the receiver */
2147: for (p = 0, idx = 0; p < nroots; p++) {
2148: for (d = 0; d < rootDegree[p]; d++) {
2149: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2150: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2151: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2152: }
2153: idx += rootDegree[p];
2154: }
2155: PetscCall(PetscFree(leafPoints));
2156: PetscCall(PetscFree(rootStrata));
2157: PetscCall(PetscSectionDestroy(&rootSection));
2158: PetscCall(PetscSFDestroy(&sfLabel));
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2163: {
2164: const PetscInt *degree;
2165: const PetscInt *points;
2166: PetscInt Nr, r, Nl, l, val, defVal;
2168: PetscFunctionBegin;
2169: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2170: /* Add in leaves */
2171: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2172: for (l = 0; l < Nl; ++l) {
2173: PetscCall(DMLabelGetValue(label, points[l], &val));
2174: if (val != defVal) valArray[points[l]] = val;
2175: }
2176: /* Add in shared roots */
2177: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2178: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2179: for (r = 0; r < Nr; ++r) {
2180: if (degree[r]) {
2181: PetscCall(DMLabelGetValue(label, r, &val));
2182: if (val != defVal) valArray[r] = val;
2183: }
2184: }
2185: PetscFunctionReturn(PETSC_SUCCESS);
2186: }
2188: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), PetscCtx ctx)
2189: {
2190: const PetscInt *degree;
2191: const PetscInt *points;
2192: PetscInt Nr, r, Nl, l, val, defVal;
2194: PetscFunctionBegin;
2195: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2196: /* Read out leaves */
2197: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2198: for (l = 0; l < Nl; ++l) {
2199: const PetscInt p = points[l];
2200: const PetscInt cval = valArray[p];
2202: if (cval != defVal) {
2203: PetscCall(DMLabelGetValue(label, p, &val));
2204: if (val == defVal) {
2205: PetscCall(DMLabelSetValue(label, p, cval));
2206: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2207: }
2208: }
2209: }
2210: /* Read out shared roots */
2211: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2212: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2213: for (r = 0; r < Nr; ++r) {
2214: if (degree[r]) {
2215: const PetscInt cval = valArray[r];
2217: if (cval != defVal) {
2218: PetscCall(DMLabelGetValue(label, r, &val));
2219: if (val == defVal) {
2220: PetscCall(DMLabelSetValue(label, r, cval));
2221: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2222: }
2223: }
2224: }
2225: }
2226: PetscFunctionReturn(PETSC_SUCCESS);
2227: }
2229: /*@
2230: DMLabelPropagateBegin - Setup a cycle of label propagation
2232: Collective
2234: Input Parameters:
2235: + label - The `DMLabel` to propagate across processes
2236: - sf - The `PetscSF` describing parallel layout of the label points
2238: Level: intermediate
2240: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2241: @*/
2242: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2243: {
2244: PetscInt Nr, r, defVal;
2245: PetscMPIInt size;
2247: PetscFunctionBegin;
2248: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2249: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2250: if (size > 1) {
2251: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2252: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2253: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2254: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2255: }
2256: PetscFunctionReturn(PETSC_SUCCESS);
2257: }
2259: /*@
2260: DMLabelPropagateEnd - Tear down a cycle of label propagation
2262: Collective
2264: Input Parameters:
2265: + label - The `DMLabel` to propagate across processes
2266: - pointSF - The `PetscSF` describing parallel layout of the label points
2268: Level: intermediate
2270: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2271: @*/
2272: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2273: {
2274: PetscFunctionBegin;
2275: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2276: PetscCall(PetscFree(label->propArray));
2277: label->propArray = NULL;
2278: PetscFunctionReturn(PETSC_SUCCESS);
2279: }
2281: /*@
2282: DMLabelPropagatePush - Execute a cycle of label propagation
2284: Collective
2286: Input Parameters:
2287: + label - The `DMLabel` to propagate across processes
2288: . pointSF - The `PetscSF` describing parallel layout of the label points
2289: . merge - The operator which merges label values
2290: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2291: - ctx - An optional application context for the callback, or `NULL`
2293: Calling sequence of `markPoint`:
2294: + label - The `DMLabel`
2295: . p - The point being marked
2296: . val - The label value for `p`
2297: - ctx - An optional application context
2299: Level: intermediate
2301: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2302: @*/
2303: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, MPI_Op merge, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, PetscCtx ctx), PetscCtx ctx)
2304: {
2305: PetscInt *valArray = label->propArray, Nr;
2306: PetscMPIInt size;
2308: PetscFunctionBegin;
2309: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2310: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2311: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2312: if (size > 1 && Nr >= 0) {
2313: /* Communicate marked edges
2314: The current implementation allocates an array the size of the number of root. We put the label values into the
2315: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2317: TODO: We could use in-place communication with a different SF
2318: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2319: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2321: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2322: 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
2323: edge to the queue.
2324: */
2325: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2326: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, merge));
2327: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, merge));
2328: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2329: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2330: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2331: }
2332: PetscFunctionReturn(PETSC_SUCCESS);
2333: }
2335: /*@
2336: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2338: Not Collective
2340: Input Parameter:
2341: . label - the `DMLabel`
2343: Output Parameters:
2344: + section - the section giving offsets for each stratum
2345: - is - An `IS` containing all the label points
2347: Level: developer
2349: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2350: @*/
2351: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2352: {
2353: IS vIS;
2354: const PetscInt *values;
2355: PetscInt *points;
2356: PetscInt nV, vS = 0, vE = 0, v, N;
2358: PetscFunctionBegin;
2360: PetscCall(DMLabelGetNumValues(label, &nV));
2361: PetscCall(DMLabelGetValueIS(label, &vIS));
2362: PetscCall(ISGetIndices(vIS, &values));
2363: if (nV) {
2364: vS = values[0];
2365: vE = values[0] + 1;
2366: }
2367: for (v = 1; v < nV; ++v) {
2368: vS = PetscMin(vS, values[v]);
2369: vE = PetscMax(vE, values[v] + 1);
2370: }
2371: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2372: PetscCall(PetscSectionSetChart(*section, vS, vE));
2373: for (v = 0; v < nV; ++v) {
2374: PetscInt n;
2376: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2377: PetscCall(PetscSectionSetDof(*section, values[v], n));
2378: }
2379: PetscCall(PetscSectionSetUp(*section));
2380: PetscCall(PetscSectionGetStorageSize(*section, &N));
2381: PetscCall(PetscMalloc1(N, &points));
2382: for (v = 0; v < nV; ++v) {
2383: IS is;
2384: const PetscInt *spoints;
2385: PetscInt dof, off, p;
2387: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2388: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2389: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2390: PetscCall(ISGetIndices(is, &spoints));
2391: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2392: PetscCall(ISRestoreIndices(is, &spoints));
2393: PetscCall(ISDestroy(&is));
2394: }
2395: PetscCall(ISRestoreIndices(vIS, &values));
2396: PetscCall(ISDestroy(&vIS));
2397: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2398: PetscFunctionReturn(PETSC_SUCCESS);
2399: }
2401: /*@
2402: DMLabelRegister - Adds a new label component implementation
2404: Not Collective
2406: Input Parameters:
2407: + name - The name of a new user-defined creation routine
2408: - create_func - The creation routine itself
2410: Notes:
2411: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2413: Example Usage:
2414: .vb
2415: DMLabelRegister("my_label", MyLabelCreate);
2416: .ve
2418: Then, your label type can be chosen with the procedural interface via
2419: .vb
2420: DMLabelCreate(MPI_Comm, DMLabel *);
2421: DMLabelSetType(DMLabel, "my_label");
2422: .ve
2423: or at runtime via the option
2424: .vb
2425: -dm_label_type my_label
2426: .ve
2428: Level: advanced
2430: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2431: @*/
2432: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2433: {
2434: PetscFunctionBegin;
2435: PetscCall(DMInitializePackage());
2436: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2437: PetscFunctionReturn(PETSC_SUCCESS);
2438: }
2440: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2441: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2443: /*@
2444: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2446: Not Collective
2448: Level: advanced
2450: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2451: @*/
2452: PetscErrorCode DMLabelRegisterAll(void)
2453: {
2454: PetscFunctionBegin;
2455: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2456: DMLabelRegisterAllCalled = PETSC_TRUE;
2458: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2459: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2460: PetscFunctionReturn(PETSC_SUCCESS);
2461: }
2463: /*@
2464: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2466: Level: developer
2468: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2469: @*/
2470: PetscErrorCode DMLabelRegisterDestroy(void)
2471: {
2472: PetscFunctionBegin;
2473: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2474: DMLabelRegisterAllCalled = PETSC_FALSE;
2475: PetscFunctionReturn(PETSC_SUCCESS);
2476: }
2478: /*@
2479: DMLabelSetType - Sets the particular implementation for a label.
2481: Collective
2483: Input Parameters:
2484: + label - The label
2485: - method - The name of the label type
2487: Options Database Key:
2488: . -dm_label_type type - Sets the label type; see `DMLabelType`
2490: Level: intermediate
2492: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`, `DMLabelType`
2493: @*/
2494: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2495: {
2496: PetscErrorCode (*r)(DMLabel);
2497: PetscBool match;
2499: PetscFunctionBegin;
2501: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2502: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2504: PetscCall(DMLabelRegisterAll());
2505: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2506: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2508: PetscTryTypeMethod(label, destroy);
2509: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2510: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2511: PetscCall((*r)(label));
2512: PetscFunctionReturn(PETSC_SUCCESS);
2513: }
2515: /*@
2516: DMLabelGetType - Gets the type name (as a string) from the label.
2518: Not Collective
2520: Input Parameter:
2521: . label - The `DMLabel`
2523: Output Parameter:
2524: . type - The `DMLabel` type name
2526: Level: intermediate
2528: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2529: @*/
2530: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2531: {
2532: PetscFunctionBegin;
2534: PetscAssertPointer(type, 2);
2535: PetscCall(DMLabelRegisterAll());
2536: *type = ((PetscObject)label)->type_name;
2537: PetscFunctionReturn(PETSC_SUCCESS);
2538: }
2540: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2541: {
2542: PetscFunctionBegin;
2543: label->ops->view = DMLabelView_Concrete;
2544: label->ops->setup = NULL;
2545: label->ops->duplicate = DMLabelDuplicate_Concrete;
2546: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2550: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2551: {
2552: PetscFunctionBegin;
2554: PetscCall(DMLabelInitialize_Concrete(label));
2555: PetscFunctionReturn(PETSC_SUCCESS);
2556: }
2558: /*@
2559: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2560: the local section and an `PetscSF` describing the section point overlap.
2562: Collective
2564: Input Parameters:
2565: + s - The `PetscSection` for the local field layout
2566: . sf - The `PetscSF` describing parallel layout of the section points
2567: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2568: . label - The label specifying the points
2569: - labelValue - The label stratum specifying the points
2571: Output Parameter:
2572: . gsection - The `PetscSection` for the global field layout
2574: Level: developer
2576: Note:
2577: This gives negative sizes and offsets to points not owned by this process
2579: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2580: @*/
2581: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2582: {
2583: PetscInt *neg = NULL, *tmpOff = NULL;
2584: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2586: PetscFunctionBegin;
2590: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2591: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2592: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2593: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2594: if (nroots >= 0) {
2595: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2596: PetscCall(PetscCalloc1(nroots, &neg));
2597: if (nroots > pEnd - pStart) {
2598: PetscCall(PetscCalloc1(nroots, &tmpOff));
2599: } else {
2600: tmpOff = &(*gsection)->atlasDof[-pStart];
2601: }
2602: }
2603: /* Mark ghost points with negative dof */
2604: for (p = pStart; p < pEnd; ++p) {
2605: PetscInt value;
2607: PetscCall(DMLabelGetValue(label, p, &value));
2608: if (value != labelValue) continue;
2609: PetscCall(PetscSectionGetDof(s, p, &dof));
2610: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2611: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2612: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2613: if (neg) neg[p] = -(dof + 1);
2614: }
2615: PetscCall(PetscSectionSetUpBC(*gsection));
2616: if (nroots >= 0) {
2617: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2618: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2619: if (nroots > pEnd - pStart) {
2620: for (p = pStart; p < pEnd; ++p) {
2621: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2622: }
2623: }
2624: }
2625: /* Calculate new sizes, get process offset, and calculate point offsets */
2626: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2627: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2628: (*gsection)->atlasOff[p] = off;
2629: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2630: }
2631: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2632: globalOff -= off;
2633: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2634: (*gsection)->atlasOff[p] += globalOff;
2635: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2636: }
2637: /* Put in negative offsets for ghost points */
2638: if (nroots >= 0) {
2639: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2640: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2641: if (nroots > pEnd - pStart) {
2642: for (p = pStart; p < pEnd; ++p) {
2643: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2644: }
2645: }
2646: }
2647: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2648: PetscCall(PetscFree(neg));
2649: PetscFunctionReturn(PETSC_SUCCESS);
2650: }
2652: typedef struct _n_PetscSectionSym_Label {
2653: DMLabel label;
2654: PetscCopyMode *modes;
2655: PetscInt *sizes;
2656: const PetscInt ***perms;
2657: const PetscScalar ***rots;
2658: PetscInt (*minMaxOrients)[2];
2659: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2660: } PetscSectionSym_Label;
2662: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2663: {
2664: PetscInt i, j;
2665: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2667: PetscFunctionBegin;
2668: for (i = 0; i <= sl->numStrata; i++) {
2669: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2670: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2671: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2672: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2673: }
2674: if (sl->perms[i]) {
2675: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2677: PetscCall(PetscFree(perms));
2678: }
2679: if (sl->rots[i]) {
2680: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2682: PetscCall(PetscFree(rots));
2683: }
2684: }
2685: }
2686: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2687: PetscCall(DMLabelDestroy(&sl->label));
2688: sl->numStrata = 0;
2689: PetscFunctionReturn(PETSC_SUCCESS);
2690: }
2692: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2693: {
2694: PetscFunctionBegin;
2695: PetscCall(PetscSectionSymLabelReset(sym));
2696: PetscCall(PetscFree(sym->data));
2697: PetscFunctionReturn(PETSC_SUCCESS);
2698: }
2700: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2701: {
2702: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2703: PetscBool isAscii;
2704: DMLabel label = sl->label;
2705: const char *name;
2707: PetscFunctionBegin;
2708: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2709: if (isAscii) {
2710: PetscInt i, j, k;
2711: PetscViewerFormat format;
2713: PetscCall(PetscViewerGetFormat(viewer, &format));
2714: if (label) {
2715: PetscCall(PetscViewerGetFormat(viewer, &format));
2716: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2717: PetscCall(PetscViewerASCIIPushTab(viewer));
2718: PetscCall(DMLabelView(label, viewer));
2719: PetscCall(PetscViewerASCIIPopTab(viewer));
2720: } else {
2721: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2722: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2723: }
2724: } else {
2725: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2726: }
2727: PetscCall(PetscViewerASCIIPushTab(viewer));
2728: for (i = 0; i <= sl->numStrata; i++) {
2729: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2731: if (!(sl->perms[i] || sl->rots[i])) {
2732: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2733: } else {
2734: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2735: PetscCall(PetscViewerASCIIPushTab(viewer));
2736: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2737: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2738: PetscCall(PetscViewerASCIIPushTab(viewer));
2739: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2740: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2741: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2742: } else {
2743: PetscInt tab;
2745: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2746: PetscCall(PetscViewerASCIIPushTab(viewer));
2747: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2748: if (sl->perms[i] && sl->perms[i][j]) {
2749: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2750: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2751: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2752: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2753: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2754: }
2755: if (sl->rots[i] && sl->rots[i][j]) {
2756: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2757: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2758: #if PetscDefined(USE_COMPLEX)
2759: 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])));
2760: #else
2761: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2762: #endif
2763: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2764: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2765: }
2766: PetscCall(PetscViewerASCIIPopTab(viewer));
2767: }
2768: }
2769: PetscCall(PetscViewerASCIIPopTab(viewer));
2770: }
2771: PetscCall(PetscViewerASCIIPopTab(viewer));
2772: }
2773: }
2774: PetscCall(PetscViewerASCIIPopTab(viewer));
2775: }
2776: PetscFunctionReturn(PETSC_SUCCESS);
2777: }
2779: /*@
2780: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2782: Logically
2784: Input Parameters:
2785: + sym - the section symmetries
2786: - label - the `DMLabel` describing the types of points
2788: Level: developer:
2790: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2791: @*/
2792: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2793: {
2794: PetscSectionSym_Label *sl;
2796: PetscFunctionBegin;
2798: sl = (PetscSectionSym_Label *)sym->data;
2799: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2800: if (label) {
2801: sl->label = label;
2802: PetscCall(PetscObjectReference((PetscObject)label));
2803: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2804: 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));
2805: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2806: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2807: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2808: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2809: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2810: }
2811: PetscFunctionReturn(PETSC_SUCCESS);
2812: }
2814: /*@
2815: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2817: Logically Collective
2819: Input Parameters:
2820: + sym - the section symmetries
2821: - stratum - the stratum value in the label that we are assigning symmetries for
2823: Output Parameters:
2824: + size - the number of dofs for points in the `stratum` of the label
2825: . minOrient - the smallest orientation for a point in this `stratum`
2826: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2827: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2828: - 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
2830: Level: developer
2832: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2833: @*/
2834: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2835: {
2836: PetscSectionSym_Label *sl;
2837: const char *name;
2838: PetscInt i;
2840: PetscFunctionBegin;
2842: sl = (PetscSectionSym_Label *)sym->data;
2843: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2844: for (i = 0; i <= sl->numStrata; i++) {
2845: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2847: if (stratum == value) break;
2848: }
2849: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2850: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2851: if (size) {
2852: PetscAssertPointer(size, 3);
2853: *size = sl->sizes[i];
2854: }
2855: if (minOrient) {
2856: PetscAssertPointer(minOrient, 4);
2857: *minOrient = sl->minMaxOrients[i][0];
2858: }
2859: if (maxOrient) {
2860: PetscAssertPointer(maxOrient, 5);
2861: *maxOrient = sl->minMaxOrients[i][1];
2862: }
2863: if (perms) {
2864: PetscAssertPointer(perms, 6);
2865: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2866: }
2867: if (rots) {
2868: PetscAssertPointer(rots, 7);
2869: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2870: }
2871: PetscFunctionReturn(PETSC_SUCCESS);
2872: }
2874: /*@
2875: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2877: Logically
2879: Input Parameters:
2880: + sym - the section symmetries
2881: . stratum - the stratum value in the label that we are assigning symmetries for
2882: . size - the number of dofs for points in the `stratum` of the label
2883: . minOrient - the smallest orientation for a point in this `stratum`
2884: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2885: . mode - how `sym` should copy the `perms` and `rots` arrays
2886: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2887: - 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
2889: Level: developer
2891: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2892: @*/
2893: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2894: {
2895: PetscSectionSym_Label *sl;
2896: const char *name;
2897: PetscInt i, j, k;
2899: PetscFunctionBegin;
2901: sl = (PetscSectionSym_Label *)sym->data;
2902: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2903: for (i = 0; i <= sl->numStrata; i++) {
2904: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2906: if (stratum == value) break;
2907: }
2908: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2909: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2910: sl->sizes[i] = size;
2911: sl->modes[i] = mode;
2912: sl->minMaxOrients[i][0] = minOrient;
2913: sl->minMaxOrients[i][1] = maxOrient;
2914: if (mode == PETSC_COPY_VALUES) {
2915: if (perms) {
2916: PetscInt **ownPerms;
2918: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2919: for (j = 0; j < maxOrient - minOrient; j++) {
2920: if (perms[j]) {
2921: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2922: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2923: }
2924: }
2925: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2926: }
2927: if (rots) {
2928: PetscScalar **ownRots;
2930: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2931: for (j = 0; j < maxOrient - minOrient; j++) {
2932: if (rots[j]) {
2933: PetscCall(PetscMalloc1(size, &ownRots[j]));
2934: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2935: }
2936: }
2937: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2938: }
2939: } else {
2940: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2941: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2942: }
2943: PetscFunctionReturn(PETSC_SUCCESS);
2944: }
2946: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2947: {
2948: PetscInt i, j, numStrata;
2949: PetscSectionSym_Label *sl;
2950: DMLabel label;
2952: PetscFunctionBegin;
2953: sl = (PetscSectionSym_Label *)sym->data;
2954: numStrata = sl->numStrata;
2955: label = sl->label;
2956: for (i = 0; i < numPoints; i++) {
2957: PetscInt point = points[2 * i];
2958: PetscInt ornt = points[2 * i + 1];
2960: for (j = 0; j < numStrata; j++) {
2961: if (label->validIS[j]) {
2962: PetscInt k;
2964: PetscCall(ISLocate(label->points[j], point, &k));
2965: if (k >= 0) break;
2966: } else {
2967: PetscBool has;
2969: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2970: if (has) break;
2971: }
2972: }
2973: 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],
2974: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2975: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2976: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2977: }
2978: PetscFunctionReturn(PETSC_SUCCESS);
2979: }
2981: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2982: {
2983: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2984: IS valIS;
2985: const PetscInt *values;
2986: PetscInt Nv;
2988: PetscFunctionBegin;
2989: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2990: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2991: PetscCall(ISGetIndices(valIS, &values));
2992: for (PetscInt v = 0; v < Nv; ++v) {
2993: const PetscInt val = values[v];
2994: PetscInt size, minOrient, maxOrient;
2995: const PetscInt **perms;
2996: const PetscScalar **rots;
2998: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2999: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
3000: }
3001: PetscCall(ISDestroy(&valIS));
3002: PetscFunctionReturn(PETSC_SUCCESS);
3003: }
3005: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3006: {
3007: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
3008: DMLabel dlabel;
3010: PetscFunctionBegin;
3011: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
3012: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
3013: PetscCall(DMLabelDestroy(&dlabel));
3014: PetscCall(PetscSectionSymCopy(sym, *dsym));
3015: PetscFunctionReturn(PETSC_SUCCESS);
3016: }
3018: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
3019: {
3020: PetscSectionSym_Label *sl;
3022: PetscFunctionBegin;
3023: PetscCall(PetscNew(&sl));
3024: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
3025: sym->ops->distribute = PetscSectionSymDistribute_Label;
3026: sym->ops->copy = PetscSectionSymCopy_Label;
3027: sym->ops->view = PetscSectionSymView_Label;
3028: sym->ops->destroy = PetscSectionSymDestroy_Label;
3029: sym->data = (void *)sl;
3030: PetscFunctionReturn(PETSC_SUCCESS);
3031: }
3033: /*@
3034: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
3036: Collective
3038: Input Parameters:
3039: + comm - the MPI communicator for the new symmetry
3040: - label - the label defining the strata
3042: Output Parameter:
3043: . sym - the section symmetries
3045: Level: developer
3047: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3048: @*/
3049: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
3050: {
3051: PetscFunctionBegin;
3052: PetscCall(DMInitializePackage());
3053: PetscCall(PetscSectionSymCreate(comm, sym));
3054: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
3055: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
3056: PetscFunctionReturn(PETSC_SUCCESS);
3057: }