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