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: /*@C
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 local_values_bt, 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, &local_values_bt));
1374: PetscCall(PetscBTCreate(value_range, &global_values_bt));
1375: for (PetscInt i = 0; i < num_values_local; i++) PetscCall(PetscBTSet(local_values_bt, values_local[i] - minmax_values[0]));
1376: PetscCallMPI(MPIU_Allreduce(local_values_bt, global_values_bt, PetscBTLength(value_range), MPI_CHAR, MPI_BOR, comm));
1377: PetscCall(PetscBTDestroy(&local_values_bt));
1378: {
1379: PetscCount num_values_global_count;
1380: num_values_global_count = PetscBTCountSet(global_values_bt, value_range);
1381: PetscCall(PetscIntCast(num_values_global_count, &num_values_global));
1382: }
1384: PetscCall(PetscMalloc1(num_values_global, &values_global));
1385: for (PetscInt i = 0, a = 0; i < value_range; i++) {
1386: if (PetscBTLookup(global_values_bt, i)) {
1387: values_global[a] = i + minmax_values[0];
1388: a++;
1389: }
1390: }
1391: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_values_global, values_global, PETSC_OWN_POINTER, values));
1393: PetscCall(PetscBTDestroy(&global_values_bt));
1394: if (is_values) {
1395: PetscCall(ISRestoreIndices(is_values, &values_local));
1396: PetscCall(ISDestroy(&is_values));
1397: }
1398: PetscFunctionReturn(PETSC_SUCCESS);
1399: }
1401: /*@
1402: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1404: Not Collective
1406: Input Parameters:
1407: + label - the `DMLabel`
1408: - value - the value
1410: Output Parameter:
1411: . index - the index of value in the list of values
1413: Level: intermediate
1415: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1416: @*/
1417: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1418: {
1419: PetscInt v;
1421: PetscFunctionBegin;
1423: PetscAssertPointer(index, 3);
1424: /* Do not assume they are sorted */
1425: for (v = 0; v < label->numStrata; ++v)
1426: if (label->stratumValues[v] == value) break;
1427: if (v >= label->numStrata) *index = -1;
1428: else *index = v;
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: /*@
1433: DMLabelHasStratum - Determine whether points exist with the given value
1435: Not Collective
1437: Input Parameters:
1438: + label - the `DMLabel`
1439: - value - the stratum value
1441: Output Parameter:
1442: . exists - Flag saying whether points exist
1444: Level: intermediate
1446: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1447: @*/
1448: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1449: {
1450: PetscInt v;
1452: PetscFunctionBegin;
1454: PetscAssertPointer(exists, 3);
1455: PetscCall(DMLabelLookupStratum(label, value, &v));
1456: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: DMLabelGetStratumSize - Get the size of a stratum
1463: Not Collective
1465: Input Parameters:
1466: + label - the `DMLabel`
1467: - value - the stratum value
1469: Output Parameter:
1470: . size - The number of points in the stratum
1472: Level: intermediate
1474: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1475: @*/
1476: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1477: {
1478: PetscInt v;
1480: PetscFunctionBegin;
1482: PetscAssertPointer(size, 3);
1483: PetscCall(DMLabelLookupStratum(label, value, &v));
1484: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1485: PetscFunctionReturn(PETSC_SUCCESS);
1486: }
1488: /*@
1489: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1491: Not Collective
1493: Input Parameters:
1494: + label - the `DMLabel`
1495: - value - the stratum value
1497: Output Parameters:
1498: + start - the smallest point in the stratum
1499: - end - the largest point in the stratum
1501: Level: intermediate
1503: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1504: @*/
1505: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1506: {
1507: IS is;
1508: PetscInt v, min, max;
1510: PetscFunctionBegin;
1512: if (start) {
1513: PetscAssertPointer(start, 3);
1514: *start = -1;
1515: }
1516: if (end) {
1517: PetscAssertPointer(end, 4);
1518: *end = -1;
1519: }
1520: PetscCall(DMLabelLookupStratum(label, value, &v));
1521: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1522: PetscCall(DMLabelMakeValid_Private(label, v));
1523: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1524: PetscUseTypeMethod(label, getstratumis, v, &is);
1525: PetscCall(ISGetMinMax(is, &min, &max));
1526: PetscCall(ISDestroy(&is));
1527: if (start) *start = min;
1528: if (end) *end = max + 1;
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1533: {
1534: PetscFunctionBegin;
1535: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1536: *pointIS = label->points[v];
1537: PetscFunctionReturn(PETSC_SUCCESS);
1538: }
1540: /*@
1541: DMLabelGetStratumIS - Get an `IS` with the stratum points
1543: Not Collective
1545: Input Parameters:
1546: + label - the `DMLabel`
1547: - value - the stratum value
1549: Output Parameter:
1550: . points - The stratum points
1552: Level: intermediate
1554: Notes:
1555: The output `IS` should be destroyed when no longer needed.
1556: Returns `NULL` if the stratum is empty.
1558: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1559: @*/
1560: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1561: {
1562: PetscInt v;
1564: PetscFunctionBegin;
1566: PetscAssertPointer(points, 3);
1567: *points = NULL;
1568: PetscCall(DMLabelLookupStratum(label, value, &v));
1569: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1570: PetscCall(DMLabelMakeValid_Private(label, v));
1571: PetscUseTypeMethod(label, getstratumis, v, points);
1572: PetscFunctionReturn(PETSC_SUCCESS);
1573: }
1575: /*@
1576: DMLabelSetStratumIS - Set the stratum points using an `IS`
1578: Not Collective
1580: Input Parameters:
1581: + label - the `DMLabel`
1582: . value - the stratum value
1583: - is - The stratum points
1585: Level: intermediate
1587: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1588: @*/
1589: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1590: {
1591: PetscInt v;
1593: PetscFunctionBegin;
1596: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1597: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1598: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1599: PetscCall(DMLabelClearStratum(label, value));
1600: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1601: PetscCall(PetscObjectReference((PetscObject)is));
1602: PetscCall(ISDestroy(&label->points[v]));
1603: label->points[v] = is;
1604: label->validIS[v] = PETSC_TRUE;
1605: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1606: if (label->bt) {
1607: const PetscInt *points;
1609: PetscCall(ISGetIndices(is, &points));
1610: for (PetscInt p = 0; p < label->stratumSizes[v]; ++p) {
1611: const PetscInt point = points[p];
1613: 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);
1614: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1615: }
1616: }
1617: PetscFunctionReturn(PETSC_SUCCESS);
1618: }
1620: /*@
1621: DMLabelClearStratum - Remove a stratum
1623: Not Collective
1625: Input Parameters:
1626: + label - the `DMLabel`
1627: - value - the stratum value
1629: Level: intermediate
1631: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1632: @*/
1633: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1634: {
1635: PetscInt v;
1637: PetscFunctionBegin;
1639: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1640: PetscCall(DMLabelLookupStratum(label, value, &v));
1641: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1642: if (label->validIS[v]) {
1643: if (label->bt) {
1644: const PetscInt *points;
1646: PetscCall(ISGetIndices(label->points[v], &points));
1647: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
1648: const PetscInt point = points[i];
1650: if (point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1651: }
1652: PetscCall(ISRestoreIndices(label->points[v], &points));
1653: }
1654: label->stratumSizes[v] = 0;
1655: PetscCall(ISDestroy(&label->points[v]));
1656: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1657: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1658: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1659: } else {
1660: PetscCall(PetscHSetIClear(label->ht[v]));
1661: }
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*@
1666: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1668: Not Collective
1670: Input Parameters:
1671: + label - The `DMLabel`
1672: . value - The label value for all points
1673: . pStart - The first point
1674: - pEnd - A point beyond all marked points
1676: Level: intermediate
1678: Note:
1679: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1681: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1682: @*/
1683: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1684: {
1685: IS pIS;
1687: PetscFunctionBegin;
1688: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1689: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1690: PetscCall(ISDestroy(&pIS));
1691: PetscFunctionReturn(PETSC_SUCCESS);
1692: }
1694: /*@
1695: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1697: Not Collective
1699: Input Parameters:
1700: + label - The `DMLabel`
1701: . value - The label value
1702: - p - A point with this value
1704: Output Parameter:
1705: . 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
1707: Level: intermediate
1709: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1710: @*/
1711: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1712: {
1713: IS pointIS;
1714: PetscInt v;
1716: PetscFunctionBegin;
1718: PetscAssertPointer(index, 4);
1719: *index = -1;
1720: PetscCall(DMLabelLookupStratum(label, value, &v));
1721: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1722: PetscCall(DMLabelMakeValid_Private(label, v));
1723: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1724: PetscCall(ISLocate(pointIS, p, index));
1725: PetscCall(ISDestroy(&pointIS));
1726: PetscFunctionReturn(PETSC_SUCCESS);
1727: }
1729: /*@
1730: DMLabelFilter - Remove all points outside of [`start`, `end`)
1732: Not Collective
1734: Input Parameters:
1735: + label - the `DMLabel`
1736: . start - the first point kept
1737: - end - one more than the last point kept
1739: Level: intermediate
1741: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1742: @*/
1743: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1744: {
1745: PetscInt v;
1747: PetscFunctionBegin;
1749: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1750: PetscCall(DMLabelDestroyIndex(label));
1751: PetscCall(DMLabelMakeAllValid_Private(label));
1752: for (v = 0; v < label->numStrata; ++v) {
1753: PetscCall(ISGeneralFilter(label->points[v], start, end));
1754: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1755: }
1756: PetscCall(DMLabelCreateIndex(label, start, end));
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: /*@
1761: DMLabelPermute - Create a new label with permuted points
1763: Not Collective
1765: Input Parameters:
1766: + label - the `DMLabel`
1767: - permutation - the point permutation
1769: Output Parameter:
1770: . labelNew - the new label containing the permuted points
1772: Level: intermediate
1774: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1775: @*/
1776: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1777: {
1778: const PetscInt *perm;
1779: PetscInt numValues, numPoints, v, q;
1781: PetscFunctionBegin;
1784: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1785: PetscCall(DMLabelMakeAllValid_Private(label));
1786: PetscCall(DMLabelDuplicate(label, labelNew));
1787: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1788: PetscCall(ISGetLocalSize(permutation, &numPoints));
1789: PetscCall(ISGetIndices(permutation, &perm));
1790: for (v = 0; v < numValues; ++v) {
1791: const PetscInt size = (*labelNew)->stratumSizes[v];
1792: const PetscInt *points;
1793: PetscInt *pointsNew;
1795: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1796: PetscCall(PetscCalloc1(size, &pointsNew));
1797: for (q = 0; q < size; ++q) {
1798: const PetscInt point = points[q];
1800: 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);
1801: pointsNew[q] = perm[point];
1802: }
1803: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1804: PetscCall(PetscSortInt(size, pointsNew));
1805: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1806: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1807: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1808: PetscCall(PetscFree(pointsNew));
1809: } else {
1810: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1811: }
1812: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1813: }
1814: PetscCall(ISRestoreIndices(permutation, &perm));
1815: if (label->bt) {
1816: PetscCall(PetscBTDestroy(&label->bt));
1817: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1818: }
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: /*@
1823: DMLabelPermuteValues - Permute the values in a label
1825: Not collective
1827: Input Parameters:
1828: + label - the `DMLabel`
1829: - permutation - the value permutation, permutation[old value] = new value
1831: Output Parameter:
1832: . label - the `DMLabel` now with permuted values
1834: Note:
1835: The modification is done in-place
1837: Level: intermediate
1839: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1840: @*/
1841: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1842: {
1843: PetscInt Nv, Np;
1845: PetscFunctionBegin;
1848: PetscCall(DMLabelGetNumValues(label, &Nv));
1849: PetscCall(ISGetLocalSize(permutation, &Np));
1850: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1851: if (PetscDefined(USE_DEBUG)) {
1852: PetscBool flg;
1853: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1854: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1855: }
1856: PetscCall(DMLabelRewriteValues(label, permutation));
1857: PetscFunctionReturn(PETSC_SUCCESS);
1858: }
1860: /*@
1861: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1863: Not collective
1865: Input Parameters:
1866: + label - the `DMLabel`
1867: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1869: Output Parameter:
1870: . label - the `DMLabel` now with permuted values
1872: Note:
1873: The modification is done in-place
1875: Level: intermediate
1877: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1878: @*/
1879: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1880: {
1881: const PetscInt *perm;
1882: PetscInt Nv, Np;
1884: PetscFunctionBegin;
1887: PetscCall(DMLabelMakeAllValid_Private(label));
1888: PetscCall(DMLabelGetNumValues(label, &Nv));
1889: PetscCall(ISGetLocalSize(permutation, &Np));
1890: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1891: PetscCall(ISGetIndices(permutation, &perm));
1892: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1893: PetscCall(ISRestoreIndices(permutation, &perm));
1894: PetscFunctionReturn(PETSC_SUCCESS);
1895: }
1897: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1898: {
1899: MPI_Comm comm;
1900: PetscInt s, l, nroots, nleaves, offset, size;
1901: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1902: PetscSection rootSection;
1903: PetscSF labelSF;
1905: PetscFunctionBegin;
1906: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1907: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1908: /* Build a section of stratum values per point, generate the according SF
1909: and distribute point-wise stratum values to leaves. */
1910: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1911: PetscCall(PetscSectionCreate(comm, &rootSection));
1912: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1913: if (label) {
1914: for (s = 0; s < label->numStrata; ++s) {
1915: const PetscInt *points;
1917: PetscCall(ISGetIndices(label->points[s], &points));
1918: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1919: PetscCall(ISRestoreIndices(label->points[s], &points));
1920: }
1921: }
1922: PetscCall(PetscSectionSetUp(rootSection));
1923: /* Create a point-wise array of stratum values */
1924: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1925: PetscCall(PetscMalloc1(size, &rootStrata));
1926: PetscCall(PetscCalloc1(nroots, &rootIdx));
1927: if (label) {
1928: for (s = 0; s < label->numStrata; ++s) {
1929: const PetscInt *points;
1931: PetscCall(ISGetIndices(label->points[s], &points));
1932: for (l = 0; l < label->stratumSizes[s]; l++) {
1933: const PetscInt p = points[l];
1934: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1935: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1936: }
1937: PetscCall(ISRestoreIndices(label->points[s], &points));
1938: }
1939: }
1940: /* Build SF that maps label points to remote processes */
1941: PetscCall(PetscSectionCreate(comm, leafSection));
1942: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1943: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1944: PetscCall(PetscFree(remoteOffsets));
1945: /* Send the strata for each point over the derived SF */
1946: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1947: PetscCall(PetscMalloc1(size, leafStrata));
1948: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1949: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1950: /* Clean up */
1951: PetscCall(PetscFree(rootStrata));
1952: PetscCall(PetscFree(rootIdx));
1953: PetscCall(PetscSectionDestroy(&rootSection));
1954: PetscCall(PetscSFDestroy(&labelSF));
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@
1959: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1961: Collective
1963: Input Parameters:
1964: + label - the `DMLabel`
1965: - sf - the map from old to new distribution
1967: Output Parameter:
1968: . labelNew - the new redistributed label
1970: Level: intermediate
1972: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1973: @*/
1974: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1975: {
1976: MPI_Comm comm;
1977: PetscSection leafSection;
1978: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1979: PetscInt *leafStrata, *strataIdx;
1980: PetscInt **points;
1981: const char *lname = NULL;
1982: char *name;
1983: PetscMPIInt nameSize;
1984: PetscHSetI stratumHash;
1985: size_t len = 0;
1986: PetscMPIInt rank;
1988: PetscFunctionBegin;
1990: if (label) {
1992: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1993: PetscCall(DMLabelMakeAllValid_Private(label));
1994: }
1995: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1996: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1997: /* Bcast name */
1998: if (rank == 0) {
1999: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2000: PetscCall(PetscStrlen(lname, &len));
2001: }
2002: PetscCall(PetscMPIIntCast(len, &nameSize));
2003: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2004: PetscCall(PetscMalloc1(nameSize + 1, &name));
2005: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2006: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2007: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2008: PetscCall(PetscFree(name));
2009: /* Bcast defaultValue */
2010: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
2011: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
2012: /* Distribute stratum values over the SF and get the point mapping on the receiver */
2013: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
2014: /* Determine received stratum values and initialise new label*/
2015: PetscCall(PetscHSetICreate(&stratumHash));
2016: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
2017: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
2018: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
2019: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
2020: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
2021: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
2022: /* Turn leafStrata into indices rather than stratum values */
2023: offset = 0;
2024: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
2025: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
2026: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
2027: for (p = 0; p < size; ++p) {
2028: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2029: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
2030: leafStrata[p] = s;
2031: break;
2032: }
2033: }
2034: }
2035: /* Rebuild the point strata on the receiver */
2036: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
2037: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2038: for (p = pStart; p < pEnd; p++) {
2039: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2040: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2041: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
2042: }
2043: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
2044: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
2045: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
2046: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2047: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
2048: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
2049: }
2050: /* Insert points into new strata */
2051: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
2052: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2053: for (p = pStart; p < pEnd; p++) {
2054: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2055: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2056: for (s = 0; s < dof; s++) {
2057: stratum = leafStrata[offset + s];
2058: points[stratum][strataIdx[stratum]++] = p;
2059: }
2060: }
2061: for (s = 0; s < (*labelNew)->numStrata; s++) {
2062: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
2063: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
2064: }
2065: PetscCall(PetscFree(points));
2066: PetscCall(PetscHSetIDestroy(&stratumHash));
2067: PetscCall(PetscFree(leafStrata));
2068: PetscCall(PetscFree(strataIdx));
2069: PetscCall(PetscSectionDestroy(&leafSection));
2070: PetscFunctionReturn(PETSC_SUCCESS);
2071: }
2073: /*@
2074: DMLabelGather - Gather all label values from leafs into roots
2076: Collective
2078: Input Parameters:
2079: + label - the `DMLabel`
2080: - sf - the `PetscSF` communication map
2082: Output Parameter:
2083: . labelNew - the new `DMLabel` with localised leaf values
2085: Level: developer
2087: Note:
2088: This is the inverse operation to `DMLabelDistribute()`.
2090: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2091: @*/
2092: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
2093: {
2094: MPI_Comm comm;
2095: PetscSection rootSection;
2096: PetscSF sfLabel;
2097: PetscSFNode *rootPoints, *leafPoints;
2098: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2099: const PetscInt *rootDegree, *ilocal;
2100: PetscInt *rootStrata;
2101: const char *lname;
2102: char *name;
2103: PetscMPIInt nameSize;
2104: size_t len = 0;
2105: PetscMPIInt rank, size;
2107: PetscFunctionBegin;
2110: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2111: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2112: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2113: PetscCallMPI(MPI_Comm_size(comm, &size));
2114: /* Bcast name */
2115: if (rank == 0) {
2116: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2117: PetscCall(PetscStrlen(lname, &len));
2118: }
2119: PetscCall(PetscMPIIntCast(len, &nameSize));
2120: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2121: PetscCall(PetscMalloc1(nameSize + 1, &name));
2122: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2123: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2124: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2125: PetscCall(PetscFree(name));
2126: /* Gather rank/index pairs of leaves into local roots to build
2127: an inverse, multi-rooted SF. Note that this ignores local leaf
2128: indexing due to the use of the multiSF in PetscSFGather. */
2129: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2130: PetscCall(PetscMalloc1(nroots, &leafPoints));
2131: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2132: for (p = 0; p < nleaves; p++) {
2133: PetscInt ilp = ilocal ? ilocal[p] : p;
2135: leafPoints[ilp].index = ilp;
2136: leafPoints[ilp].rank = rank;
2137: }
2138: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2139: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2140: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2141: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2142: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2143: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2144: PetscCall(PetscSFCreate(comm, &sfLabel));
2145: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2146: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2147: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2148: /* Rebuild the point strata on the receiver */
2149: for (p = 0, idx = 0; p < nroots; p++) {
2150: for (d = 0; d < rootDegree[p]; d++) {
2151: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2152: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2153: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2154: }
2155: idx += rootDegree[p];
2156: }
2157: PetscCall(PetscFree(leafPoints));
2158: PetscCall(PetscFree(rootStrata));
2159: PetscCall(PetscSectionDestroy(&rootSection));
2160: PetscCall(PetscSFDestroy(&sfLabel));
2161: PetscFunctionReturn(PETSC_SUCCESS);
2162: }
2164: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2165: {
2166: const PetscInt *degree;
2167: const PetscInt *points;
2168: PetscInt Nr, r, Nl, l, val, defVal;
2170: PetscFunctionBegin;
2171: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2172: /* Add in leaves */
2173: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2174: for (l = 0; l < Nl; ++l) {
2175: PetscCall(DMLabelGetValue(label, points[l], &val));
2176: if (val != defVal) valArray[points[l]] = val;
2177: }
2178: /* Add in shared roots */
2179: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2180: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2181: for (r = 0; r < Nr; ++r) {
2182: if (degree[r]) {
2183: PetscCall(DMLabelGetValue(label, r, &val));
2184: if (val != defVal) valArray[r] = val;
2185: }
2186: }
2187: PetscFunctionReturn(PETSC_SUCCESS);
2188: }
2190: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), PetscCtx ctx)
2191: {
2192: const PetscInt *degree;
2193: const PetscInt *points;
2194: PetscInt Nr, r, Nl, l, val, defVal;
2196: PetscFunctionBegin;
2197: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2198: /* Read out leaves */
2199: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2200: for (l = 0; l < Nl; ++l) {
2201: const PetscInt p = points[l];
2202: const PetscInt cval = valArray[p];
2204: if (cval != defVal) {
2205: PetscCall(DMLabelGetValue(label, p, &val));
2206: if (val == defVal) {
2207: PetscCall(DMLabelSetValue(label, p, cval));
2208: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2209: }
2210: }
2211: }
2212: /* Read out shared roots */
2213: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2214: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2215: for (r = 0; r < Nr; ++r) {
2216: if (degree[r]) {
2217: const PetscInt cval = valArray[r];
2219: if (cval != defVal) {
2220: PetscCall(DMLabelGetValue(label, r, &val));
2221: if (val == defVal) {
2222: PetscCall(DMLabelSetValue(label, r, cval));
2223: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2224: }
2225: }
2226: }
2227: }
2228: PetscFunctionReturn(PETSC_SUCCESS);
2229: }
2231: /*@
2232: DMLabelPropagateBegin - Setup a cycle of label propagation
2234: Collective
2236: Input Parameters:
2237: + label - The `DMLabel` to propagate across processes
2238: - sf - The `PetscSF` describing parallel layout of the label points
2240: Level: intermediate
2242: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2243: @*/
2244: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2245: {
2246: PetscInt Nr, r, defVal;
2247: PetscMPIInt size;
2249: PetscFunctionBegin;
2250: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2251: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2252: if (size > 1) {
2253: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2254: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2255: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2256: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2257: }
2258: PetscFunctionReturn(PETSC_SUCCESS);
2259: }
2261: /*@
2262: DMLabelPropagateEnd - Tear down a cycle of label propagation
2264: Collective
2266: Input Parameters:
2267: + label - The `DMLabel` to propagate across processes
2268: - pointSF - The `PetscSF` describing parallel layout of the label points
2270: Level: intermediate
2272: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2273: @*/
2274: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2275: {
2276: PetscFunctionBegin;
2277: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2278: PetscCall(PetscFree(label->propArray));
2279: label->propArray = NULL;
2280: PetscFunctionReturn(PETSC_SUCCESS);
2281: }
2283: /*@C
2284: DMLabelPropagatePush - Tear down a cycle of label propagation
2286: Collective
2288: Input Parameters:
2289: + label - The `DMLabel` to propagate across processes
2290: . pointSF - The `PetscSF` describing parallel layout of the label points
2291: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2292: - ctx - An optional user context for the callback, or `NULL`
2294: Calling sequence of `markPoint`:
2295: + label - The `DMLabel`
2296: . p - The point being marked
2297: . val - The label value for `p`
2298: - ctx - An optional user context
2300: Level: intermediate
2302: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2303: @*/
2304: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, PetscCtx ctx), PetscCtx ctx)
2305: {
2306: PetscInt *valArray = label->propArray, Nr;
2307: PetscMPIInt size;
2309: PetscFunctionBegin;
2310: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2311: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2312: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2313: if (size > 1 && Nr >= 0) {
2314: /* Communicate marked edges
2315: The current implementation allocates an array the size of the number of root. We put the label values into the
2316: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2318: TODO: We could use in-place communication with a different SF
2319: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2320: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2322: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2323: 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
2324: edge to the queue.
2325: */
2326: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2327: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2328: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2329: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2330: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2331: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2332: }
2333: PetscFunctionReturn(PETSC_SUCCESS);
2334: }
2336: /*@
2337: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2339: Not Collective
2341: Input Parameter:
2342: . label - the `DMLabel`
2344: Output Parameters:
2345: + section - the section giving offsets for each stratum
2346: - is - An `IS` containing all the label points
2348: Level: developer
2350: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2351: @*/
2352: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2353: {
2354: IS vIS;
2355: const PetscInt *values;
2356: PetscInt *points;
2357: PetscInt nV, vS = 0, vE = 0, v, N;
2359: PetscFunctionBegin;
2361: PetscCall(DMLabelGetNumValues(label, &nV));
2362: PetscCall(DMLabelGetValueIS(label, &vIS));
2363: PetscCall(ISGetIndices(vIS, &values));
2364: if (nV) {
2365: vS = values[0];
2366: vE = values[0] + 1;
2367: }
2368: for (v = 1; v < nV; ++v) {
2369: vS = PetscMin(vS, values[v]);
2370: vE = PetscMax(vE, values[v] + 1);
2371: }
2372: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2373: PetscCall(PetscSectionSetChart(*section, vS, vE));
2374: for (v = 0; v < nV; ++v) {
2375: PetscInt n;
2377: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2378: PetscCall(PetscSectionSetDof(*section, values[v], n));
2379: }
2380: PetscCall(PetscSectionSetUp(*section));
2381: PetscCall(PetscSectionGetStorageSize(*section, &N));
2382: PetscCall(PetscMalloc1(N, &points));
2383: for (v = 0; v < nV; ++v) {
2384: IS is;
2385: const PetscInt *spoints;
2386: PetscInt dof, off, p;
2388: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2389: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2390: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2391: PetscCall(ISGetIndices(is, &spoints));
2392: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2393: PetscCall(ISRestoreIndices(is, &spoints));
2394: PetscCall(ISDestroy(&is));
2395: }
2396: PetscCall(ISRestoreIndices(vIS, &values));
2397: PetscCall(ISDestroy(&vIS));
2398: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2399: PetscFunctionReturn(PETSC_SUCCESS);
2400: }
2402: /*@C
2403: DMLabelRegister - Adds a new label component implementation
2405: Not Collective
2407: Input Parameters:
2408: + name - The name of a new user-defined creation routine
2409: - create_func - The creation routine itself
2411: Notes:
2412: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2414: Example Usage:
2415: .vb
2416: DMLabelRegister("my_label", MyLabelCreate);
2417: .ve
2419: Then, your label type can be chosen with the procedural interface via
2420: .vb
2421: DMLabelCreate(MPI_Comm, DMLabel *);
2422: DMLabelSetType(DMLabel, "my_label");
2423: .ve
2424: or at runtime via the option
2425: .vb
2426: -dm_label_type my_label
2427: .ve
2429: Level: advanced
2431: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2432: @*/
2433: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2434: {
2435: PetscFunctionBegin;
2436: PetscCall(DMInitializePackage());
2437: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2438: PetscFunctionReturn(PETSC_SUCCESS);
2439: }
2441: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2442: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2444: /*@C
2445: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2447: Not Collective
2449: Level: advanced
2451: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2452: @*/
2453: PetscErrorCode DMLabelRegisterAll(void)
2454: {
2455: PetscFunctionBegin;
2456: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2457: DMLabelRegisterAllCalled = PETSC_TRUE;
2459: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2460: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2461: PetscFunctionReturn(PETSC_SUCCESS);
2462: }
2464: /*@C
2465: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2467: Level: developer
2469: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2470: @*/
2471: PetscErrorCode DMLabelRegisterDestroy(void)
2472: {
2473: PetscFunctionBegin;
2474: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2475: DMLabelRegisterAllCalled = PETSC_FALSE;
2476: PetscFunctionReturn(PETSC_SUCCESS);
2477: }
2479: /*@
2480: DMLabelSetType - Sets the particular implementation for a label.
2482: Collective
2484: Input Parameters:
2485: + label - The label
2486: - method - The name of the label type
2488: Options Database Key:
2489: . -dm_label_type type - Sets the label type; see `DMLabelType`
2491: Level: intermediate
2493: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`, `DMLabelType`
2494: @*/
2495: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2496: {
2497: PetscErrorCode (*r)(DMLabel);
2498: PetscBool match;
2500: PetscFunctionBegin;
2502: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2503: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2505: PetscCall(DMLabelRegisterAll());
2506: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2507: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2509: PetscTryTypeMethod(label, destroy);
2510: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2511: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2512: PetscCall((*r)(label));
2513: PetscFunctionReturn(PETSC_SUCCESS);
2514: }
2516: /*@
2517: DMLabelGetType - Gets the type name (as a string) from the label.
2519: Not Collective
2521: Input Parameter:
2522: . label - The `DMLabel`
2524: Output Parameter:
2525: . type - The `DMLabel` type name
2527: Level: intermediate
2529: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2530: @*/
2531: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2532: {
2533: PetscFunctionBegin;
2535: PetscAssertPointer(type, 2);
2536: PetscCall(DMLabelRegisterAll());
2537: *type = ((PetscObject)label)->type_name;
2538: PetscFunctionReturn(PETSC_SUCCESS);
2539: }
2541: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2542: {
2543: PetscFunctionBegin;
2544: label->ops->view = DMLabelView_Concrete;
2545: label->ops->setup = NULL;
2546: label->ops->duplicate = DMLabelDuplicate_Concrete;
2547: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2548: PetscFunctionReturn(PETSC_SUCCESS);
2549: }
2551: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2552: {
2553: PetscFunctionBegin;
2555: PetscCall(DMLabelInitialize_Concrete(label));
2556: PetscFunctionReturn(PETSC_SUCCESS);
2557: }
2559: /*@
2560: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2561: the local section and an `PetscSF` describing the section point overlap.
2563: Collective
2565: Input Parameters:
2566: + s - The `PetscSection` for the local field layout
2567: . sf - The `PetscSF` describing parallel layout of the section points
2568: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2569: . label - The label specifying the points
2570: - labelValue - The label stratum specifying the points
2572: Output Parameter:
2573: . gsection - The `PetscSection` for the global field layout
2575: Level: developer
2577: Note:
2578: This gives negative sizes and offsets to points not owned by this process
2580: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2581: @*/
2582: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2583: {
2584: PetscInt *neg = NULL, *tmpOff = NULL;
2585: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2587: PetscFunctionBegin;
2591: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2592: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2593: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2594: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2595: if (nroots >= 0) {
2596: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2597: PetscCall(PetscCalloc1(nroots, &neg));
2598: if (nroots > pEnd - pStart) {
2599: PetscCall(PetscCalloc1(nroots, &tmpOff));
2600: } else {
2601: tmpOff = &(*gsection)->atlasDof[-pStart];
2602: }
2603: }
2604: /* Mark ghost points with negative dof */
2605: for (p = pStart; p < pEnd; ++p) {
2606: PetscInt value;
2608: PetscCall(DMLabelGetValue(label, p, &value));
2609: if (value != labelValue) continue;
2610: PetscCall(PetscSectionGetDof(s, p, &dof));
2611: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2612: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2613: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2614: if (neg) neg[p] = -(dof + 1);
2615: }
2616: PetscCall(PetscSectionSetUpBC(*gsection));
2617: if (nroots >= 0) {
2618: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2619: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2620: if (nroots > pEnd - pStart) {
2621: for (p = pStart; p < pEnd; ++p) {
2622: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2623: }
2624: }
2625: }
2626: /* Calculate new sizes, get process offset, and calculate point offsets */
2627: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2628: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2629: (*gsection)->atlasOff[p] = off;
2630: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2631: }
2632: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2633: globalOff -= off;
2634: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2635: (*gsection)->atlasOff[p] += globalOff;
2636: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2637: }
2638: /* Put in negative offsets for ghost points */
2639: if (nroots >= 0) {
2640: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2641: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2642: if (nroots > pEnd - pStart) {
2643: for (p = pStart; p < pEnd; ++p) {
2644: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2645: }
2646: }
2647: }
2648: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2649: PetscCall(PetscFree(neg));
2650: PetscFunctionReturn(PETSC_SUCCESS);
2651: }
2653: typedef struct _n_PetscSectionSym_Label {
2654: DMLabel label;
2655: PetscCopyMode *modes;
2656: PetscInt *sizes;
2657: const PetscInt ***perms;
2658: const PetscScalar ***rots;
2659: PetscInt (*minMaxOrients)[2];
2660: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2661: } PetscSectionSym_Label;
2663: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2664: {
2665: PetscInt i, j;
2666: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2668: PetscFunctionBegin;
2669: for (i = 0; i <= sl->numStrata; i++) {
2670: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2671: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2672: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2673: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2674: }
2675: if (sl->perms[i]) {
2676: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2678: PetscCall(PetscFree(perms));
2679: }
2680: if (sl->rots[i]) {
2681: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2683: PetscCall(PetscFree(rots));
2684: }
2685: }
2686: }
2687: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2688: PetscCall(DMLabelDestroy(&sl->label));
2689: sl->numStrata = 0;
2690: PetscFunctionReturn(PETSC_SUCCESS);
2691: }
2693: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2694: {
2695: PetscFunctionBegin;
2696: PetscCall(PetscSectionSymLabelReset(sym));
2697: PetscCall(PetscFree(sym->data));
2698: PetscFunctionReturn(PETSC_SUCCESS);
2699: }
2701: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2702: {
2703: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2704: PetscBool isAscii;
2705: DMLabel label = sl->label;
2706: const char *name;
2708: PetscFunctionBegin;
2709: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2710: if (isAscii) {
2711: PetscInt i, j, k;
2712: PetscViewerFormat format;
2714: PetscCall(PetscViewerGetFormat(viewer, &format));
2715: if (label) {
2716: PetscCall(PetscViewerGetFormat(viewer, &format));
2717: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2718: PetscCall(PetscViewerASCIIPushTab(viewer));
2719: PetscCall(DMLabelView(label, viewer));
2720: PetscCall(PetscViewerASCIIPopTab(viewer));
2721: } else {
2722: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2723: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2724: }
2725: } else {
2726: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2727: }
2728: PetscCall(PetscViewerASCIIPushTab(viewer));
2729: for (i = 0; i <= sl->numStrata; i++) {
2730: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2732: if (!(sl->perms[i] || sl->rots[i])) {
2733: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2734: } else {
2735: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2736: PetscCall(PetscViewerASCIIPushTab(viewer));
2737: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2738: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2739: PetscCall(PetscViewerASCIIPushTab(viewer));
2740: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2741: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2742: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2743: } else {
2744: PetscInt tab;
2746: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2747: PetscCall(PetscViewerASCIIPushTab(viewer));
2748: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2749: if (sl->perms[i] && sl->perms[i][j]) {
2750: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2751: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2752: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2753: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2754: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2755: }
2756: if (sl->rots[i] && sl->rots[i][j]) {
2757: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2758: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2759: #if defined(PETSC_USE_COMPLEX)
2760: 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])));
2761: #else
2762: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2763: #endif
2764: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2765: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2766: }
2767: PetscCall(PetscViewerASCIIPopTab(viewer));
2768: }
2769: }
2770: PetscCall(PetscViewerASCIIPopTab(viewer));
2771: }
2772: PetscCall(PetscViewerASCIIPopTab(viewer));
2773: }
2774: }
2775: PetscCall(PetscViewerASCIIPopTab(viewer));
2776: }
2777: PetscFunctionReturn(PETSC_SUCCESS);
2778: }
2780: /*@
2781: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2783: Logically
2785: Input Parameters:
2786: + sym - the section symmetries
2787: - label - the `DMLabel` describing the types of points
2789: Level: developer:
2791: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2792: @*/
2793: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2794: {
2795: PetscSectionSym_Label *sl;
2797: PetscFunctionBegin;
2799: sl = (PetscSectionSym_Label *)sym->data;
2800: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2801: if (label) {
2802: sl->label = label;
2803: PetscCall(PetscObjectReference((PetscObject)label));
2804: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2805: 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));
2806: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2807: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2808: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2809: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2810: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2811: }
2812: PetscFunctionReturn(PETSC_SUCCESS);
2813: }
2815: /*@C
2816: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2818: Logically Collective
2820: Input Parameters:
2821: + sym - the section symmetries
2822: - stratum - the stratum value in the label that we are assigning symmetries for
2824: Output Parameters:
2825: + size - the number of dofs for points in the `stratum` of the label
2826: . minOrient - the smallest orientation for a point in this `stratum`
2827: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2828: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2829: - 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
2831: Level: developer
2833: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2834: @*/
2835: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2836: {
2837: PetscSectionSym_Label *sl;
2838: const char *name;
2839: PetscInt i;
2841: PetscFunctionBegin;
2843: sl = (PetscSectionSym_Label *)sym->data;
2844: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2845: for (i = 0; i <= sl->numStrata; i++) {
2846: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2848: if (stratum == value) break;
2849: }
2850: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2851: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2852: if (size) {
2853: PetscAssertPointer(size, 3);
2854: *size = sl->sizes[i];
2855: }
2856: if (minOrient) {
2857: PetscAssertPointer(minOrient, 4);
2858: *minOrient = sl->minMaxOrients[i][0];
2859: }
2860: if (maxOrient) {
2861: PetscAssertPointer(maxOrient, 5);
2862: *maxOrient = sl->minMaxOrients[i][1];
2863: }
2864: if (perms) {
2865: PetscAssertPointer(perms, 6);
2866: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2867: }
2868: if (rots) {
2869: PetscAssertPointer(rots, 7);
2870: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2871: }
2872: PetscFunctionReturn(PETSC_SUCCESS);
2873: }
2875: /*@C
2876: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2878: Logically
2880: Input Parameters:
2881: + sym - the section symmetries
2882: . stratum - the stratum value in the label that we are assigning symmetries for
2883: . size - the number of dofs for points in the `stratum` of the label
2884: . minOrient - the smallest orientation for a point in this `stratum`
2885: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2886: . mode - how `sym` should copy the `perms` and `rots` arrays
2887: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2888: - 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
2890: Level: developer
2892: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2893: @*/
2894: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2895: {
2896: PetscSectionSym_Label *sl;
2897: const char *name;
2898: PetscInt i, j, k;
2900: PetscFunctionBegin;
2902: sl = (PetscSectionSym_Label *)sym->data;
2903: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2904: for (i = 0; i <= sl->numStrata; i++) {
2905: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2907: if (stratum == value) break;
2908: }
2909: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2910: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2911: sl->sizes[i] = size;
2912: sl->modes[i] = mode;
2913: sl->minMaxOrients[i][0] = minOrient;
2914: sl->minMaxOrients[i][1] = maxOrient;
2915: if (mode == PETSC_COPY_VALUES) {
2916: if (perms) {
2917: PetscInt **ownPerms;
2919: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2920: for (j = 0; j < maxOrient - minOrient; j++) {
2921: if (perms[j]) {
2922: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2923: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2924: }
2925: }
2926: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2927: }
2928: if (rots) {
2929: PetscScalar **ownRots;
2931: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2932: for (j = 0; j < maxOrient - minOrient; j++) {
2933: if (rots[j]) {
2934: PetscCall(PetscMalloc1(size, &ownRots[j]));
2935: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2936: }
2937: }
2938: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2939: }
2940: } else {
2941: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2942: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2943: }
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2948: {
2949: PetscInt i, j, numStrata;
2950: PetscSectionSym_Label *sl;
2951: DMLabel label;
2953: PetscFunctionBegin;
2954: sl = (PetscSectionSym_Label *)sym->data;
2955: numStrata = sl->numStrata;
2956: label = sl->label;
2957: for (i = 0; i < numPoints; i++) {
2958: PetscInt point = points[2 * i];
2959: PetscInt ornt = points[2 * i + 1];
2961: for (j = 0; j < numStrata; j++) {
2962: if (label->validIS[j]) {
2963: PetscInt k;
2965: PetscCall(ISLocate(label->points[j], point, &k));
2966: if (k >= 0) break;
2967: } else {
2968: PetscBool has;
2970: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2971: if (has) break;
2972: }
2973: }
2974: 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],
2975: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2976: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2977: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2978: }
2979: PetscFunctionReturn(PETSC_SUCCESS);
2980: }
2982: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2983: {
2984: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2985: IS valIS;
2986: const PetscInt *values;
2987: PetscInt Nv;
2989: PetscFunctionBegin;
2990: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2991: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2992: PetscCall(ISGetIndices(valIS, &values));
2993: for (PetscInt v = 0; v < Nv; ++v) {
2994: const PetscInt val = values[v];
2995: PetscInt size, minOrient, maxOrient;
2996: const PetscInt **perms;
2997: const PetscScalar **rots;
2999: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
3000: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
3001: }
3002: PetscCall(ISDestroy(&valIS));
3003: PetscFunctionReturn(PETSC_SUCCESS);
3004: }
3006: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3007: {
3008: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
3009: DMLabel dlabel;
3011: PetscFunctionBegin;
3012: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
3013: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
3014: PetscCall(DMLabelDestroy(&dlabel));
3015: PetscCall(PetscSectionSymCopy(sym, *dsym));
3016: PetscFunctionReturn(PETSC_SUCCESS);
3017: }
3019: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
3020: {
3021: PetscSectionSym_Label *sl;
3023: PetscFunctionBegin;
3024: PetscCall(PetscNew(&sl));
3025: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
3026: sym->ops->distribute = PetscSectionSymDistribute_Label;
3027: sym->ops->copy = PetscSectionSymCopy_Label;
3028: sym->ops->view = PetscSectionSymView_Label;
3029: sym->ops->destroy = PetscSectionSymDestroy_Label;
3030: sym->data = (void *)sl;
3031: PetscFunctionReturn(PETSC_SUCCESS);
3032: }
3034: /*@
3035: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
3037: Collective
3039: Input Parameters:
3040: + comm - the MPI communicator for the new symmetry
3041: - label - the label defining the strata
3043: Output Parameter:
3044: . sym - the section symmetries
3046: Level: developer
3048: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3049: @*/
3050: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
3051: {
3052: PetscFunctionBegin;
3053: PetscCall(DMInitializePackage());
3054: PetscCall(PetscSectionSymCreate(comm, sym));
3055: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
3056: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
3057: PetscFunctionReturn(PETSC_SUCCESS);
3058: }