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 viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
492: Level: intermediate
494: Note:
495: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
496: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
498: .seealso: [](ch_dmbase), `DMLabel`, `DMLabelView()`, `PetscObjectViewFromOptions()`, `DMLabelCreate()`, `PetscOptionsCreateViewer()`
499: @*/
500: PetscErrorCode DMLabelViewFromOptions(DMLabel label, PeOp PetscObject obj, const char name[])
501: {
502: PetscFunctionBegin;
504: PetscCall(PetscObjectViewFromOptions((PetscObject)label, obj, name));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: /*@
509: DMLabelReset - Destroys internal data structures in a `DMLabel`
511: Not Collective
513: Input Parameter:
514: . label - The `DMLabel`
516: Level: beginner
518: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
519: @*/
520: PetscErrorCode DMLabelReset(DMLabel label)
521: {
522: PetscFunctionBegin;
524: for (PetscInt v = 0; v < label->numStrata; ++v) {
525: PetscCall(PetscHSetIDestroy(&label->ht[v]));
526: PetscCall(ISDestroy(&label->points[v]));
527: }
528: label->numStrata = 0;
529: PetscCall(PetscFree(label->stratumValues));
530: PetscCall(PetscFree(label->stratumSizes));
531: PetscCall(PetscFree(label->ht));
532: PetscCall(PetscFree(label->points));
533: PetscCall(PetscFree(label->validIS));
534: PetscCall(PetscHMapIReset(label->hmap));
535: label->pStart = -1;
536: label->pEnd = -1;
537: PetscCall(PetscBTDestroy(&label->bt));
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: DMLabelDestroy - Destroys a `DMLabel`
544: Collective
546: Input Parameter:
547: . label - The `DMLabel`
549: Level: beginner
551: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
552: @*/
553: PetscErrorCode DMLabelDestroy(DMLabel *label)
554: {
555: PetscFunctionBegin;
556: if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
558: if (--((PetscObject)*label)->refct > 0) {
559: *label = NULL;
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
562: PetscCall(DMLabelReset(*label));
563: PetscCall(PetscHMapIDestroy(&(*label)->hmap));
564: PetscCall(PetscHeaderDestroy(label));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
569: {
570: PetscFunctionBegin;
571: for (PetscInt v = 0; v < label->numStrata; ++v) {
572: PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
573: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
574: (*labelnew)->points[v] = label->points[v];
575: }
576: PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
577: PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@
582: DMLabelDuplicate - Duplicates a `DMLabel`
584: Collective
586: Input Parameter:
587: . label - The `DMLabel`
589: Output Parameter:
590: . labelnew - new label
592: Level: intermediate
594: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
595: @*/
596: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
597: {
598: const char *name;
600: PetscFunctionBegin;
602: PetscCall(DMLabelMakeAllValid_Private(label));
603: PetscCall(PetscObjectGetName((PetscObject)label, &name));
604: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));
606: (*labelnew)->numStrata = label->numStrata;
607: (*labelnew)->defaultValue = label->defaultValue;
608: (*labelnew)->readonly = label->readonly;
609: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
610: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
611: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
612: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
613: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
614: for (PetscInt v = 0; v < label->numStrata; ++v) {
615: (*labelnew)->stratumValues[v] = label->stratumValues[v];
616: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
617: (*labelnew)->validIS[v] = PETSC_TRUE;
618: }
619: (*labelnew)->pStart = -1;
620: (*labelnew)->pEnd = -1;
621: (*labelnew)->bt = NULL;
622: PetscUseTypeMethod(label, duplicate, labelnew);
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: /*@
627: DMLabelCompare - Compare two `DMLabel` objects
629: Collective; No Fortran Support
631: Input Parameters:
632: + comm - Comm over which to compare labels
633: . l0 - First `DMLabel`
634: - l1 - Second `DMLabel`
636: Output Parameters:
637: + equal - (Optional) Flag whether the two labels are equal
638: - message - (Optional) Message describing the difference
640: Level: intermediate
642: Notes:
643: The output flag equal is the same on all processes.
644: If it is passed as `NULL` and difference is found, an error is thrown on all processes.
645: Make sure to pass `NULL` on all processes.
647: The output message is set independently on each rank.
648: It is set to `NULL` if no difference was found on the current rank. It must be freed by user.
649: If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner.
650: Make sure to pass `NULL` on all processes.
652: For the comparison, we ignore the order of stratum values, and strata with no points.
654: The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel.
656: Developer Note:
657: Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()`
659: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
660: @*/
661: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
662: {
663: const char *name0, *name1;
664: char msg[PETSC_MAX_PATH_LEN] = "";
665: PetscBool eq;
666: PetscMPIInt rank;
668: PetscFunctionBegin;
671: if (equal) PetscAssertPointer(equal, 4);
672: if (message) PetscAssertPointer(message, 5);
673: PetscCallMPI(MPI_Comm_rank(comm, &rank));
674: PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
675: PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
676: {
677: PetscInt v0, v1;
679: PetscCall(DMLabelGetDefaultValue(l0, &v0));
680: PetscCall(DMLabelGetDefaultValue(l1, &v1));
681: eq = (PetscBool)(v0 == v1);
682: 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));
683: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
684: if (!eq) goto finish;
685: }
686: {
687: IS is0, is1;
689: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
690: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
691: PetscCall(ISEqual(is0, is1, &eq));
692: PetscCall(ISDestroy(&is0));
693: PetscCall(ISDestroy(&is1));
694: if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
695: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
696: if (!eq) goto finish;
697: }
698: {
699: PetscInt nValues;
701: PetscCall(DMLabelGetNumValues(l0, &nValues));
702: for (PetscInt i = 0; i < nValues; i++) {
703: const PetscInt v = l0->stratumValues[i];
704: PetscInt n;
705: IS is0, is1;
707: PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
708: if (!n) continue;
709: PetscCall(DMLabelGetStratumIS(l0, v, &is0));
710: PetscCall(DMLabelGetStratumIS(l1, v, &is1));
711: PetscCall(ISEqualUnsorted(is0, is1, &eq));
712: PetscCall(ISDestroy(&is0));
713: PetscCall(ISDestroy(&is1));
714: if (!eq) {
715: 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));
716: break;
717: }
718: }
719: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPI_C_BOOL, MPI_LAND, comm));
720: }
721: finish:
722: /* If message output arg not set, print to stderr */
723: if (message) {
724: *message = NULL;
725: if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
726: } else {
727: if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
728: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
729: }
730: /* If same output arg not ser and labels are not equal, throw error */
731: if (equal) *equal = eq;
732: else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*@
737: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
739: Not Collective
741: Input Parameter:
742: . label - The `DMLabel`
744: Level: intermediate
746: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
747: @*/
748: PetscErrorCode DMLabelComputeIndex(DMLabel label)
749: {
750: PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v;
752: PetscFunctionBegin;
754: PetscCall(DMLabelMakeAllValid_Private(label));
755: for (v = 0; v < label->numStrata; ++v) {
756: const PetscInt *points;
758: PetscCall(ISGetIndices(label->points[v], &points));
759: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
760: const PetscInt point = points[i];
762: pStart = PetscMin(point, pStart);
763: pEnd = PetscMax(point + 1, pEnd);
764: }
765: PetscCall(ISRestoreIndices(label->points[v], &points));
766: }
767: label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
768: label->pEnd = pEnd;
769: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@
774: DMLabelCreateIndex - Create an index structure for membership determination
776: Not Collective
778: Input Parameters:
779: + label - The `DMLabel`
780: . pStart - The smallest point
781: - pEnd - The largest point + 1
783: Level: intermediate
785: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
786: @*/
787: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
788: {
789: PetscFunctionBegin;
791: PetscCall(DMLabelDestroyIndex(label));
792: PetscCall(DMLabelMakeAllValid_Private(label));
793: label->pStart = pStart;
794: label->pEnd = pEnd;
795: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
796: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
797: for (PetscInt v = 0; v < label->numStrata; ++v) {
798: IS pointIS;
799: const PetscInt *points;
801: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
802: PetscCall(ISGetIndices(pointIS, &points));
803: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
804: const PetscInt point = points[i];
806: 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);
807: PetscCall(PetscBTSet(label->bt, point - pStart));
808: }
809: PetscCall(ISRestoreIndices(label->points[v], &points));
810: PetscCall(ISDestroy(&pointIS));
811: }
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /*@
816: DMLabelDestroyIndex - Destroy the index structure
818: Not Collective
820: Input Parameter:
821: . label - the `DMLabel`
823: Level: intermediate
825: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
826: @*/
827: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
828: {
829: PetscFunctionBegin;
831: label->pStart = -1;
832: label->pEnd = -1;
833: PetscCall(PetscBTDestroy(&label->bt));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*@
838: DMLabelGetBounds - Return the smallest and largest point in the label
840: Not Collective
842: Input Parameter:
843: . label - the `DMLabel`
845: Output Parameters:
846: + pStart - The smallest point
847: - pEnd - The largest point + 1
849: Level: intermediate
851: Note:
852: This will compute an index for the label if one does not exist.
854: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
855: @*/
856: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
857: {
858: PetscFunctionBegin;
860: if (label->pStart == -1 && label->pEnd == -1) PetscCall(DMLabelComputeIndex(label));
861: if (pStart) {
862: PetscAssertPointer(pStart, 2);
863: *pStart = label->pStart;
864: }
865: if (pEnd) {
866: PetscAssertPointer(pEnd, 3);
867: *pEnd = label->pEnd;
868: }
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: /*@
873: DMLabelHasValue - Determine whether a label assigns the value to any point
875: Not Collective
877: Input Parameters:
878: + label - the `DMLabel`
879: - value - the value
881: Output Parameter:
882: . contains - Flag indicating whether the label maps this value to any point
884: Level: developer
886: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
887: @*/
888: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
889: {
890: PetscInt v;
892: PetscFunctionBegin;
894: PetscAssertPointer(contains, 3);
895: PetscCall(DMLabelLookupStratum(label, value, &v));
896: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: /*@
901: DMLabelHasPoint - Determine whether a label assigns a value to a point
903: Not Collective
905: Input Parameters:
906: + label - the `DMLabel`
907: - point - the point
909: Output Parameter:
910: . contains - Flag indicating whether the label maps this point to a value
912: Level: developer
914: Note:
915: The user must call `DMLabelCreateIndex()` before this function.
917: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
918: @*/
919: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
920: {
921: PetscInt pStart, pEnd;
923: PetscFunctionBeginHot;
925: PetscAssertPointer(contains, 3);
926: /* DMLabelGetBounds() calls DMLabelCreateIndex() only if needed */
927: PetscCall(DMLabelGetBounds(label, &pStart, &pEnd));
928: PetscCall(DMLabelMakeAllValid_Private(label));
929: *contains = point >= pStart && point < pEnd && (PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE);
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: /*@
934: DMLabelStratumHasPoint - Return true if the stratum contains a point
936: Not Collective
938: Input Parameters:
939: + label - the `DMLabel`
940: . value - the stratum value
941: - point - the point
943: Output Parameter:
944: . contains - true if the stratum contains the point
946: Level: intermediate
948: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
949: @*/
950: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
951: {
952: PetscFunctionBeginHot;
954: PetscAssertPointer(contains, 4);
955: if (value == label->defaultValue) {
956: PetscInt pointVal;
958: PetscCall(DMLabelGetValue(label, point, &pointVal));
959: *contains = (PetscBool)(pointVal == value);
960: } else {
961: PetscInt v;
963: PetscCall(DMLabelLookupStratum(label, value, &v));
964: if (v >= 0) {
965: if (label->validIS[v] || label->readonly) {
966: IS is;
967: PetscInt i;
969: PetscUseTypeMethod(label, getstratumis, v, &is);
970: PetscCall(ISLocate(is, point, &i));
971: PetscCall(ISDestroy(&is));
972: *contains = (PetscBool)(i >= 0);
973: } else {
974: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
975: }
976: } else { // value is not present
977: *contains = PETSC_FALSE;
978: }
979: }
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
985: When a label is created, it is initialized to -1.
987: Not Collective
989: Input Parameter:
990: . label - a `DMLabel` object
992: Output Parameter:
993: . defaultValue - the default value
995: Level: beginner
997: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
998: @*/
999: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
1000: {
1001: PetscFunctionBegin;
1003: *defaultValue = label->defaultValue;
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: /*@
1008: DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
1009: When a label is created, it is initialized to -1.
1011: Not Collective
1013: Input Parameter:
1014: . label - a `DMLabel` object
1016: Output Parameter:
1017: . defaultValue - the default value
1019: Level: beginner
1021: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1022: @*/
1023: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1024: {
1025: PetscFunctionBegin;
1027: label->defaultValue = defaultValue;
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: /*@
1032: 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
1033: `DMLabelSetDefaultValue()`)
1035: Not Collective
1037: Input Parameters:
1038: + label - the `DMLabel`
1039: - point - the point
1041: Output Parameter:
1042: . value - The point value, or the default value (-1 by default)
1044: Level: intermediate
1046: Note:
1047: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1048: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1050: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1051: @*/
1052: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1053: {
1054: PetscInt v;
1056: PetscFunctionBeginHot;
1058: PetscAssertPointer(value, 3);
1059: *value = label->defaultValue;
1060: for (v = 0; v < label->numStrata; ++v) {
1061: if (label->validIS[v] || label->readonly) {
1062: IS is;
1063: PetscInt i;
1065: PetscUseTypeMethod(label, getstratumis, v, &is);
1066: PetscCall(ISLocate(label->points[v], point, &i));
1067: PetscCall(ISDestroy(&is));
1068: if (i >= 0) {
1069: *value = label->stratumValues[v];
1070: break;
1071: }
1072: } else {
1073: PetscBool has;
1075: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1076: if (has) {
1077: *value = label->stratumValues[v];
1078: break;
1079: }
1080: }
1081: }
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: /*@
1086: 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
1087: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1089: Not Collective
1091: Input Parameters:
1092: + label - the `DMLabel`
1093: . point - the point
1094: - value - The point value
1096: Level: intermediate
1098: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1099: @*/
1100: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1101: {
1102: PetscInt v;
1104: PetscFunctionBegin;
1106: /* Find label value, add new entry if needed */
1107: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1108: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1109: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1110: /* Set key */
1111: PetscCall(DMLabelMakeInvalid_Private(label, v));
1112: PetscCall(PetscHSetIAdd(label->ht[v], point));
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: /*@
1117: DMLabelClearValue - Clear the value a label assigns to a point
1119: Not Collective
1121: Input Parameters:
1122: + label - the `DMLabel`
1123: . point - the point
1124: - value - The point value
1126: Level: intermediate
1128: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1129: @*/
1130: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1131: {
1132: PetscInt v;
1134: PetscFunctionBegin;
1136: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1137: /* Find label value */
1138: PetscCall(DMLabelLookupStratum(label, value, &v));
1139: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1141: if (label->bt && point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1143: /* Delete key */
1144: PetscCall(DMLabelMakeInvalid_Private(label, v));
1145: PetscCall(PetscHSetIDel(label->ht[v], point));
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: /*@
1150: DMLabelInsertIS - Set all points in the `IS` to a value
1152: Not Collective
1154: Input Parameters:
1155: + label - the `DMLabel`
1156: . is - the point `IS`
1157: - value - The point value
1159: Level: intermediate
1161: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1162: @*/
1163: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1164: {
1165: PetscInt v, n, p;
1166: const PetscInt *points;
1168: PetscFunctionBegin;
1171: /* Find label value, add new entry if needed */
1172: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1173: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1174: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1175: /* Set keys */
1176: PetscCall(DMLabelMakeInvalid_Private(label, v));
1177: PetscCall(ISGetLocalSize(is, &n));
1178: PetscCall(ISGetIndices(is, &points));
1179: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1180: PetscCall(ISRestoreIndices(is, &points));
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: /*@
1185: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1187: Not Collective
1189: Input Parameter:
1190: . label - the `DMLabel`
1192: Output Parameter:
1193: . numValues - the number of values
1195: Level: intermediate
1197: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1198: @*/
1199: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1200: {
1201: PetscFunctionBegin;
1203: PetscAssertPointer(numValues, 2);
1204: *numValues = label->numStrata;
1205: PetscFunctionReturn(PETSC_SUCCESS);
1206: }
1208: /*@
1209: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1211: Not Collective
1213: Input Parameter:
1214: . label - the `DMLabel`
1216: Output Parameter:
1217: . values - the value `IS`
1219: Level: intermediate
1221: Notes:
1222: The `values` should be destroyed when no longer needed.
1224: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1226: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1228: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1229: @*/
1230: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1231: {
1232: PetscFunctionBegin;
1234: PetscAssertPointer(values, 2);
1235: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1236: PetscFunctionReturn(PETSC_SUCCESS);
1237: }
1239: /*@
1240: DMLabelGetValueBounds - Return the smallest and largest value in the label
1242: Not Collective
1244: Input Parameter:
1245: . label - the `DMLabel`
1247: Output Parameters:
1248: + minValue - The smallest value
1249: - maxValue - The largest value
1251: Level: intermediate
1253: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1254: @*/
1255: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
1256: {
1257: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;
1259: PetscFunctionBegin;
1261: for (PetscInt v = 0; v < label->numStrata; ++v) {
1262: min = PetscMin(min, label->stratumValues[v]);
1263: max = PetscMax(max, label->stratumValues[v]);
1264: }
1265: if (minValue) {
1266: PetscAssertPointer(minValue, 2);
1267: *minValue = min;
1268: }
1269: if (maxValue) {
1270: PetscAssertPointer(maxValue, 3);
1271: *maxValue = max;
1272: }
1273: PetscFunctionReturn(PETSC_SUCCESS);
1274: }
1276: /*@
1277: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1279: Not Collective
1281: Input Parameter:
1282: . label - the `DMLabel`
1284: Output Parameter:
1285: . values - the value `IS`
1287: Level: intermediate
1289: Notes:
1290: The `values` should be destroyed when no longer needed.
1292: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1294: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1295: @*/
1296: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1297: {
1298: PetscInt i, j;
1299: PetscInt *valuesArr;
1301: PetscFunctionBegin;
1303: PetscAssertPointer(values, 2);
1304: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1305: for (i = 0, j = 0; i < label->numStrata; i++) {
1306: PetscInt n;
1308: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1309: if (n) valuesArr[j++] = label->stratumValues[i];
1310: }
1311: if (j == label->numStrata) {
1312: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1313: } else {
1314: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1315: }
1316: PetscCall(PetscFree(valuesArr));
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1320: /*@
1321: DMLabelGetValueISGlobal - Get an `IS` of all values that the `DMlabel` takes across all ranks
1323: Collective
1325: Input Parameters:
1326: + comm - MPI communicator to collect values
1327: . label - the `DMLabel`, may be `NULL` for ranks in `comm` which do not have the corresponding `DMLabel`
1328: - get_nonempty - whether to get nonempty stratum values (akin to `DMLabelGetNonEmptyStratumValuesIS()`)
1330: Output Parameter:
1331: . values - the value `IS`
1333: Level: intermediate
1335: Notes:
1336: The `values` should be destroyed when no longer needed. If no rank contributes label values, this routine returns an empty `IS`.
1338: This is similar to `DMLabelGetValueIS()` and `DMLabelGetNonEmptyStratumValuesIS()`, but gets the (nonempty) values across all ranks in `comm`.
1340: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1341: @*/
1342: PetscErrorCode DMLabelGetValueISGlobal(MPI_Comm comm, DMLabel label, PetscBool get_nonempty, IS *values)
1343: {
1344: PetscInt num_values_local = 0, num_values_global, value_range, minmax_values[2], minmax_values_loc[2] = {PETSC_INT_MAX, PETSC_INT_MIN};
1345: PetscInt *values_global;
1346: const PetscInt *values_local = NULL;
1347: IS is_values = NULL;
1348: PetscBT global_values_bt;
1350: PetscFunctionBegin;
1352: if (PetscDefined(USE_DEBUG)) {
1354: IS dummy;
1355: PetscCall(ISCreate(comm, &dummy));
1357: PetscCall(ISDestroy(&dummy));
1358: }
1359: PetscAssertPointer(values, 4);
1361: if (label) {
1362: if (get_nonempty) PetscCall(DMLabelGetNonEmptyStratumValuesIS(label, &is_values));
1363: else PetscCall(DMLabelGetValueIS(label, &is_values));
1364: PetscCall(ISGetIndices(is_values, &values_local));
1365: PetscCall(ISGetLocalSize(is_values, &num_values_local));
1366: }
1367: for (PetscInt i = 0; i < num_values_local; i++) {
1368: minmax_values_loc[0] = PetscMin(minmax_values_loc[0], values_local[i]);
1369: minmax_values_loc[1] = PetscMax(minmax_values_loc[1], values_local[i]);
1370: }
1372: PetscCall(PetscGlobalMinMaxInt(comm, minmax_values_loc, minmax_values));
1373: // The global range is empty when no rank has any values.
1374: if (minmax_values[0] > minmax_values[1]) value_range = 0;
1375: else {
1376: PetscCheck(minmax_values[0] >= 0 || minmax_values[1] <= minmax_values[0] + (PETSC_INT_MAX - 1), comm, PETSC_ERR_SUP, "Global label value range [%" PetscInt_FMT ", %" PetscInt_FMT "] is too large", minmax_values[0], minmax_values[1]);
1377: PetscCheck(minmax_values[0] < 0 || minmax_values[1] - minmax_values[0] < PETSC_INT_MAX, comm, PETSC_ERR_SUP, "Global label value range [%" PetscInt_FMT ", %" PetscInt_FMT "] is too large", minmax_values[0], minmax_values[1]);
1378: value_range = minmax_values[1] - minmax_values[0] + 1;
1379: }
1381: // Create a "ballot" where each rank marks which values they have into the PetscBT.
1382: // An Allreduce using bitwise-OR over the ranks then communicates which values are owned by a rank in comm
1383: PetscCall(PetscBTCreate(value_range, &global_values_bt));
1384: for (PetscInt i = 0; i < num_values_local; i++) PetscCall(PetscBTSet(global_values_bt, values_local[i] - minmax_values[0]));
1385: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, global_values_bt, PetscBTLength(value_range), MPI_CHAR, MPI_BOR, comm));
1386: {
1387: PetscCount num_values_global_count;
1388: num_values_global_count = PetscBTCountSet(global_values_bt, value_range);
1389: PetscCall(PetscIntCast(num_values_global_count, &num_values_global));
1390: }
1392: PetscCall(PetscMalloc1(num_values_global, &values_global));
1393: for (PetscInt i = 0, a = 0; i < value_range; i++) {
1394: if (PetscBTLookup(global_values_bt, i)) {
1395: values_global[a] = i + minmax_values[0];
1396: a++;
1397: }
1398: }
1399: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_values_global, values_global, PETSC_OWN_POINTER, values));
1401: PetscCall(PetscBTDestroy(&global_values_bt));
1402: if (is_values) {
1403: PetscCall(ISRestoreIndices(is_values, &values_local));
1404: PetscCall(ISDestroy(&is_values));
1405: }
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: /*@
1410: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1412: Not Collective
1414: Input Parameters:
1415: + label - the `DMLabel`
1416: - value - the value
1418: Output Parameter:
1419: . index - the index of value in the list of values
1421: Level: intermediate
1423: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1424: @*/
1425: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1426: {
1427: PetscInt v;
1429: PetscFunctionBegin;
1431: PetscAssertPointer(index, 3);
1432: /* Do not assume they are sorted */
1433: for (v = 0; v < label->numStrata; ++v)
1434: if (label->stratumValues[v] == value) break;
1435: if (v >= label->numStrata) *index = -1;
1436: else *index = v;
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: /*@
1441: DMLabelHasStratum - Determine whether points exist with the given value
1443: Not Collective
1445: Input Parameters:
1446: + label - the `DMLabel`
1447: - value - the stratum value
1449: Output Parameter:
1450: . exists - Flag saying whether points exist
1452: Level: intermediate
1454: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1455: @*/
1456: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1457: {
1458: PetscInt v;
1460: PetscFunctionBegin;
1462: PetscAssertPointer(exists, 3);
1463: PetscCall(DMLabelLookupStratum(label, value, &v));
1464: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: /*@
1469: DMLabelGetStratumSize - Get the size of a stratum
1471: Not Collective
1473: Input Parameters:
1474: + label - the `DMLabel`
1475: - value - the stratum value
1477: Output Parameter:
1478: . size - The number of points in the stratum
1480: Level: intermediate
1482: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1483: @*/
1484: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1485: {
1486: PetscInt v;
1488: PetscFunctionBegin;
1490: PetscAssertPointer(size, 3);
1491: PetscCall(DMLabelLookupStratum(label, value, &v));
1492: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1493: PetscFunctionReturn(PETSC_SUCCESS);
1494: }
1496: /*@
1497: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1499: Not Collective
1501: Input Parameters:
1502: + label - the `DMLabel`
1503: - value - the stratum value
1505: Output Parameters:
1506: + start - the smallest point in the stratum
1507: - end - the largest point in the stratum
1509: Level: intermediate
1511: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1512: @*/
1513: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1514: {
1515: IS is;
1516: PetscInt v, min, max;
1518: PetscFunctionBegin;
1520: if (start) {
1521: PetscAssertPointer(start, 3);
1522: *start = -1;
1523: }
1524: if (end) {
1525: PetscAssertPointer(end, 4);
1526: *end = -1;
1527: }
1528: PetscCall(DMLabelLookupStratum(label, value, &v));
1529: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1530: PetscCall(DMLabelMakeValid_Private(label, v));
1531: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1532: PetscUseTypeMethod(label, getstratumis, v, &is);
1533: PetscCall(ISGetMinMax(is, &min, &max));
1534: PetscCall(ISDestroy(&is));
1535: if (start) *start = min;
1536: if (end) *end = max + 1;
1537: PetscFunctionReturn(PETSC_SUCCESS);
1538: }
1540: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1541: {
1542: PetscFunctionBegin;
1543: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1544: *pointIS = label->points[v];
1545: PetscFunctionReturn(PETSC_SUCCESS);
1546: }
1548: /*@
1549: DMLabelGetStratumIS - Get an `IS` with the stratum points
1551: Not Collective
1553: Input Parameters:
1554: + label - the `DMLabel`
1555: - value - the stratum value
1557: Output Parameter:
1558: . points - The stratum points
1560: Level: intermediate
1562: Notes:
1563: The output `IS` should be destroyed when no longer needed.
1564: Returns `NULL` if the stratum is empty.
1566: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1567: @*/
1568: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1569: {
1570: PetscInt v;
1572: PetscFunctionBegin;
1574: PetscAssertPointer(points, 3);
1575: *points = NULL;
1576: PetscCall(DMLabelLookupStratum(label, value, &v));
1577: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1578: PetscCall(DMLabelMakeValid_Private(label, v));
1579: PetscUseTypeMethod(label, getstratumis, v, points);
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: DMLabelSetStratumIS - Set the stratum points using an `IS`
1586: Not Collective
1588: Input Parameters:
1589: + label - the `DMLabel`
1590: . value - the stratum value
1591: - is - The stratum points
1593: Level: intermediate
1595: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1596: @*/
1597: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1598: {
1599: PetscInt v;
1601: PetscFunctionBegin;
1604: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1605: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1606: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1607: PetscCall(DMLabelClearStratum(label, value));
1608: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1609: PetscCall(PetscObjectReference((PetscObject)is));
1610: PetscCall(ISDestroy(&label->points[v]));
1611: label->points[v] = is;
1612: label->validIS[v] = PETSC_TRUE;
1613: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1614: if (label->bt) {
1615: const PetscInt *points;
1617: PetscCall(ISGetIndices(is, &points));
1618: for (PetscInt p = 0; p < label->stratumSizes[v]; ++p) {
1619: const PetscInt point = points[p];
1621: 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);
1622: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1623: }
1624: }
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: /*@
1629: DMLabelClearStratum - Remove a stratum
1631: Not Collective
1633: Input Parameters:
1634: + label - the `DMLabel`
1635: - value - the stratum value
1637: Level: intermediate
1639: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1640: @*/
1641: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1642: {
1643: PetscInt v;
1645: PetscFunctionBegin;
1647: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1648: PetscCall(DMLabelLookupStratum(label, value, &v));
1649: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1650: if (label->validIS[v]) {
1651: if (label->bt) {
1652: const PetscInt *points;
1654: PetscCall(ISGetIndices(label->points[v], &points));
1655: for (PetscInt i = 0; i < label->stratumSizes[v]; ++i) {
1656: const PetscInt point = points[i];
1658: if (point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1659: }
1660: PetscCall(ISRestoreIndices(label->points[v], &points));
1661: }
1662: label->stratumSizes[v] = 0;
1663: PetscCall(ISDestroy(&label->points[v]));
1664: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1665: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1666: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1667: } else {
1668: PetscCall(PetscHSetIClear(label->ht[v]));
1669: }
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: /*@
1674: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1676: Not Collective
1678: Input Parameters:
1679: + label - The `DMLabel`
1680: . value - The label value for all points
1681: . pStart - The first point
1682: - pEnd - A point beyond all marked points
1684: Level: intermediate
1686: Note:
1687: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1689: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1690: @*/
1691: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1692: {
1693: IS pIS;
1695: PetscFunctionBegin;
1696: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1697: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1698: PetscCall(ISDestroy(&pIS));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: /*@
1703: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1705: Not Collective
1707: Input Parameters:
1708: + label - The `DMLabel`
1709: . value - The label value
1710: - p - A point with this value
1712: Output Parameter:
1713: . 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
1715: Level: intermediate
1717: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1718: @*/
1719: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1720: {
1721: IS pointIS;
1722: PetscInt v;
1724: PetscFunctionBegin;
1726: PetscAssertPointer(index, 4);
1727: *index = -1;
1728: PetscCall(DMLabelLookupStratum(label, value, &v));
1729: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1730: PetscCall(DMLabelMakeValid_Private(label, v));
1731: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1732: PetscCall(ISLocate(pointIS, p, index));
1733: PetscCall(ISDestroy(&pointIS));
1734: PetscFunctionReturn(PETSC_SUCCESS);
1735: }
1737: /*@
1738: DMLabelFilter - Remove all points outside of [`start`, `end`)
1740: Not Collective
1742: Input Parameters:
1743: + label - the `DMLabel`
1744: . start - the first point kept
1745: - end - one more than the last point kept
1747: Level: intermediate
1749: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1750: @*/
1751: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1752: {
1753: PetscInt v;
1755: PetscFunctionBegin;
1757: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1758: PetscCall(DMLabelDestroyIndex(label));
1759: PetscCall(DMLabelMakeAllValid_Private(label));
1760: for (v = 0; v < label->numStrata; ++v) {
1761: PetscCall(ISGeneralFilter(label->points[v], start, end));
1762: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1763: }
1764: PetscCall(DMLabelCreateIndex(label, start, end));
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: /*@
1769: DMLabelPermute - Create a new label with permuted points
1771: Not Collective
1773: Input Parameters:
1774: + label - the `DMLabel`
1775: - permutation - the point permutation
1777: Output Parameter:
1778: . labelNew - the new label containing the permuted points
1780: Level: intermediate
1782: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1783: @*/
1784: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1785: {
1786: const PetscInt *perm;
1787: PetscInt numValues, numPoints, v, q;
1789: PetscFunctionBegin;
1792: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1793: PetscCall(DMLabelMakeAllValid_Private(label));
1794: PetscCall(DMLabelDuplicate(label, labelNew));
1795: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1796: PetscCall(ISGetLocalSize(permutation, &numPoints));
1797: PetscCall(ISGetIndices(permutation, &perm));
1798: for (v = 0; v < numValues; ++v) {
1799: const PetscInt size = (*labelNew)->stratumSizes[v];
1800: const PetscInt *points;
1801: PetscInt *pointsNew;
1803: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1804: PetscCall(PetscCalloc1(size, &pointsNew));
1805: for (q = 0; q < size; ++q) {
1806: const PetscInt point = points[q];
1808: 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);
1809: pointsNew[q] = perm[point];
1810: }
1811: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1812: PetscCall(PetscSortInt(size, pointsNew));
1813: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1814: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1815: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1816: PetscCall(PetscFree(pointsNew));
1817: } else {
1818: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1819: }
1820: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1821: }
1822: PetscCall(ISRestoreIndices(permutation, &perm));
1823: if (label->bt) {
1824: PetscCall(PetscBTDestroy(&label->bt));
1825: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1826: }
1827: PetscFunctionReturn(PETSC_SUCCESS);
1828: }
1830: /*@
1831: DMLabelPermuteValues - Permute the values in a label
1833: Not collective
1835: Input Parameters:
1836: + label - the `DMLabel`
1837: - permutation - the value permutation, permutation[old value] = new value
1839: Output Parameter:
1840: . label - the `DMLabel` now with permuted values
1842: Note:
1843: The modification is done in-place
1845: Level: intermediate
1847: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1848: @*/
1849: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1850: {
1851: PetscInt Nv, Np;
1853: PetscFunctionBegin;
1856: PetscCall(DMLabelGetNumValues(label, &Nv));
1857: PetscCall(ISGetLocalSize(permutation, &Np));
1858: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1859: if (PetscDefined(USE_DEBUG)) {
1860: PetscBool flg;
1861: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1862: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1863: }
1864: PetscCall(DMLabelRewriteValues(label, permutation));
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: /*@
1869: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1871: Not collective
1873: Input Parameters:
1874: + label - the `DMLabel`
1875: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1877: Output Parameter:
1878: . label - the `DMLabel` now with permuted values
1880: Note:
1881: The modification is done in-place
1883: Level: intermediate
1885: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1886: @*/
1887: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1888: {
1889: const PetscInt *perm;
1890: PetscInt Nv, Np;
1892: PetscFunctionBegin;
1895: PetscCall(DMLabelMakeAllValid_Private(label));
1896: PetscCall(DMLabelGetNumValues(label, &Nv));
1897: PetscCall(ISGetLocalSize(permutation, &Np));
1898: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1899: PetscCall(ISGetIndices(permutation, &perm));
1900: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1901: PetscCall(ISRestoreIndices(permutation, &perm));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1906: {
1907: MPI_Comm comm;
1908: PetscInt s, l, nroots, nleaves, offset, size;
1909: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1910: PetscSection rootSection;
1911: PetscSF labelSF;
1913: PetscFunctionBegin;
1914: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1915: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1916: /* Build a section of stratum values per point, generate the according SF
1917: and distribute point-wise stratum values to leaves. */
1918: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1919: PetscCall(PetscSectionCreate(comm, &rootSection));
1920: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1921: if (label) {
1922: for (s = 0; s < label->numStrata; ++s) {
1923: const PetscInt *points;
1925: PetscCall(ISGetIndices(label->points[s], &points));
1926: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1927: PetscCall(ISRestoreIndices(label->points[s], &points));
1928: }
1929: }
1930: PetscCall(PetscSectionSetUp(rootSection));
1931: /* Create a point-wise array of stratum values */
1932: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1933: PetscCall(PetscMalloc1(size, &rootStrata));
1934: PetscCall(PetscCalloc1(nroots, &rootIdx));
1935: if (label) {
1936: for (s = 0; s < label->numStrata; ++s) {
1937: const PetscInt *points;
1939: PetscCall(ISGetIndices(label->points[s], &points));
1940: for (l = 0; l < label->stratumSizes[s]; l++) {
1941: const PetscInt p = points[l];
1942: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1943: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1944: }
1945: PetscCall(ISRestoreIndices(label->points[s], &points));
1946: }
1947: }
1948: /* Build SF that maps label points to remote processes */
1949: PetscCall(PetscSectionCreate(comm, leafSection));
1950: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1951: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1952: PetscCall(PetscFree(remoteOffsets));
1953: /* Send the strata for each point over the derived SF */
1954: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1955: PetscCall(PetscMalloc1(size, leafStrata));
1956: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1957: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1958: /* Clean up */
1959: PetscCall(PetscFree(rootStrata));
1960: PetscCall(PetscFree(rootIdx));
1961: PetscCall(PetscSectionDestroy(&rootSection));
1962: PetscCall(PetscSFDestroy(&labelSF));
1963: PetscFunctionReturn(PETSC_SUCCESS);
1964: }
1966: /*@
1967: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1969: Collective
1971: Input Parameters:
1972: + label - the `DMLabel`
1973: - sf - the map from old to new distribution
1975: Output Parameter:
1976: . labelNew - the new redistributed label
1978: Level: intermediate
1980: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1981: @*/
1982: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1983: {
1984: MPI_Comm comm;
1985: PetscSection leafSection;
1986: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1987: PetscInt *leafStrata, *strataIdx;
1988: PetscInt **points;
1989: const char *lname = NULL;
1990: char *name;
1991: PetscMPIInt nameSize;
1992: PetscHSetI stratumHash;
1993: size_t len = 0;
1994: PetscMPIInt rank;
1996: PetscFunctionBegin;
1998: if (label) {
2000: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2001: PetscCall(DMLabelMakeAllValid_Private(label));
2002: }
2003: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2004: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2005: /* Bcast name */
2006: if (rank == 0) {
2007: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2008: PetscCall(PetscStrlen(lname, &len));
2009: }
2010: PetscCall(PetscMPIIntCast(len, &nameSize));
2011: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2012: PetscCall(PetscMalloc1(nameSize + 1, &name));
2013: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2014: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2015: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2016: PetscCall(PetscFree(name));
2017: /* Bcast defaultValue */
2018: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
2019: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
2020: /* Distribute stratum values over the SF and get the point mapping on the receiver */
2021: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
2022: /* Determine received stratum values and initialise new label*/
2023: PetscCall(PetscHSetICreate(&stratumHash));
2024: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
2025: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
2026: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
2027: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
2028: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
2029: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
2030: /* Turn leafStrata into indices rather than stratum values */
2031: offset = 0;
2032: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
2033: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
2034: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
2035: for (p = 0; p < size; ++p) {
2036: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2037: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
2038: leafStrata[p] = s;
2039: break;
2040: }
2041: }
2042: }
2043: /* Rebuild the point strata on the receiver */
2044: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
2045: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2046: for (p = pStart; p < pEnd; p++) {
2047: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2048: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2049: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
2050: }
2051: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
2052: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
2053: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
2054: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2055: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
2056: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
2057: }
2058: /* Insert points into new strata */
2059: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
2060: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2061: for (p = pStart; p < pEnd; p++) {
2062: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2063: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2064: for (s = 0; s < dof; s++) {
2065: stratum = leafStrata[offset + s];
2066: points[stratum][strataIdx[stratum]++] = p;
2067: }
2068: }
2069: for (s = 0; s < (*labelNew)->numStrata; s++) {
2070: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
2071: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
2072: }
2073: PetscCall(PetscFree(points));
2074: PetscCall(PetscHSetIDestroy(&stratumHash));
2075: PetscCall(PetscFree(leafStrata));
2076: PetscCall(PetscFree(strataIdx));
2077: PetscCall(PetscSectionDestroy(&leafSection));
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: /*@
2082: DMLabelGather - Gather all label values from leafs into roots
2084: Collective
2086: Input Parameters:
2087: + label - the `DMLabel`
2088: - sf - the `PetscSF` communication map
2090: Output Parameter:
2091: . labelNew - the new `DMLabel` with localised leaf values
2093: Level: developer
2095: Note:
2096: This is the inverse operation to `DMLabelDistribute()`.
2098: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2099: @*/
2100: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
2101: {
2102: MPI_Comm comm;
2103: PetscSection rootSection;
2104: PetscSF sfLabel;
2105: PetscSFNode *rootPoints, *leafPoints;
2106: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2107: const PetscInt *rootDegree, *ilocal;
2108: PetscInt *rootStrata;
2109: const char *lname;
2110: char *name;
2111: PetscMPIInt nameSize;
2112: size_t len = 0;
2113: PetscMPIInt rank, size;
2115: PetscFunctionBegin;
2118: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2119: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2120: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2121: PetscCallMPI(MPI_Comm_size(comm, &size));
2122: /* Bcast name */
2123: if (rank == 0) {
2124: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2125: PetscCall(PetscStrlen(lname, &len));
2126: }
2127: PetscCall(PetscMPIIntCast(len, &nameSize));
2128: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2129: PetscCall(PetscMalloc1(nameSize + 1, &name));
2130: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2131: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2132: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2133: PetscCall(PetscFree(name));
2134: /* Gather rank/index pairs of leaves into local roots to build
2135: an inverse, multi-rooted SF. Note that this ignores local leaf
2136: indexing due to the use of the multiSF in PetscSFGather. */
2137: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2138: PetscCall(PetscMalloc1(nroots, &leafPoints));
2139: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2140: for (p = 0; p < nleaves; p++) {
2141: PetscInt ilp = ilocal ? ilocal[p] : p;
2143: leafPoints[ilp].index = ilp;
2144: leafPoints[ilp].rank = rank;
2145: }
2146: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2147: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2148: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2149: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2150: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2151: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2152: PetscCall(PetscSFCreate(comm, &sfLabel));
2153: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2154: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2155: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2156: /* Rebuild the point strata on the receiver */
2157: for (p = 0, idx = 0; p < nroots; p++) {
2158: for (d = 0; d < rootDegree[p]; d++) {
2159: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2160: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2161: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2162: }
2163: idx += rootDegree[p];
2164: }
2165: PetscCall(PetscFree(leafPoints));
2166: PetscCall(PetscFree(rootStrata));
2167: PetscCall(PetscSectionDestroy(&rootSection));
2168: PetscCall(PetscSFDestroy(&sfLabel));
2169: PetscFunctionReturn(PETSC_SUCCESS);
2170: }
2172: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2173: {
2174: const PetscInt *degree;
2175: const PetscInt *points;
2176: PetscInt Nr, r, Nl, l, val, defVal;
2178: PetscFunctionBegin;
2179: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2180: /* Add in leaves */
2181: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2182: for (l = 0; l < Nl; ++l) {
2183: PetscCall(DMLabelGetValue(label, points[l], &val));
2184: if (val != defVal) valArray[points[l]] = val;
2185: }
2186: /* Add in shared roots */
2187: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2188: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2189: for (r = 0; r < Nr; ++r) {
2190: if (degree[r]) {
2191: PetscCall(DMLabelGetValue(label, r, &val));
2192: if (val != defVal) valArray[r] = val;
2193: }
2194: }
2195: PetscFunctionReturn(PETSC_SUCCESS);
2196: }
2198: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), PetscCtx ctx)
2199: {
2200: const PetscInt *degree;
2201: const PetscInt *points;
2202: PetscInt Nr, r, Nl, l, val, defVal;
2204: PetscFunctionBegin;
2205: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2206: /* Read out leaves */
2207: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2208: for (l = 0; l < Nl; ++l) {
2209: const PetscInt p = points[l];
2210: const PetscInt cval = valArray[p];
2212: if (cval != defVal) {
2213: PetscCall(DMLabelGetValue(label, p, &val));
2214: if (val == defVal) {
2215: PetscCall(DMLabelSetValue(label, p, cval));
2216: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2217: }
2218: }
2219: }
2220: /* Read out shared roots */
2221: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2222: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2223: for (r = 0; r < Nr; ++r) {
2224: if (degree[r]) {
2225: const PetscInt cval = valArray[r];
2227: if (cval != defVal) {
2228: PetscCall(DMLabelGetValue(label, r, &val));
2229: if (val == defVal) {
2230: PetscCall(DMLabelSetValue(label, r, cval));
2231: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2232: }
2233: }
2234: }
2235: }
2236: PetscFunctionReturn(PETSC_SUCCESS);
2237: }
2239: /*@
2240: DMLabelPropagateBegin - Setup a cycle of label propagation
2242: Collective
2244: Input Parameters:
2245: + label - The `DMLabel` to propagate across processes
2246: - sf - The `PetscSF` describing parallel layout of the label points
2248: Level: intermediate
2250: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2251: @*/
2252: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2253: {
2254: PetscInt Nr, r, defVal;
2255: PetscMPIInt size;
2257: PetscFunctionBegin;
2258: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2259: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2260: if (size > 1) {
2261: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2262: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2263: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2264: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2265: }
2266: PetscFunctionReturn(PETSC_SUCCESS);
2267: }
2269: /*@
2270: DMLabelPropagateEnd - Tear down a cycle of label propagation
2272: Collective
2274: Input Parameters:
2275: + label - The `DMLabel` to propagate across processes
2276: - pointSF - The `PetscSF` describing parallel layout of the label points
2278: Level: intermediate
2280: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2281: @*/
2282: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2283: {
2284: PetscFunctionBegin;
2285: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2286: PetscCall(PetscFree(label->propArray));
2287: label->propArray = NULL;
2288: PetscFunctionReturn(PETSC_SUCCESS);
2289: }
2291: /*@
2292: DMLabelPropagatePush - Execute a cycle of label propagation
2294: Collective
2296: Input Parameters:
2297: + label - The `DMLabel` to propagate across processes
2298: . pointSF - The `PetscSF` describing parallel layout of the label points
2299: . merge - The operator which merges label values
2300: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2301: - ctx - An optional application context for the callback, or `NULL`
2303: Calling sequence of `markPoint`:
2304: + label - The `DMLabel`
2305: . p - The point being marked
2306: . val - The label value for `p`
2307: - ctx - An optional application context
2309: Level: intermediate
2311: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2312: @*/
2313: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, MPI_Op merge, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, PetscCtx ctx), PetscCtx ctx)
2314: {
2315: PetscInt *valArray = label->propArray, Nr;
2316: PetscMPIInt size;
2318: PetscFunctionBegin;
2319: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2320: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2321: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2322: if (size > 1 && Nr >= 0) {
2323: /* Communicate marked edges
2324: The current implementation allocates an array the size of the number of root. We put the label values into the
2325: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2327: TODO: We could use in-place communication with a different SF
2328: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2329: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2331: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2332: 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
2333: edge to the queue.
2334: */
2335: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2336: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, merge));
2337: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, merge));
2338: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2339: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2340: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2341: }
2342: PetscFunctionReturn(PETSC_SUCCESS);
2343: }
2345: /*@
2346: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2348: Not Collective
2350: Input Parameter:
2351: . label - the `DMLabel`
2353: Output Parameters:
2354: + section - the section giving offsets for each stratum
2355: - is - An `IS` containing all the label points
2357: Level: developer
2359: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2360: @*/
2361: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2362: {
2363: IS vIS;
2364: const PetscInt *values;
2365: PetscInt *points;
2366: PetscInt nV, vS = 0, vE = 0, v, N;
2368: PetscFunctionBegin;
2370: PetscCall(DMLabelGetNumValues(label, &nV));
2371: PetscCall(DMLabelGetValueIS(label, &vIS));
2372: PetscCall(ISGetIndices(vIS, &values));
2373: if (nV) {
2374: vS = values[0];
2375: vE = values[0] + 1;
2376: }
2377: for (v = 1; v < nV; ++v) {
2378: vS = PetscMin(vS, values[v]);
2379: vE = PetscMax(vE, values[v] + 1);
2380: }
2381: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2382: PetscCall(PetscSectionSetChart(*section, vS, vE));
2383: for (v = 0; v < nV; ++v) {
2384: PetscInt n;
2386: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2387: PetscCall(PetscSectionSetDof(*section, values[v], n));
2388: }
2389: PetscCall(PetscSectionSetUp(*section));
2390: PetscCall(PetscSectionGetStorageSize(*section, &N));
2391: PetscCall(PetscMalloc1(N, &points));
2392: for (v = 0; v < nV; ++v) {
2393: IS is;
2394: const PetscInt *spoints;
2395: PetscInt dof, off, p;
2397: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2398: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2399: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2400: PetscCall(ISGetIndices(is, &spoints));
2401: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2402: PetscCall(ISRestoreIndices(is, &spoints));
2403: PetscCall(ISDestroy(&is));
2404: }
2405: PetscCall(ISRestoreIndices(vIS, &values));
2406: PetscCall(ISDestroy(&vIS));
2407: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2408: PetscFunctionReturn(PETSC_SUCCESS);
2409: }
2411: /*@
2412: DMLabelRegister - Adds a new label component implementation
2414: Not Collective
2416: Input Parameters:
2417: + name - The name of a new user-defined creation routine
2418: - create_func - The creation routine itself
2420: Notes:
2421: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2423: Example Usage:
2424: .vb
2425: DMLabelRegister("my_label", MyLabelCreate);
2426: .ve
2428: Then, your label type can be chosen with the procedural interface via
2429: .vb
2430: DMLabelCreate(MPI_Comm, DMLabel *);
2431: DMLabelSetType(DMLabel, "my_label");
2432: .ve
2433: or at runtime via the option
2434: .vb
2435: -dm_label_type my_label
2436: .ve
2438: Level: advanced
2440: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2441: @*/
2442: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2443: {
2444: PetscFunctionBegin;
2445: PetscCall(DMInitializePackage());
2446: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2447: PetscFunctionReturn(PETSC_SUCCESS);
2448: }
2450: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2451: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2453: /*@
2454: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2456: Not Collective
2458: Level: advanced
2460: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2461: @*/
2462: PetscErrorCode DMLabelRegisterAll(void)
2463: {
2464: PetscFunctionBegin;
2465: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2466: DMLabelRegisterAllCalled = PETSC_TRUE;
2468: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2469: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2470: PetscFunctionReturn(PETSC_SUCCESS);
2471: }
2473: /*@
2474: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2476: Level: developer
2478: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2479: @*/
2480: PetscErrorCode DMLabelRegisterDestroy(void)
2481: {
2482: PetscFunctionBegin;
2483: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2484: DMLabelRegisterAllCalled = PETSC_FALSE;
2485: PetscFunctionReturn(PETSC_SUCCESS);
2486: }
2488: /*@
2489: DMLabelSetType - Sets the particular implementation for a label.
2491: Collective
2493: Input Parameters:
2494: + label - The label
2495: - method - The name of the label type
2497: Options Database Key:
2498: . -dm_label_type type - Sets the label type; see `DMLabelType`
2500: Level: intermediate
2502: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`, `DMLabelType`
2503: @*/
2504: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2505: {
2506: PetscErrorCode (*r)(DMLabel);
2507: PetscBool match;
2509: PetscFunctionBegin;
2511: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2512: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2514: PetscCall(DMLabelRegisterAll());
2515: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2516: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2518: PetscTryTypeMethod(label, destroy);
2519: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2520: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2521: PetscCall((*r)(label));
2522: PetscFunctionReturn(PETSC_SUCCESS);
2523: }
2525: /*@
2526: DMLabelGetType - Gets the type name (as a string) from the label.
2528: Not Collective
2530: Input Parameter:
2531: . label - The `DMLabel`
2533: Output Parameter:
2534: . type - The `DMLabel` type name
2536: Level: intermediate
2538: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2539: @*/
2540: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2541: {
2542: PetscFunctionBegin;
2544: PetscAssertPointer(type, 2);
2545: PetscCall(DMLabelRegisterAll());
2546: *type = ((PetscObject)label)->type_name;
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2550: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2551: {
2552: PetscFunctionBegin;
2553: label->ops->view = DMLabelView_Concrete;
2554: label->ops->setup = NULL;
2555: label->ops->duplicate = DMLabelDuplicate_Concrete;
2556: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2557: PetscFunctionReturn(PETSC_SUCCESS);
2558: }
2560: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2561: {
2562: PetscFunctionBegin;
2564: PetscCall(DMLabelInitialize_Concrete(label));
2565: PetscFunctionReturn(PETSC_SUCCESS);
2566: }
2568: /*@
2569: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2570: the local section and an `PetscSF` describing the section point overlap.
2572: Collective
2574: Input Parameters:
2575: + s - The `PetscSection` for the local field layout
2576: . sf - The `PetscSF` describing parallel layout of the section points
2577: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2578: . label - The label specifying the points
2579: - labelValue - The label stratum specifying the points
2581: Output Parameter:
2582: . gsection - The `PetscSection` for the global field layout
2584: Level: developer
2586: Note:
2587: This gives negative sizes and offsets to points not owned by this process
2589: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2590: @*/
2591: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2592: {
2593: PetscInt *neg = NULL, *tmpOff = NULL;
2594: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2596: PetscFunctionBegin;
2600: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2601: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2602: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2603: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2604: if (nroots >= 0) {
2605: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2606: PetscCall(PetscCalloc1(nroots, &neg));
2607: if (nroots > pEnd - pStart) {
2608: PetscCall(PetscCalloc1(nroots, &tmpOff));
2609: } else {
2610: tmpOff = &(*gsection)->atlasDof[-pStart];
2611: }
2612: }
2613: /* Mark ghost points with negative dof */
2614: for (p = pStart; p < pEnd; ++p) {
2615: PetscInt value;
2617: PetscCall(DMLabelGetValue(label, p, &value));
2618: if (value != labelValue) continue;
2619: PetscCall(PetscSectionGetDof(s, p, &dof));
2620: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2621: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2622: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2623: if (neg) neg[p] = -(dof + 1);
2624: }
2625: PetscCall(PetscSectionSetUpBC(*gsection));
2626: if (nroots >= 0) {
2627: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2628: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2629: if (nroots > pEnd - pStart) {
2630: for (p = pStart; p < pEnd; ++p) {
2631: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2632: }
2633: }
2634: }
2635: /* Calculate new sizes, get process offset, and calculate point offsets */
2636: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2637: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2638: (*gsection)->atlasOff[p] = off;
2639: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2640: }
2641: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2642: globalOff -= off;
2643: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2644: (*gsection)->atlasOff[p] += globalOff;
2645: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2646: }
2647: /* Put in negative offsets for ghost points */
2648: if (nroots >= 0) {
2649: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2650: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2651: if (nroots > pEnd - pStart) {
2652: for (p = pStart; p < pEnd; ++p) {
2653: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2654: }
2655: }
2656: }
2657: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2658: PetscCall(PetscFree(neg));
2659: PetscFunctionReturn(PETSC_SUCCESS);
2660: }
2662: typedef struct _n_PetscSectionSym_Label {
2663: DMLabel label;
2664: PetscCopyMode *modes;
2665: PetscInt *sizes;
2666: const PetscInt ***perms;
2667: const PetscScalar ***rots;
2668: PetscInt (*minMaxOrients)[2];
2669: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2670: } PetscSectionSym_Label;
2672: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2673: {
2674: PetscInt i, j;
2675: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2677: PetscFunctionBegin;
2678: for (i = 0; i <= sl->numStrata; i++) {
2679: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2680: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2681: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2682: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2683: }
2684: if (sl->perms[i]) {
2685: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2687: PetscCall(PetscFree(perms));
2688: }
2689: if (sl->rots[i]) {
2690: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2692: PetscCall(PetscFree(rots));
2693: }
2694: }
2695: }
2696: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2697: PetscCall(DMLabelDestroy(&sl->label));
2698: sl->numStrata = 0;
2699: PetscFunctionReturn(PETSC_SUCCESS);
2700: }
2702: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2703: {
2704: PetscFunctionBegin;
2705: PetscCall(PetscSectionSymLabelReset(sym));
2706: PetscCall(PetscFree(sym->data));
2707: PetscFunctionReturn(PETSC_SUCCESS);
2708: }
2710: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2711: {
2712: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2713: PetscBool isAscii;
2714: DMLabel label = sl->label;
2715: const char *name;
2717: PetscFunctionBegin;
2718: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2719: if (isAscii) {
2720: PetscInt i, j, k;
2721: PetscViewerFormat format;
2723: PetscCall(PetscViewerGetFormat(viewer, &format));
2724: if (label) {
2725: PetscCall(PetscViewerGetFormat(viewer, &format));
2726: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2727: PetscCall(PetscViewerASCIIPushTab(viewer));
2728: PetscCall(DMLabelView(label, viewer));
2729: PetscCall(PetscViewerASCIIPopTab(viewer));
2730: } else {
2731: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2732: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2733: }
2734: } else {
2735: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2736: }
2737: PetscCall(PetscViewerASCIIPushTab(viewer));
2738: for (i = 0; i <= sl->numStrata; i++) {
2739: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2741: if (!(sl->perms[i] || sl->rots[i])) {
2742: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2743: } else {
2744: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2745: PetscCall(PetscViewerASCIIPushTab(viewer));
2746: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2747: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2748: PetscCall(PetscViewerASCIIPushTab(viewer));
2749: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2750: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2751: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2752: } else {
2753: PetscInt tab;
2755: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2756: PetscCall(PetscViewerASCIIPushTab(viewer));
2757: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2758: if (sl->perms[i] && sl->perms[i][j]) {
2759: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2760: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2761: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2762: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2763: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2764: }
2765: if (sl->rots[i] && sl->rots[i][j]) {
2766: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2767: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2768: #if PetscDefined(USE_COMPLEX)
2769: 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])));
2770: #else
2771: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2772: #endif
2773: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2774: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2775: }
2776: PetscCall(PetscViewerASCIIPopTab(viewer));
2777: }
2778: }
2779: PetscCall(PetscViewerASCIIPopTab(viewer));
2780: }
2781: PetscCall(PetscViewerASCIIPopTab(viewer));
2782: }
2783: }
2784: PetscCall(PetscViewerASCIIPopTab(viewer));
2785: }
2786: PetscFunctionReturn(PETSC_SUCCESS);
2787: }
2789: /*@
2790: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2792: Logically
2794: Input Parameters:
2795: + sym - the section symmetries
2796: - label - the `DMLabel` describing the types of points
2798: Level: developer:
2800: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2801: @*/
2802: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2803: {
2804: PetscSectionSym_Label *sl;
2806: PetscFunctionBegin;
2808: sl = (PetscSectionSym_Label *)sym->data;
2809: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2810: if (label) {
2811: sl->label = label;
2812: PetscCall(PetscObjectReference((PetscObject)label));
2813: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2814: 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));
2815: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2816: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2817: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2818: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2819: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2820: }
2821: PetscFunctionReturn(PETSC_SUCCESS);
2822: }
2824: /*@
2825: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2827: Logically Collective
2829: Input Parameters:
2830: + sym - the section symmetries
2831: - stratum - the stratum value in the label that we are assigning symmetries for
2833: Output Parameters:
2834: + size - the number of dofs for points in the `stratum` of the label
2835: . minOrient - the smallest orientation for a point in this `stratum`
2836: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2837: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2838: - 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
2840: Level: developer
2842: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2843: @*/
2844: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2845: {
2846: PetscSectionSym_Label *sl;
2847: const char *name;
2848: PetscInt i;
2850: PetscFunctionBegin;
2852: sl = (PetscSectionSym_Label *)sym->data;
2853: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2854: for (i = 0; i <= sl->numStrata; i++) {
2855: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2857: if (stratum == value) break;
2858: }
2859: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2860: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2861: if (size) {
2862: PetscAssertPointer(size, 3);
2863: *size = sl->sizes[i];
2864: }
2865: if (minOrient) {
2866: PetscAssertPointer(minOrient, 4);
2867: *minOrient = sl->minMaxOrients[i][0];
2868: }
2869: if (maxOrient) {
2870: PetscAssertPointer(maxOrient, 5);
2871: *maxOrient = sl->minMaxOrients[i][1];
2872: }
2873: if (perms) {
2874: PetscAssertPointer(perms, 6);
2875: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2876: }
2877: if (rots) {
2878: PetscAssertPointer(rots, 7);
2879: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2880: }
2881: PetscFunctionReturn(PETSC_SUCCESS);
2882: }
2884: /*@
2885: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2887: Logically
2889: Input Parameters:
2890: + sym - the section symmetries
2891: . stratum - the stratum value in the label that we are assigning symmetries for
2892: . size - the number of dofs for points in the `stratum` of the label
2893: . minOrient - the smallest orientation for a point in this `stratum`
2894: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2895: . mode - how `sym` should copy the `perms` and `rots` arrays
2896: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2897: - 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
2899: Level: developer
2901: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2902: @*/
2903: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2904: {
2905: PetscSectionSym_Label *sl;
2906: const char *name;
2907: PetscInt i, j, k;
2909: PetscFunctionBegin;
2911: sl = (PetscSectionSym_Label *)sym->data;
2912: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2913: for (i = 0; i <= sl->numStrata; i++) {
2914: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2916: if (stratum == value) break;
2917: }
2918: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2919: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2920: sl->sizes[i] = size;
2921: sl->modes[i] = mode;
2922: sl->minMaxOrients[i][0] = minOrient;
2923: sl->minMaxOrients[i][1] = maxOrient;
2924: if (mode == PETSC_COPY_VALUES) {
2925: if (perms) {
2926: PetscInt **ownPerms;
2928: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2929: for (j = 0; j < maxOrient - minOrient; j++) {
2930: if (perms[j]) {
2931: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2932: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2933: }
2934: }
2935: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2936: }
2937: if (rots) {
2938: PetscScalar **ownRots;
2940: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2941: for (j = 0; j < maxOrient - minOrient; j++) {
2942: if (rots[j]) {
2943: PetscCall(PetscMalloc1(size, &ownRots[j]));
2944: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2945: }
2946: }
2947: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2948: }
2949: } else {
2950: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2951: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2952: }
2953: PetscFunctionReturn(PETSC_SUCCESS);
2954: }
2956: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2957: {
2958: PetscInt i, j, numStrata;
2959: PetscSectionSym_Label *sl;
2960: DMLabel label;
2962: PetscFunctionBegin;
2963: sl = (PetscSectionSym_Label *)sym->data;
2964: numStrata = sl->numStrata;
2965: label = sl->label;
2966: for (i = 0; i < numPoints; i++) {
2967: PetscInt point = points[2 * i];
2968: PetscInt ornt = points[2 * i + 1];
2970: for (j = 0; j < numStrata; j++) {
2971: if (label->validIS[j]) {
2972: PetscInt k;
2974: PetscCall(ISLocate(label->points[j], point, &k));
2975: if (k >= 0) break;
2976: } else {
2977: PetscBool has;
2979: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2980: if (has) break;
2981: }
2982: }
2983: 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],
2984: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2985: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2986: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2987: }
2988: PetscFunctionReturn(PETSC_SUCCESS);
2989: }
2991: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2992: {
2993: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2994: IS valIS;
2995: const PetscInt *values;
2996: PetscInt Nv;
2998: PetscFunctionBegin;
2999: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
3000: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
3001: PetscCall(ISGetIndices(valIS, &values));
3002: for (PetscInt v = 0; v < Nv; ++v) {
3003: const PetscInt val = values[v];
3004: PetscInt size, minOrient, maxOrient;
3005: const PetscInt **perms;
3006: const PetscScalar **rots;
3008: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
3009: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
3010: }
3011: PetscCall(ISDestroy(&valIS));
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3016: {
3017: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
3018: DMLabel dlabel;
3020: PetscFunctionBegin;
3021: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
3022: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
3023: PetscCall(DMLabelDestroy(&dlabel));
3024: PetscCall(PetscSectionSymCopy(sym, *dsym));
3025: PetscFunctionReturn(PETSC_SUCCESS);
3026: }
3028: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
3029: {
3030: PetscSectionSym_Label *sl;
3032: PetscFunctionBegin;
3033: PetscCall(PetscNew(&sl));
3034: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
3035: sym->ops->distribute = PetscSectionSymDistribute_Label;
3036: sym->ops->copy = PetscSectionSymCopy_Label;
3037: sym->ops->view = PetscSectionSymView_Label;
3038: sym->ops->destroy = PetscSectionSymDestroy_Label;
3039: sym->data = (void *)sl;
3040: PetscFunctionReturn(PETSC_SUCCESS);
3041: }
3043: /*@
3044: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
3046: Collective
3048: Input Parameters:
3049: + comm - the MPI communicator for the new symmetry
3050: - label - the label defining the strata
3052: Output Parameter:
3053: . sym - the section symmetries
3055: Level: developer
3057: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3058: @*/
3059: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
3060: {
3061: PetscFunctionBegin;
3062: PetscCall(DMInitializePackage());
3063: PetscCall(PetscSectionSymCreate(comm, sym));
3064: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
3065: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
3066: PetscFunctionReturn(PETSC_SUCCESS);
3067: }