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: PetscInt p;
166: const PetscInt *points;
168: PetscFunctionBegin;
169: if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
170: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeInvalid_Private", v);
171: if (label->points[v]) {
172: PetscCall(ISGetIndices(label->points[v], &points));
173: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
174: PetscCall(ISRestoreIndices(label->points[v], &points));
175: PetscCall(ISDestroy(&label->points[v]));
176: }
177: label->validIS[v] = PETSC_FALSE;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
182: {
183: PetscInt v;
185: PetscFunctionBegin;
186: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeInvalid_Private(label, v));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
191: #define DMLABEL_LOOKUP_THRESHOLD 16
192: #endif
194: PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
195: {
196: PetscInt v;
198: PetscFunctionBegin;
199: *index = -1;
200: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD || label->readonly) {
201: for (v = 0; v < label->numStrata; ++v)
202: if (label->stratumValues[v] == value) {
203: *index = v;
204: break;
205: }
206: } else {
207: PetscCall(PetscHMapIGet(label->hmap, value, index));
208: }
209: if (PetscDefined(USE_DEBUG) && !label->readonly) { /* Check strata hash map consistency */
210: PetscInt len, loc = -1;
211: PetscCall(PetscHMapIGetSize(label->hmap, &len));
212: PetscCheck(len == label->numStrata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
213: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
214: PetscCall(PetscHMapIGet(label->hmap, value, &loc));
215: } else {
216: for (v = 0; v < label->numStrata; ++v)
217: if (label->stratumValues[v] == value) {
218: loc = v;
219: break;
220: }
221: }
222: PetscCheck(loc == *index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
223: }
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
228: {
229: PetscInt v;
230: PetscInt *tmpV;
231: PetscInt *tmpS;
232: PetscHSetI *tmpH, ht;
233: IS *tmpP, is;
234: PetscBool *tmpB;
235: PetscHMapI hmap = label->hmap;
237: PetscFunctionBegin;
238: v = label->numStrata;
239: tmpV = label->stratumValues;
240: tmpS = label->stratumSizes;
241: tmpH = label->ht;
242: tmpP = label->points;
243: tmpB = label->validIS;
244: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
245: PetscInt *oldV = tmpV;
246: PetscInt *oldS = tmpS;
247: PetscHSetI *oldH = tmpH;
248: IS *oldP = tmpP;
249: PetscBool *oldB = tmpB;
250: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV));
251: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS));
252: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpH), &tmpH));
253: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpP), &tmpP));
254: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB));
255: PetscCall(PetscArraycpy(tmpV, oldV, v));
256: PetscCall(PetscArraycpy(tmpS, oldS, v));
257: PetscCall(PetscArraycpy(tmpH, oldH, v));
258: PetscCall(PetscArraycpy(tmpP, oldP, v));
259: PetscCall(PetscArraycpy(tmpB, oldB, v));
260: PetscCall(PetscFree(oldV));
261: PetscCall(PetscFree(oldS));
262: PetscCall(PetscFree(oldH));
263: PetscCall(PetscFree(oldP));
264: PetscCall(PetscFree(oldB));
265: }
266: label->numStrata = v + 1;
267: label->stratumValues = tmpV;
268: label->stratumSizes = tmpS;
269: label->ht = tmpH;
270: label->points = tmpP;
271: label->validIS = tmpB;
272: PetscCall(PetscHSetICreate(&ht));
273: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
274: PetscCall(PetscHMapISet(hmap, value, v));
275: tmpV[v] = value;
276: tmpS[v] = 0;
277: tmpH[v] = ht;
278: tmpP[v] = is;
279: tmpB[v] = PETSC_TRUE;
280: PetscCall(PetscObjectStateIncrease((PetscObject)label));
281: *index = v;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
286: {
287: PetscFunctionBegin;
288: PetscCall(DMLabelLookupStratum(label, value, index));
289: if (*index < 0) PetscCall(DMLabelNewStratum(label, value, index));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
294: {
295: PetscFunctionBegin;
296: *size = 0;
297: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
298: if (label->readonly || label->validIS[v]) {
299: *size = label->stratumSizes[v];
300: } else {
301: PetscCall(PetscHSetIGetSize(label->ht[v], size));
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: DMLabelAddStratum - Adds a new stratum value in a `DMLabel`
309: Input Parameters:
310: + label - The `DMLabel`
311: - value - The stratum value
313: Level: beginner
315: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
316: @*/
317: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
318: {
319: PetscInt v;
321: PetscFunctionBegin;
323: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
324: PetscCall(DMLabelLookupAddStratum(label, value, &v));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: DMLabelAddStrata - Adds new stratum values in a `DMLabel`
331: Not Collective
333: Input Parameters:
334: + label - The `DMLabel`
335: . numStrata - The number of stratum values
336: - stratumValues - The stratum values
338: Level: beginner
340: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
341: @*/
342: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
343: {
344: PetscInt *values, v;
346: PetscFunctionBegin;
348: if (numStrata) PetscAssertPointer(stratumValues, 3);
349: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
350: PetscCall(PetscMalloc1(numStrata, &values));
351: PetscCall(PetscArraycpy(values, stratumValues, numStrata));
352: PetscCall(PetscSortRemoveDupsInt(&numStrata, values));
353: if (!label->numStrata) { /* Fast preallocation */
354: PetscInt *tmpV;
355: PetscInt *tmpS;
356: PetscHSetI *tmpH, ht;
357: IS *tmpP, is;
358: PetscBool *tmpB;
359: PetscHMapI hmap = label->hmap;
361: PetscCall(PetscMalloc1(numStrata, &tmpV));
362: PetscCall(PetscMalloc1(numStrata, &tmpS));
363: PetscCall(PetscCalloc1(numStrata, &tmpH));
364: PetscCall(PetscCalloc1(numStrata, &tmpP));
365: PetscCall(PetscMalloc1(numStrata, &tmpB));
366: label->numStrata = numStrata;
367: label->stratumValues = tmpV;
368: label->stratumSizes = tmpS;
369: label->ht = tmpH;
370: label->points = tmpP;
371: label->validIS = tmpB;
372: for (v = 0; v < numStrata; ++v) {
373: PetscCall(PetscHSetICreate(&ht));
374: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
375: PetscCall(PetscHMapISet(hmap, values[v], v));
376: tmpV[v] = values[v];
377: tmpS[v] = 0;
378: tmpH[v] = ht;
379: tmpP[v] = is;
380: tmpB[v] = PETSC_TRUE;
381: }
382: PetscCall(PetscObjectStateIncrease((PetscObject)label));
383: } else {
384: for (v = 0; v < numStrata; ++v) PetscCall(DMLabelAddStratum(label, values[v]));
385: }
386: PetscCall(PetscFree(values));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@
391: DMLabelAddStrataIS - Adds new stratum values in a `DMLabel`
393: Not Collective
395: Input Parameters:
396: + label - The `DMLabel`
397: - valueIS - Index set with stratum values
399: Level: beginner
401: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
402: @*/
403: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
404: {
405: PetscInt numStrata;
406: const PetscInt *stratumValues;
408: PetscFunctionBegin;
411: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
412: PetscCall(ISGetLocalSize(valueIS, &numStrata));
413: PetscCall(ISGetIndices(valueIS, &stratumValues));
414: PetscCall(DMLabelAddStrata(label, numStrata, stratumValues));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode DMLabelView_Concrete_Ascii(DMLabel label, PetscViewer viewer)
419: {
420: PetscInt v;
421: PetscMPIInt rank;
423: PetscFunctionBegin;
424: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
425: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
426: if (label) {
427: const char *name;
429: PetscCall(PetscObjectGetName((PetscObject)label, &name));
430: PetscCall(PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name));
431: if (label->bt) PetscCall(PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd));
432: for (v = 0; v < label->numStrata; ++v) {
433: const PetscInt value = label->stratumValues[v];
434: const PetscInt *points;
435: PetscInt p;
437: PetscCall(ISGetIndices(label->points[v], &points));
438: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
439: PetscCall(ISRestoreIndices(label->points[v], &points));
440: }
441: }
442: PetscCall(PetscViewerFlush(viewer));
443: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
448: {
449: PetscBool isascii;
451: PetscFunctionBegin;
452: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
453: if (isascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: DMLabelView - View the label
460: Collective
462: Input Parameters:
463: + label - The `DMLabel`
464: - viewer - The `PetscViewer`
466: Level: intermediate
468: .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
469: @*/
470: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
471: {
472: PetscFunctionBegin;
474: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer));
476: PetscCall(DMLabelMakeAllValid_Private(label));
477: PetscUseTypeMethod(label, view, viewer);
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: DMLabelViewFromOptions - View a `DMLabel` in a particular way based on a request in the options database
484: Collective
486: Input Parameters:
487: + label - the `DMLabel` object
488: . obj - optional object that provides the prefix for the options database (if `NULL` then the prefix in `obj` is used)
489: - name - option string that is used to activate viewing
491: Level: intermediate
493: Note:
494: See `PetscObjectViewFromOptions()` for a list of values that can be provided in the options database to determine how the `DMLabel` is viewed
496: .seealso: [](ch_dmbase), `DMLabel`, `DMLabelView()`, `PetscObjectViewFromOptions()`, `DMLabelCreate()`
497: @*/
498: PetscErrorCode DMLabelViewFromOptions(DMLabel label, PeOp PetscObject obj, const char name[])
499: {
500: PetscFunctionBegin;
502: PetscCall(PetscObjectViewFromOptions((PetscObject)label, obj, name));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: /*@
507: DMLabelReset - Destroys internal data structures in a `DMLabel`
509: Not Collective
511: Input Parameter:
512: . label - The `DMLabel`
514: Level: beginner
516: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
517: @*/
518: PetscErrorCode DMLabelReset(DMLabel label)
519: {
520: PetscInt v;
522: PetscFunctionBegin;
524: for (v = 0; v < label->numStrata; ++v) {
525: if (label->ht[v]) 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: /*@C
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 i, nValues;
701: PetscCall(DMLabelGetNumValues(l0, &nValues));
702: for (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;
757: PetscInt i;
759: PetscCall(ISGetIndices(label->points[v], &points));
760: for (i = 0; i < label->stratumSizes[v]; ++i) {
761: const PetscInt point = points[i];
763: pStart = PetscMin(point, pStart);
764: pEnd = PetscMax(point + 1, pEnd);
765: }
766: PetscCall(ISRestoreIndices(label->points[v], &points));
767: }
768: label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
769: label->pEnd = pEnd;
770: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: DMLabelCreateIndex - Create an index structure for membership determination
777: Not Collective
779: Input Parameters:
780: + label - The `DMLabel`
781: . pStart - The smallest point
782: - pEnd - The largest point + 1
784: Level: intermediate
786: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
787: @*/
788: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
789: {
790: PetscInt v;
792: PetscFunctionBegin;
794: PetscCall(DMLabelDestroyIndex(label));
795: PetscCall(DMLabelMakeAllValid_Private(label));
796: label->pStart = pStart;
797: label->pEnd = pEnd;
798: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
799: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
800: for (v = 0; v < label->numStrata; ++v) {
801: IS pointIS;
802: const PetscInt *points;
803: PetscInt i;
805: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
806: PetscCall(ISGetIndices(pointIS, &points));
807: for (i = 0; i < label->stratumSizes[v]; ++i) {
808: const PetscInt point = points[i];
810: 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);
811: PetscCall(PetscBTSet(label->bt, point - pStart));
812: }
813: PetscCall(ISRestoreIndices(label->points[v], &points));
814: PetscCall(ISDestroy(&pointIS));
815: }
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: /*@
820: DMLabelDestroyIndex - Destroy the index structure
822: Not Collective
824: Input Parameter:
825: . label - the `DMLabel`
827: Level: intermediate
829: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
830: @*/
831: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
832: {
833: PetscFunctionBegin;
835: label->pStart = -1;
836: label->pEnd = -1;
837: PetscCall(PetscBTDestroy(&label->bt));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: /*@
842: DMLabelGetBounds - Return the smallest and largest point in the label
844: Not Collective
846: Input Parameter:
847: . label - the `DMLabel`
849: Output Parameters:
850: + pStart - The smallest point
851: - pEnd - The largest point + 1
853: Level: intermediate
855: Note:
856: This will compute an index for the label if one does not exist.
858: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
859: @*/
860: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
861: {
862: PetscFunctionBegin;
864: if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
865: if (pStart) {
866: PetscAssertPointer(pStart, 2);
867: *pStart = label->pStart;
868: }
869: if (pEnd) {
870: PetscAssertPointer(pEnd, 3);
871: *pEnd = label->pEnd;
872: }
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@
877: DMLabelHasValue - Determine whether a label assigns the value to any point
879: Not Collective
881: Input Parameters:
882: + label - the `DMLabel`
883: - value - the value
885: Output Parameter:
886: . contains - Flag indicating whether the label maps this value to any point
888: Level: developer
890: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
891: @*/
892: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
893: {
894: PetscInt v;
896: PetscFunctionBegin;
898: PetscAssertPointer(contains, 3);
899: PetscCall(DMLabelLookupStratum(label, value, &v));
900: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
901: PetscFunctionReturn(PETSC_SUCCESS);
902: }
904: /*@
905: DMLabelHasPoint - Determine whether a label assigns a value to a point
907: Not Collective
909: Input Parameters:
910: + label - the `DMLabel`
911: - point - the point
913: Output Parameter:
914: . contains - Flag indicating whether the label maps this point to a value
916: Level: developer
918: Note:
919: The user must call `DMLabelCreateIndex()` before this function.
921: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
922: @*/
923: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
924: {
925: PetscInt pStart, pEnd;
927: PetscFunctionBeginHot;
929: PetscAssertPointer(contains, 3);
930: /* DMLabelGetBounds() calls DMLabelCreateIndex() only if needed */
931: PetscCall(DMLabelGetBounds(label, &pStart, &pEnd));
932: PetscCall(DMLabelMakeAllValid_Private(label));
933: *contains = point >= pStart && point < pEnd && (PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE);
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: /*@
938: DMLabelStratumHasPoint - Return true if the stratum contains a point
940: Not Collective
942: Input Parameters:
943: + label - the `DMLabel`
944: . value - the stratum value
945: - point - the point
947: Output Parameter:
948: . contains - true if the stratum contains the point
950: Level: intermediate
952: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
953: @*/
954: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
955: {
956: PetscFunctionBeginHot;
958: PetscAssertPointer(contains, 4);
959: if (value == label->defaultValue) {
960: PetscInt pointVal;
962: PetscCall(DMLabelGetValue(label, point, &pointVal));
963: *contains = (PetscBool)(pointVal == value);
964: } else {
965: PetscInt v;
967: PetscCall(DMLabelLookupStratum(label, value, &v));
968: if (v >= 0) {
969: if (label->validIS[v] || label->readonly) {
970: IS is;
971: PetscInt i;
973: PetscUseTypeMethod(label, getstratumis, v, &is);
974: PetscCall(ISLocate(is, point, &i));
975: PetscCall(ISDestroy(&is));
976: *contains = (PetscBool)(i >= 0);
977: } else {
978: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
979: }
980: } else { // value is not present
981: *contains = PETSC_FALSE;
982: }
983: }
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
989: When a label is created, it is initialized to -1.
991: Not Collective
993: Input Parameter:
994: . label - a `DMLabel` object
996: Output Parameter:
997: . defaultValue - the default value
999: Level: beginner
1001: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1002: @*/
1003: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
1004: {
1005: PetscFunctionBegin;
1007: *defaultValue = label->defaultValue;
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@
1012: DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
1013: When a label is created, it is initialized to -1.
1015: Not Collective
1017: Input Parameter:
1018: . label - a `DMLabel` object
1020: Output Parameter:
1021: . defaultValue - the default value
1023: Level: beginner
1025: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1026: @*/
1027: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1028: {
1029: PetscFunctionBegin;
1031: label->defaultValue = defaultValue;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: /*@
1036: 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
1037: `DMLabelSetDefaultValue()`)
1039: Not Collective
1041: Input Parameters:
1042: + label - the `DMLabel`
1043: - point - the point
1045: Output Parameter:
1046: . value - The point value, or the default value (-1 by default)
1048: Level: intermediate
1050: Note:
1051: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1052: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1054: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1055: @*/
1056: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1057: {
1058: PetscInt v;
1060: PetscFunctionBeginHot;
1062: PetscAssertPointer(value, 3);
1063: *value = label->defaultValue;
1064: for (v = 0; v < label->numStrata; ++v) {
1065: if (label->validIS[v] || label->readonly) {
1066: IS is;
1067: PetscInt i;
1069: PetscUseTypeMethod(label, getstratumis, v, &is);
1070: PetscCall(ISLocate(label->points[v], point, &i));
1071: PetscCall(ISDestroy(&is));
1072: if (i >= 0) {
1073: *value = label->stratumValues[v];
1074: break;
1075: }
1076: } else {
1077: PetscBool has;
1079: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1080: if (has) {
1081: *value = label->stratumValues[v];
1082: break;
1083: }
1084: }
1085: }
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: /*@
1090: 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
1091: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1093: Not Collective
1095: Input Parameters:
1096: + label - the `DMLabel`
1097: . point - the point
1098: - value - The point value
1100: Level: intermediate
1102: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1103: @*/
1104: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1105: {
1106: PetscInt v;
1108: PetscFunctionBegin;
1110: /* Find label value, add new entry if needed */
1111: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1112: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1113: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1114: /* Set key */
1115: PetscCall(DMLabelMakeInvalid_Private(label, v));
1116: PetscCall(PetscHSetIAdd(label->ht[v], point));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*@
1121: DMLabelClearValue - Clear the value a label assigns to a point
1123: Not Collective
1125: Input Parameters:
1126: + label - the `DMLabel`
1127: . point - the point
1128: - value - The point value
1130: Level: intermediate
1132: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1133: @*/
1134: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1135: {
1136: PetscInt v;
1138: PetscFunctionBegin;
1140: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1141: /* Find label value */
1142: PetscCall(DMLabelLookupStratum(label, value, &v));
1143: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1145: if (label->bt && point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1147: /* Delete key */
1148: PetscCall(DMLabelMakeInvalid_Private(label, v));
1149: PetscCall(PetscHSetIDel(label->ht[v], point));
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: /*@
1154: DMLabelInsertIS - Set all points in the `IS` to a value
1156: Not Collective
1158: Input Parameters:
1159: + label - the `DMLabel`
1160: . is - the point `IS`
1161: - value - The point value
1163: Level: intermediate
1165: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1166: @*/
1167: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1168: {
1169: PetscInt v, n, p;
1170: const PetscInt *points;
1172: PetscFunctionBegin;
1175: /* Find label value, add new entry if needed */
1176: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1177: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1178: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1179: /* Set keys */
1180: PetscCall(DMLabelMakeInvalid_Private(label, v));
1181: PetscCall(ISGetLocalSize(is, &n));
1182: PetscCall(ISGetIndices(is, &points));
1183: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1184: PetscCall(ISRestoreIndices(is, &points));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: /*@
1189: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1191: Not Collective
1193: Input Parameter:
1194: . label - the `DMLabel`
1196: Output Parameter:
1197: . numValues - the number of values
1199: Level: intermediate
1201: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1202: @*/
1203: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1204: {
1205: PetscFunctionBegin;
1207: PetscAssertPointer(numValues, 2);
1208: *numValues = label->numStrata;
1209: PetscFunctionReturn(PETSC_SUCCESS);
1210: }
1212: /*@
1213: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1215: Not Collective
1217: Input Parameter:
1218: . label - the `DMLabel`
1220: Output Parameter:
1221: . values - the value `IS`
1223: Level: intermediate
1225: Notes:
1226: The `values` should be destroyed when no longer needed.
1228: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1230: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1232: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1233: @*/
1234: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1235: {
1236: PetscFunctionBegin;
1238: PetscAssertPointer(values, 2);
1239: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1240: PetscFunctionReturn(PETSC_SUCCESS);
1241: }
1243: /*@
1244: DMLabelGetValueBounds - Return the smallest and largest value in the label
1246: Not Collective
1248: Input Parameter:
1249: . label - the `DMLabel`
1251: Output Parameters:
1252: + minValue - The smallest value
1253: - maxValue - The largest value
1255: Level: intermediate
1257: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1258: @*/
1259: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
1260: {
1261: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;
1263: PetscFunctionBegin;
1265: for (PetscInt v = 0; v < label->numStrata; ++v) {
1266: min = PetscMin(min, label->stratumValues[v]);
1267: max = PetscMax(max, label->stratumValues[v]);
1268: }
1269: if (minValue) {
1270: PetscAssertPointer(minValue, 2);
1271: *minValue = min;
1272: }
1273: if (maxValue) {
1274: PetscAssertPointer(maxValue, 3);
1275: *maxValue = max;
1276: }
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: /*@
1281: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1283: Not Collective
1285: Input Parameter:
1286: . label - the `DMLabel`
1288: Output Parameter:
1289: . values - the value `IS`
1291: Level: intermediate
1293: Notes:
1294: The `values` should be destroyed when no longer needed.
1296: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1298: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetValueISGlobal()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1299: @*/
1300: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1301: {
1302: PetscInt i, j;
1303: PetscInt *valuesArr;
1305: PetscFunctionBegin;
1307: PetscAssertPointer(values, 2);
1308: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1309: for (i = 0, j = 0; i < label->numStrata; i++) {
1310: PetscInt n;
1312: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1313: if (n) valuesArr[j++] = label->stratumValues[i];
1314: }
1315: if (j == label->numStrata) {
1316: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1317: } else {
1318: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1319: }
1320: PetscCall(PetscFree(valuesArr));
1321: PetscFunctionReturn(PETSC_SUCCESS);
1322: }
1324: /*@
1325: DMLabelGetValueISGlobal - Get an `IS` of all values that the `DMlabel` takes across all ranks
1327: Collective
1329: Input Parameter:
1330: + comm - MPI communicator to collect values
1331: . label - the `DMLabel`, may be `NULL` for ranks in `comm` which do not have the corresponding `DMLabel`
1332: - get_nonempty - whether to get nonempty stratum values (akin to `DMLabelGetNonEmptyStratumValuesIS()`)
1334: Output Parameter:
1335: . values - the value `IS`
1337: Level: intermediate
1339: Notes:
1340: The `values` should be destroyed when no longer needed.
1342: This is similar to `DMLabelGetValueIS()` and `DMLabelGetNonEmptyStratumValuesIS()`, but gets the (nonempty) values across all ranks in `comm`.
1344: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1345: @*/
1346: PetscErrorCode DMLabelGetValueISGlobal(MPI_Comm comm, DMLabel label, PetscBool get_nonempty, IS *values)
1347: {
1348: PetscInt num_values_local = 0, num_values_global, minmax_values[2], minmax_values_loc[2] = {PETSC_INT_MAX, PETSC_INT_MIN};
1349: IS is_values = NULL;
1350: const PetscInt *values_local = NULL;
1351: PetscInt *values_global;
1353: PetscFunctionBegin;
1355: if (PetscDefined(USE_DEBUG)) {
1357: IS dummy;
1358: PetscCall(ISCreate(comm, &dummy));
1360: PetscCall(ISDestroy(&dummy));
1361: }
1362: PetscAssertPointer(values, 4);
1364: if (label) {
1365: if (get_nonempty) PetscCall(DMLabelGetNonEmptyStratumValuesIS(label, &is_values));
1366: else PetscCall(DMLabelGetValueIS(label, &is_values));
1367: PetscCall(ISGetIndices(is_values, &values_local));
1368: PetscCall(ISGetLocalSize(is_values, &num_values_local));
1369: }
1370: for (PetscInt i = 0; i < num_values_local; i++) {
1371: minmax_values_loc[0] = PetscMin(minmax_values_loc[0], values_local[i]);
1372: minmax_values_loc[1] = PetscMax(minmax_values_loc[1], values_local[i]);
1373: }
1375: PetscCall(PetscGlobalMinMaxInt(comm, minmax_values_loc, minmax_values));
1376: PetscInt value_range = minmax_values[1] - minmax_values[0] + 1;
1377: PetscBT local_values_bt, global_values_bt;
1379: // Create a "ballot" where each rank marks which values they have into the PetscBT.
1380: // An Allreduce using bitwise-OR over the ranks then communicates which values are owned by a rank in comm
1381: PetscCall(PetscBTCreate(value_range, &local_values_bt));
1382: PetscCall(PetscBTCreate(value_range, &global_values_bt));
1383: for (PetscInt i = 0; i < num_values_local; i++) PetscCall(PetscBTSet(local_values_bt, values_local[i] - minmax_values[0]));
1384: PetscCallMPI(MPIU_Allreduce(local_values_bt, global_values_bt, PetscBTLength(value_range), MPI_CHAR, MPI_BOR, comm));
1385: PetscCall(PetscBTDestroy(&local_values_bt));
1386: {
1387: PetscCount num_values_global_count;
1388: num_values_global_count = PetscBTCountSet(global_values_bt, value_range);
1389: PetscCall(PetscIntCast(num_values_global_count, &num_values_global));
1390: }
1392: PetscCall(PetscMalloc1(num_values_global, &values_global));
1393: for (PetscInt i = 0, a = 0; i < value_range; i++) {
1394: if (PetscBTLookup(global_values_bt, i)) {
1395: values_global[a] = i + minmax_values[0];
1396: a++;
1397: }
1398: }
1399: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_values_global, values_global, PETSC_OWN_POINTER, values));
1401: PetscCall(PetscBTDestroy(&global_values_bt));
1402: if (is_values) {
1403: PetscCall(ISRestoreIndices(is_values, &values_local));
1404: PetscCall(ISDestroy(&is_values));
1405: }
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: /*@
1410: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1412: Not Collective
1414: Input Parameters:
1415: + label - the `DMLabel`
1416: - value - the value
1418: Output Parameter:
1419: . index - the index of value in the list of values
1421: Level: intermediate
1423: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1424: @*/
1425: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1426: {
1427: PetscInt v;
1429: PetscFunctionBegin;
1431: PetscAssertPointer(index, 3);
1432: /* Do not assume they are sorted */
1433: for (v = 0; v < label->numStrata; ++v)
1434: if (label->stratumValues[v] == value) break;
1435: if (v >= label->numStrata) *index = -1;
1436: else *index = v;
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: /*@
1441: DMLabelHasStratum - Determine whether points exist with the given value
1443: Not Collective
1445: Input Parameters:
1446: + label - the `DMLabel`
1447: - value - the stratum value
1449: Output Parameter:
1450: . exists - Flag saying whether points exist
1452: Level: intermediate
1454: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1455: @*/
1456: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1457: {
1458: PetscInt v;
1460: PetscFunctionBegin;
1462: PetscAssertPointer(exists, 3);
1463: PetscCall(DMLabelLookupStratum(label, value, &v));
1464: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: /*@
1469: DMLabelGetStratumSize - Get the size of a stratum
1471: Not Collective
1473: Input Parameters:
1474: + label - the `DMLabel`
1475: - value - the stratum value
1477: Output Parameter:
1478: . size - The number of points in the stratum
1480: Level: intermediate
1482: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1483: @*/
1484: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1485: {
1486: PetscInt v;
1488: PetscFunctionBegin;
1490: PetscAssertPointer(size, 3);
1491: PetscCall(DMLabelLookupStratum(label, value, &v));
1492: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1493: PetscFunctionReturn(PETSC_SUCCESS);
1494: }
1496: /*@
1497: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1499: Not Collective
1501: Input Parameters:
1502: + label - the `DMLabel`
1503: - value - the stratum value
1505: Output Parameters:
1506: + start - the smallest point in the stratum
1507: - end - the largest point in the stratum
1509: Level: intermediate
1511: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1512: @*/
1513: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1514: {
1515: IS is;
1516: PetscInt v, min, max;
1518: PetscFunctionBegin;
1520: if (start) {
1521: PetscAssertPointer(start, 3);
1522: *start = -1;
1523: }
1524: if (end) {
1525: PetscAssertPointer(end, 4);
1526: *end = -1;
1527: }
1528: PetscCall(DMLabelLookupStratum(label, value, &v));
1529: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1530: PetscCall(DMLabelMakeValid_Private(label, v));
1531: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1532: PetscUseTypeMethod(label, getstratumis, v, &is);
1533: PetscCall(ISGetMinMax(is, &min, &max));
1534: PetscCall(ISDestroy(&is));
1535: if (start) *start = min;
1536: if (end) *end = max + 1;
1537: PetscFunctionReturn(PETSC_SUCCESS);
1538: }
1540: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1541: {
1542: PetscFunctionBegin;
1543: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1544: *pointIS = label->points[v];
1545: PetscFunctionReturn(PETSC_SUCCESS);
1546: }
1548: /*@
1549: DMLabelGetStratumIS - Get an `IS` with the stratum points
1551: Not Collective
1553: Input Parameters:
1554: + label - the `DMLabel`
1555: - value - the stratum value
1557: Output Parameter:
1558: . points - The stratum points
1560: Level: intermediate
1562: Notes:
1563: The output `IS` should be destroyed when no longer needed.
1564: Returns `NULL` if the stratum is empty.
1566: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1567: @*/
1568: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1569: {
1570: PetscInt v;
1572: PetscFunctionBegin;
1574: PetscAssertPointer(points, 3);
1575: *points = NULL;
1576: PetscCall(DMLabelLookupStratum(label, value, &v));
1577: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1578: PetscCall(DMLabelMakeValid_Private(label, v));
1579: PetscUseTypeMethod(label, getstratumis, v, points);
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: DMLabelSetStratumIS - Set the stratum points using an `IS`
1586: Not Collective
1588: Input Parameters:
1589: + label - the `DMLabel`
1590: . value - the stratum value
1591: - is - The stratum points
1593: Level: intermediate
1595: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1596: @*/
1597: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1598: {
1599: PetscInt v;
1601: PetscFunctionBegin;
1604: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1605: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1606: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1607: PetscCall(DMLabelClearStratum(label, value));
1608: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1609: PetscCall(PetscObjectReference((PetscObject)is));
1610: PetscCall(ISDestroy(&label->points[v]));
1611: label->points[v] = is;
1612: label->validIS[v] = PETSC_TRUE;
1613: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1614: if (label->bt) {
1615: const PetscInt *points;
1616: PetscInt p;
1618: PetscCall(ISGetIndices(is, &points));
1619: for (p = 0; p < label->stratumSizes[v]; ++p) {
1620: const PetscInt point = points[p];
1622: 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);
1623: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1624: }
1625: }
1626: PetscFunctionReturn(PETSC_SUCCESS);
1627: }
1629: /*@
1630: DMLabelClearStratum - Remove a stratum
1632: Not Collective
1634: Input Parameters:
1635: + label - the `DMLabel`
1636: - value - the stratum value
1638: Level: intermediate
1640: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1641: @*/
1642: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1643: {
1644: PetscInt v;
1646: PetscFunctionBegin;
1648: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1649: PetscCall(DMLabelLookupStratum(label, value, &v));
1650: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1651: if (label->validIS[v]) {
1652: if (label->bt) {
1653: PetscInt i;
1654: const PetscInt *points;
1656: PetscCall(ISGetIndices(label->points[v], &points));
1657: for (i = 0; i < label->stratumSizes[v]; ++i) {
1658: const PetscInt point = points[i];
1660: if (point >= label->pStart && point < label->pEnd) PetscCall(PetscBTClear(label->bt, point - label->pStart));
1661: }
1662: PetscCall(ISRestoreIndices(label->points[v], &points));
1663: }
1664: label->stratumSizes[v] = 0;
1665: PetscCall(ISDestroy(&label->points[v]));
1666: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1667: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1668: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1669: } else {
1670: PetscCall(PetscHSetIClear(label->ht[v]));
1671: }
1672: PetscFunctionReturn(PETSC_SUCCESS);
1673: }
1675: /*@
1676: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1678: Not Collective
1680: Input Parameters:
1681: + label - The `DMLabel`
1682: . value - The label value for all points
1683: . pStart - The first point
1684: - pEnd - A point beyond all marked points
1686: Level: intermediate
1688: Note:
1689: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1691: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1692: @*/
1693: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1694: {
1695: IS pIS;
1697: PetscFunctionBegin;
1698: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1699: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1700: PetscCall(ISDestroy(&pIS));
1701: PetscFunctionReturn(PETSC_SUCCESS);
1702: }
1704: /*@
1705: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1707: Not Collective
1709: Input Parameters:
1710: + label - The `DMLabel`
1711: . value - The label value
1712: - p - A point with this value
1714: Output Parameter:
1715: . 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
1717: Level: intermediate
1719: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1720: @*/
1721: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1722: {
1723: IS pointIS;
1724: PetscInt v;
1726: PetscFunctionBegin;
1728: PetscAssertPointer(index, 4);
1729: *index = -1;
1730: PetscCall(DMLabelLookupStratum(label, value, &v));
1731: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1732: PetscCall(DMLabelMakeValid_Private(label, v));
1733: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1734: PetscCall(ISLocate(pointIS, p, index));
1735: PetscCall(ISDestroy(&pointIS));
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: /*@
1740: DMLabelFilter - Remove all points outside of [`start`, `end`)
1742: Not Collective
1744: Input Parameters:
1745: + label - the `DMLabel`
1746: . start - the first point kept
1747: - end - one more than the last point kept
1749: Level: intermediate
1751: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1752: @*/
1753: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1754: {
1755: PetscInt v;
1757: PetscFunctionBegin;
1759: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1760: PetscCall(DMLabelDestroyIndex(label));
1761: PetscCall(DMLabelMakeAllValid_Private(label));
1762: for (v = 0; v < label->numStrata; ++v) {
1763: PetscCall(ISGeneralFilter(label->points[v], start, end));
1764: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1765: }
1766: PetscCall(DMLabelCreateIndex(label, start, end));
1767: PetscFunctionReturn(PETSC_SUCCESS);
1768: }
1770: /*@
1771: DMLabelPermute - Create a new label with permuted points
1773: Not Collective
1775: Input Parameters:
1776: + label - the `DMLabel`
1777: - permutation - the point permutation
1779: Output Parameter:
1780: . labelNew - the new label containing the permuted points
1782: Level: intermediate
1784: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1785: @*/
1786: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1787: {
1788: const PetscInt *perm;
1789: PetscInt numValues, numPoints, v, q;
1791: PetscFunctionBegin;
1794: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1795: PetscCall(DMLabelMakeAllValid_Private(label));
1796: PetscCall(DMLabelDuplicate(label, labelNew));
1797: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1798: PetscCall(ISGetLocalSize(permutation, &numPoints));
1799: PetscCall(ISGetIndices(permutation, &perm));
1800: for (v = 0; v < numValues; ++v) {
1801: const PetscInt size = (*labelNew)->stratumSizes[v];
1802: const PetscInt *points;
1803: PetscInt *pointsNew;
1805: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1806: PetscCall(PetscCalloc1(size, &pointsNew));
1807: for (q = 0; q < size; ++q) {
1808: const PetscInt point = points[q];
1810: 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);
1811: pointsNew[q] = perm[point];
1812: }
1813: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1814: PetscCall(PetscSortInt(size, pointsNew));
1815: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1816: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1817: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1818: PetscCall(PetscFree(pointsNew));
1819: } else {
1820: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1821: }
1822: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1823: }
1824: PetscCall(ISRestoreIndices(permutation, &perm));
1825: if (label->bt) {
1826: PetscCall(PetscBTDestroy(&label->bt));
1827: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1828: }
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: /*@
1833: DMLabelPermuteValues - Permute the values in a label
1835: Not collective
1837: Input Parameters:
1838: + label - the `DMLabel`
1839: - permutation - the value permutation, permutation[old value] = new value
1841: Output Parameter:
1842: . label - the `DMLabel` now with permuted values
1844: Note:
1845: The modification is done in-place
1847: Level: intermediate
1849: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1850: @*/
1851: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1852: {
1853: PetscInt Nv, Np;
1855: PetscFunctionBegin;
1858: PetscCall(DMLabelGetNumValues(label, &Nv));
1859: PetscCall(ISGetLocalSize(permutation, &Np));
1860: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1861: if (PetscDefined(USE_DEBUG)) {
1862: PetscBool flg;
1863: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1864: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1865: }
1866: PetscCall(DMLabelRewriteValues(label, permutation));
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: /*@
1871: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1873: Not collective
1875: Input Parameters:
1876: + label - the `DMLabel`
1877: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1879: Output Parameter:
1880: . label - the `DMLabel` now with permuted values
1882: Note:
1883: The modification is done in-place
1885: Level: intermediate
1887: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1888: @*/
1889: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1890: {
1891: const PetscInt *perm;
1892: PetscInt Nv, Np;
1894: PetscFunctionBegin;
1897: PetscCall(DMLabelMakeAllValid_Private(label));
1898: PetscCall(DMLabelGetNumValues(label, &Nv));
1899: PetscCall(ISGetLocalSize(permutation, &Np));
1900: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1901: PetscCall(ISGetIndices(permutation, &perm));
1902: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1903: PetscCall(ISRestoreIndices(permutation, &perm));
1904: PetscFunctionReturn(PETSC_SUCCESS);
1905: }
1907: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1908: {
1909: MPI_Comm comm;
1910: PetscInt s, l, nroots, nleaves, offset, size;
1911: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1912: PetscSection rootSection;
1913: PetscSF labelSF;
1915: PetscFunctionBegin;
1916: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1917: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1918: /* Build a section of stratum values per point, generate the according SF
1919: and distribute point-wise stratum values to leaves. */
1920: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1921: PetscCall(PetscSectionCreate(comm, &rootSection));
1922: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1923: if (label) {
1924: for (s = 0; s < label->numStrata; ++s) {
1925: const PetscInt *points;
1927: PetscCall(ISGetIndices(label->points[s], &points));
1928: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1929: PetscCall(ISRestoreIndices(label->points[s], &points));
1930: }
1931: }
1932: PetscCall(PetscSectionSetUp(rootSection));
1933: /* Create a point-wise array of stratum values */
1934: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1935: PetscCall(PetscMalloc1(size, &rootStrata));
1936: PetscCall(PetscCalloc1(nroots, &rootIdx));
1937: if (label) {
1938: for (s = 0; s < label->numStrata; ++s) {
1939: const PetscInt *points;
1941: PetscCall(ISGetIndices(label->points[s], &points));
1942: for (l = 0; l < label->stratumSizes[s]; l++) {
1943: const PetscInt p = points[l];
1944: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1945: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1946: }
1947: PetscCall(ISRestoreIndices(label->points[s], &points));
1948: }
1949: }
1950: /* Build SF that maps label points to remote processes */
1951: PetscCall(PetscSectionCreate(comm, leafSection));
1952: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1953: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1954: PetscCall(PetscFree(remoteOffsets));
1955: /* Send the strata for each point over the derived SF */
1956: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1957: PetscCall(PetscMalloc1(size, leafStrata));
1958: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1959: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1960: /* Clean up */
1961: PetscCall(PetscFree(rootStrata));
1962: PetscCall(PetscFree(rootIdx));
1963: PetscCall(PetscSectionDestroy(&rootSection));
1964: PetscCall(PetscSFDestroy(&labelSF));
1965: PetscFunctionReturn(PETSC_SUCCESS);
1966: }
1968: /*@
1969: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1971: Collective
1973: Input Parameters:
1974: + label - the `DMLabel`
1975: - sf - the map from old to new distribution
1977: Output Parameter:
1978: . labelNew - the new redistributed label
1980: Level: intermediate
1982: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1983: @*/
1984: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1985: {
1986: MPI_Comm comm;
1987: PetscSection leafSection;
1988: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1989: PetscInt *leafStrata, *strataIdx;
1990: PetscInt **points;
1991: const char *lname = NULL;
1992: char *name;
1993: PetscMPIInt nameSize;
1994: PetscHSetI stratumHash;
1995: size_t len = 0;
1996: PetscMPIInt rank;
1998: PetscFunctionBegin;
2000: if (label) {
2002: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2003: PetscCall(DMLabelMakeAllValid_Private(label));
2004: }
2005: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2006: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2007: /* Bcast name */
2008: if (rank == 0) {
2009: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2010: PetscCall(PetscStrlen(lname, &len));
2011: }
2012: PetscCall(PetscMPIIntCast(len, &nameSize));
2013: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2014: PetscCall(PetscMalloc1(nameSize + 1, &name));
2015: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2016: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2017: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2018: PetscCall(PetscFree(name));
2019: /* Bcast defaultValue */
2020: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
2021: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
2022: /* Distribute stratum values over the SF and get the point mapping on the receiver */
2023: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
2024: /* Determine received stratum values and initialise new label*/
2025: PetscCall(PetscHSetICreate(&stratumHash));
2026: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
2027: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
2028: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
2029: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
2030: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
2031: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
2032: /* Turn leafStrata into indices rather than stratum values */
2033: offset = 0;
2034: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
2035: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
2036: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
2037: for (p = 0; p < size; ++p) {
2038: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2039: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
2040: leafStrata[p] = s;
2041: break;
2042: }
2043: }
2044: }
2045: /* Rebuild the point strata on the receiver */
2046: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
2047: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2048: for (p = pStart; p < pEnd; p++) {
2049: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2050: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2051: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
2052: }
2053: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
2054: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
2055: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
2056: for (s = 0; s < (*labelNew)->numStrata; ++s) {
2057: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
2058: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
2059: }
2060: /* Insert points into new strata */
2061: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
2062: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2063: for (p = pStart; p < pEnd; p++) {
2064: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2065: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
2066: for (s = 0; s < dof; s++) {
2067: stratum = leafStrata[offset + s];
2068: points[stratum][strataIdx[stratum]++] = p;
2069: }
2070: }
2071: for (s = 0; s < (*labelNew)->numStrata; s++) {
2072: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
2073: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
2074: }
2075: PetscCall(PetscFree(points));
2076: PetscCall(PetscHSetIDestroy(&stratumHash));
2077: PetscCall(PetscFree(leafStrata));
2078: PetscCall(PetscFree(strataIdx));
2079: PetscCall(PetscSectionDestroy(&leafSection));
2080: PetscFunctionReturn(PETSC_SUCCESS);
2081: }
2083: /*@
2084: DMLabelGather - Gather all label values from leafs into roots
2086: Collective
2088: Input Parameters:
2089: + label - the `DMLabel`
2090: - sf - the `PetscSF` communication map
2092: Output Parameter:
2093: . labelNew - the new `DMLabel` with localised leaf values
2095: Level: developer
2097: Note:
2098: This is the inverse operation to `DMLabelDistribute()`.
2100: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2101: @*/
2102: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
2103: {
2104: MPI_Comm comm;
2105: PetscSection rootSection;
2106: PetscSF sfLabel;
2107: PetscSFNode *rootPoints, *leafPoints;
2108: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2109: const PetscInt *rootDegree, *ilocal;
2110: PetscInt *rootStrata;
2111: const char *lname;
2112: char *name;
2113: PetscMPIInt nameSize;
2114: size_t len = 0;
2115: PetscMPIInt rank, size;
2117: PetscFunctionBegin;
2120: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2121: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2122: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2123: PetscCallMPI(MPI_Comm_size(comm, &size));
2124: /* Bcast name */
2125: if (rank == 0) {
2126: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2127: PetscCall(PetscStrlen(lname, &len));
2128: }
2129: PetscCall(PetscMPIIntCast(len, &nameSize));
2130: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2131: PetscCall(PetscMalloc1(nameSize + 1, &name));
2132: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2133: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2134: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2135: PetscCall(PetscFree(name));
2136: /* Gather rank/index pairs of leaves into local roots to build
2137: an inverse, multi-rooted SF. Note that this ignores local leaf
2138: indexing due to the use of the multiSF in PetscSFGather. */
2139: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2140: PetscCall(PetscMalloc1(nroots, &leafPoints));
2141: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2142: for (p = 0; p < nleaves; p++) {
2143: PetscInt ilp = ilocal ? ilocal[p] : p;
2145: leafPoints[ilp].index = ilp;
2146: leafPoints[ilp].rank = rank;
2147: }
2148: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2149: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2150: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2151: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2152: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2153: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2154: PetscCall(PetscSFCreate(comm, &sfLabel));
2155: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2156: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2157: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2158: /* Rebuild the point strata on the receiver */
2159: for (p = 0, idx = 0; p < nroots; p++) {
2160: for (d = 0; d < rootDegree[p]; d++) {
2161: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2162: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2163: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2164: }
2165: idx += rootDegree[p];
2166: }
2167: PetscCall(PetscFree(leafPoints));
2168: PetscCall(PetscFree(rootStrata));
2169: PetscCall(PetscSectionDestroy(&rootSection));
2170: PetscCall(PetscSFDestroy(&sfLabel));
2171: PetscFunctionReturn(PETSC_SUCCESS);
2172: }
2174: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2175: {
2176: const PetscInt *degree;
2177: const PetscInt *points;
2178: PetscInt Nr, r, Nl, l, val, defVal;
2180: PetscFunctionBegin;
2181: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2182: /* Add in leaves */
2183: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2184: for (l = 0; l < Nl; ++l) {
2185: PetscCall(DMLabelGetValue(label, points[l], &val));
2186: if (val != defVal) valArray[points[l]] = val;
2187: }
2188: /* Add in shared roots */
2189: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2190: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2191: for (r = 0; r < Nr; ++r) {
2192: if (degree[r]) {
2193: PetscCall(DMLabelGetValue(label, r, &val));
2194: if (val != defVal) valArray[r] = val;
2195: }
2196: }
2197: PetscFunctionReturn(PETSC_SUCCESS);
2198: }
2200: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), PetscCtx ctx)
2201: {
2202: const PetscInt *degree;
2203: const PetscInt *points;
2204: PetscInt Nr, r, Nl, l, val, defVal;
2206: PetscFunctionBegin;
2207: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2208: /* Read out leaves */
2209: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2210: for (l = 0; l < Nl; ++l) {
2211: const PetscInt p = points[l];
2212: const PetscInt cval = valArray[p];
2214: if (cval != defVal) {
2215: PetscCall(DMLabelGetValue(label, p, &val));
2216: if (val == defVal) {
2217: PetscCall(DMLabelSetValue(label, p, cval));
2218: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2219: }
2220: }
2221: }
2222: /* Read out shared roots */
2223: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2224: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2225: for (r = 0; r < Nr; ++r) {
2226: if (degree[r]) {
2227: const PetscInt cval = valArray[r];
2229: if (cval != defVal) {
2230: PetscCall(DMLabelGetValue(label, r, &val));
2231: if (val == defVal) {
2232: PetscCall(DMLabelSetValue(label, r, cval));
2233: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2234: }
2235: }
2236: }
2237: }
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: }
2241: /*@
2242: DMLabelPropagateBegin - Setup a cycle of label propagation
2244: Collective
2246: Input Parameters:
2247: + label - The `DMLabel` to propagate across processes
2248: - sf - The `PetscSF` describing parallel layout of the label points
2250: Level: intermediate
2252: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2253: @*/
2254: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2255: {
2256: PetscInt Nr, r, defVal;
2257: PetscMPIInt size;
2259: PetscFunctionBegin;
2260: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2261: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2262: if (size > 1) {
2263: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2264: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2265: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2266: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2267: }
2268: PetscFunctionReturn(PETSC_SUCCESS);
2269: }
2271: /*@
2272: DMLabelPropagateEnd - Tear down a cycle of label propagation
2274: Collective
2276: Input Parameters:
2277: + label - The `DMLabel` to propagate across processes
2278: - pointSF - The `PetscSF` describing parallel layout of the label points
2280: Level: intermediate
2282: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2283: @*/
2284: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2285: {
2286: PetscFunctionBegin;
2287: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2288: PetscCall(PetscFree(label->propArray));
2289: label->propArray = NULL;
2290: PetscFunctionReturn(PETSC_SUCCESS);
2291: }
2293: /*@C
2294: DMLabelPropagatePush - Tear down a cycle of label propagation
2296: Collective
2298: Input Parameters:
2299: + label - The `DMLabel` to propagate across processes
2300: . pointSF - The `PetscSF` describing parallel layout of the label points
2301: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2302: - ctx - An optional user context for the callback, or `NULL`
2304: Calling sequence of `markPoint`:
2305: + label - The `DMLabel`
2306: . p - The point being marked
2307: . val - The label value for `p`
2308: - ctx - An optional user context
2310: Level: intermediate
2312: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2313: @*/
2314: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, PetscCtx ctx), PetscCtx ctx)
2315: {
2316: PetscInt *valArray = label->propArray, Nr;
2317: PetscMPIInt size;
2319: PetscFunctionBegin;
2320: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2321: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2322: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2323: if (size > 1 && Nr >= 0) {
2324: /* Communicate marked edges
2325: The current implementation allocates an array the size of the number of root. We put the label values into the
2326: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2328: TODO: We could use in-place communication with a different SF
2329: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2330: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2332: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2333: 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
2334: edge to the queue.
2335: */
2336: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2337: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2338: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2339: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2340: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2341: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2342: }
2343: PetscFunctionReturn(PETSC_SUCCESS);
2344: }
2346: /*@
2347: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2349: Not Collective
2351: Input Parameter:
2352: . label - the `DMLabel`
2354: Output Parameters:
2355: + section - the section giving offsets for each stratum
2356: - is - An `IS` containing all the label points
2358: Level: developer
2360: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2361: @*/
2362: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2363: {
2364: IS vIS;
2365: const PetscInt *values;
2366: PetscInt *points;
2367: PetscInt nV, vS = 0, vE = 0, v, N;
2369: PetscFunctionBegin;
2371: PetscCall(DMLabelGetNumValues(label, &nV));
2372: PetscCall(DMLabelGetValueIS(label, &vIS));
2373: PetscCall(ISGetIndices(vIS, &values));
2374: if (nV) {
2375: vS = values[0];
2376: vE = values[0] + 1;
2377: }
2378: for (v = 1; v < nV; ++v) {
2379: vS = PetscMin(vS, values[v]);
2380: vE = PetscMax(vE, values[v] + 1);
2381: }
2382: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2383: PetscCall(PetscSectionSetChart(*section, vS, vE));
2384: for (v = 0; v < nV; ++v) {
2385: PetscInt n;
2387: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2388: PetscCall(PetscSectionSetDof(*section, values[v], n));
2389: }
2390: PetscCall(PetscSectionSetUp(*section));
2391: PetscCall(PetscSectionGetStorageSize(*section, &N));
2392: PetscCall(PetscMalloc1(N, &points));
2393: for (v = 0; v < nV; ++v) {
2394: IS is;
2395: const PetscInt *spoints;
2396: PetscInt dof, off, p;
2398: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2399: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2400: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2401: PetscCall(ISGetIndices(is, &spoints));
2402: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2403: PetscCall(ISRestoreIndices(is, &spoints));
2404: PetscCall(ISDestroy(&is));
2405: }
2406: PetscCall(ISRestoreIndices(vIS, &values));
2407: PetscCall(ISDestroy(&vIS));
2408: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2409: PetscFunctionReturn(PETSC_SUCCESS);
2410: }
2412: /*@C
2413: DMLabelRegister - Adds a new label component implementation
2415: Not Collective
2417: Input Parameters:
2418: + name - The name of a new user-defined creation routine
2419: - create_func - The creation routine itself
2421: Notes:
2422: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2424: Example Usage:
2425: .vb
2426: DMLabelRegister("my_label", MyLabelCreate);
2427: .ve
2429: Then, your label type can be chosen with the procedural interface via
2430: .vb
2431: DMLabelCreate(MPI_Comm, DMLabel *);
2432: DMLabelSetType(DMLabel, "my_label");
2433: .ve
2434: or at runtime via the option
2435: .vb
2436: -dm_label_type my_label
2437: .ve
2439: Level: advanced
2441: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2442: @*/
2443: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2444: {
2445: PetscFunctionBegin;
2446: PetscCall(DMInitializePackage());
2447: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2448: PetscFunctionReturn(PETSC_SUCCESS);
2449: }
2451: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2452: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2454: /*@C
2455: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2457: Not Collective
2459: Level: advanced
2461: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2462: @*/
2463: PetscErrorCode DMLabelRegisterAll(void)
2464: {
2465: PetscFunctionBegin;
2466: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2467: DMLabelRegisterAllCalled = PETSC_TRUE;
2469: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2470: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2471: PetscFunctionReturn(PETSC_SUCCESS);
2472: }
2474: /*@C
2475: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2477: Level: developer
2479: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2480: @*/
2481: PetscErrorCode DMLabelRegisterDestroy(void)
2482: {
2483: PetscFunctionBegin;
2484: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2485: DMLabelRegisterAllCalled = PETSC_FALSE;
2486: PetscFunctionReturn(PETSC_SUCCESS);
2487: }
2489: /*@
2490: DMLabelSetType - Sets the particular implementation for a label.
2492: Collective
2494: Input Parameters:
2495: + label - The label
2496: - method - The name of the label type
2498: Options Database Key:
2499: . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType`
2501: Level: intermediate
2503: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2504: @*/
2505: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2506: {
2507: PetscErrorCode (*r)(DMLabel);
2508: PetscBool match;
2510: PetscFunctionBegin;
2512: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2513: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2515: PetscCall(DMLabelRegisterAll());
2516: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2517: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2519: PetscTryTypeMethod(label, destroy);
2520: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2521: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2522: PetscCall((*r)(label));
2523: PetscFunctionReturn(PETSC_SUCCESS);
2524: }
2526: /*@
2527: DMLabelGetType - Gets the type name (as a string) from the label.
2529: Not Collective
2531: Input Parameter:
2532: . label - The `DMLabel`
2534: Output Parameter:
2535: . type - The `DMLabel` type name
2537: Level: intermediate
2539: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2540: @*/
2541: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2542: {
2543: PetscFunctionBegin;
2545: PetscAssertPointer(type, 2);
2546: PetscCall(DMLabelRegisterAll());
2547: *type = ((PetscObject)label)->type_name;
2548: PetscFunctionReturn(PETSC_SUCCESS);
2549: }
2551: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2552: {
2553: PetscFunctionBegin;
2554: label->ops->view = DMLabelView_Concrete;
2555: label->ops->setup = NULL;
2556: label->ops->duplicate = DMLabelDuplicate_Concrete;
2557: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2558: PetscFunctionReturn(PETSC_SUCCESS);
2559: }
2561: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2562: {
2563: PetscFunctionBegin;
2565: PetscCall(DMLabelInitialize_Concrete(label));
2566: PetscFunctionReturn(PETSC_SUCCESS);
2567: }
2569: /*@
2570: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2571: the local section and an `PetscSF` describing the section point overlap.
2573: Collective
2575: Input Parameters:
2576: + s - The `PetscSection` for the local field layout
2577: . sf - The `PetscSF` describing parallel layout of the section points
2578: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2579: . label - The label specifying the points
2580: - labelValue - The label stratum specifying the points
2582: Output Parameter:
2583: . gsection - The `PetscSection` for the global field layout
2585: Level: developer
2587: Note:
2588: This gives negative sizes and offsets to points not owned by this process
2590: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2591: @*/
2592: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2593: {
2594: PetscInt *neg = NULL, *tmpOff = NULL;
2595: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2597: PetscFunctionBegin;
2601: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2602: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2603: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2604: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2605: if (nroots >= 0) {
2606: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2607: PetscCall(PetscCalloc1(nroots, &neg));
2608: if (nroots > pEnd - pStart) {
2609: PetscCall(PetscCalloc1(nroots, &tmpOff));
2610: } else {
2611: tmpOff = &(*gsection)->atlasDof[-pStart];
2612: }
2613: }
2614: /* Mark ghost points with negative dof */
2615: for (p = pStart; p < pEnd; ++p) {
2616: PetscInt value;
2618: PetscCall(DMLabelGetValue(label, p, &value));
2619: if (value != labelValue) continue;
2620: PetscCall(PetscSectionGetDof(s, p, &dof));
2621: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2622: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2623: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2624: if (neg) neg[p] = -(dof + 1);
2625: }
2626: PetscCall(PetscSectionSetUpBC(*gsection));
2627: if (nroots >= 0) {
2628: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2629: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2630: if (nroots > pEnd - pStart) {
2631: for (p = pStart; p < pEnd; ++p) {
2632: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2633: }
2634: }
2635: }
2636: /* Calculate new sizes, get process offset, and calculate point offsets */
2637: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2638: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2639: (*gsection)->atlasOff[p] = off;
2640: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2641: }
2642: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2643: globalOff -= off;
2644: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2645: (*gsection)->atlasOff[p] += globalOff;
2646: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2647: }
2648: /* Put in negative offsets for ghost points */
2649: if (nroots >= 0) {
2650: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2651: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2652: if (nroots > pEnd - pStart) {
2653: for (p = pStart; p < pEnd; ++p) {
2654: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2655: }
2656: }
2657: }
2658: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2659: PetscCall(PetscFree(neg));
2660: PetscFunctionReturn(PETSC_SUCCESS);
2661: }
2663: typedef struct _n_PetscSectionSym_Label {
2664: DMLabel label;
2665: PetscCopyMode *modes;
2666: PetscInt *sizes;
2667: const PetscInt ***perms;
2668: const PetscScalar ***rots;
2669: PetscInt (*minMaxOrients)[2];
2670: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2671: } PetscSectionSym_Label;
2673: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2674: {
2675: PetscInt i, j;
2676: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2678: PetscFunctionBegin;
2679: for (i = 0; i <= sl->numStrata; i++) {
2680: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2681: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2682: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2683: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2684: }
2685: if (sl->perms[i]) {
2686: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2688: PetscCall(PetscFree(perms));
2689: }
2690: if (sl->rots[i]) {
2691: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2693: PetscCall(PetscFree(rots));
2694: }
2695: }
2696: }
2697: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2698: PetscCall(DMLabelDestroy(&sl->label));
2699: sl->numStrata = 0;
2700: PetscFunctionReturn(PETSC_SUCCESS);
2701: }
2703: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2704: {
2705: PetscFunctionBegin;
2706: PetscCall(PetscSectionSymLabelReset(sym));
2707: PetscCall(PetscFree(sym->data));
2708: PetscFunctionReturn(PETSC_SUCCESS);
2709: }
2711: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2712: {
2713: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2714: PetscBool isAscii;
2715: DMLabel label = sl->label;
2716: const char *name;
2718: PetscFunctionBegin;
2719: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2720: if (isAscii) {
2721: PetscInt i, j, k;
2722: PetscViewerFormat format;
2724: PetscCall(PetscViewerGetFormat(viewer, &format));
2725: if (label) {
2726: PetscCall(PetscViewerGetFormat(viewer, &format));
2727: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2728: PetscCall(PetscViewerASCIIPushTab(viewer));
2729: PetscCall(DMLabelView(label, viewer));
2730: PetscCall(PetscViewerASCIIPopTab(viewer));
2731: } else {
2732: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2733: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2734: }
2735: } else {
2736: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2737: }
2738: PetscCall(PetscViewerASCIIPushTab(viewer));
2739: for (i = 0; i <= sl->numStrata; i++) {
2740: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2742: if (!(sl->perms[i] || sl->rots[i])) {
2743: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2744: } else {
2745: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2746: PetscCall(PetscViewerASCIIPushTab(viewer));
2747: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2748: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2749: PetscCall(PetscViewerASCIIPushTab(viewer));
2750: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2751: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2752: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2753: } else {
2754: PetscInt tab;
2756: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2757: PetscCall(PetscViewerASCIIPushTab(viewer));
2758: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2759: if (sl->perms[i] && sl->perms[i][j]) {
2760: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2761: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2762: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2763: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2764: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2765: }
2766: if (sl->rots[i] && sl->rots[i][j]) {
2767: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2768: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2769: #if defined(PETSC_USE_COMPLEX)
2770: 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])));
2771: #else
2772: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2773: #endif
2774: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2775: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2776: }
2777: PetscCall(PetscViewerASCIIPopTab(viewer));
2778: }
2779: }
2780: PetscCall(PetscViewerASCIIPopTab(viewer));
2781: }
2782: PetscCall(PetscViewerASCIIPopTab(viewer));
2783: }
2784: }
2785: PetscCall(PetscViewerASCIIPopTab(viewer));
2786: }
2787: PetscFunctionReturn(PETSC_SUCCESS);
2788: }
2790: /*@
2791: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2793: Logically
2795: Input Parameters:
2796: + sym - the section symmetries
2797: - label - the `DMLabel` describing the types of points
2799: Level: developer:
2801: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2802: @*/
2803: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2804: {
2805: PetscSectionSym_Label *sl;
2807: PetscFunctionBegin;
2809: sl = (PetscSectionSym_Label *)sym->data;
2810: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2811: if (label) {
2812: sl->label = label;
2813: PetscCall(PetscObjectReference((PetscObject)label));
2814: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2815: 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));
2816: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2817: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2818: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2819: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2820: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2821: }
2822: PetscFunctionReturn(PETSC_SUCCESS);
2823: }
2825: /*@C
2826: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2828: Logically Collective
2830: Input Parameters:
2831: + sym - the section symmetries
2832: - stratum - the stratum value in the label that we are assigning symmetries for
2834: Output Parameters:
2835: + size - the number of dofs for points in the `stratum` of the label
2836: . minOrient - the smallest orientation for a point in this `stratum`
2837: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2838: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2839: - 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
2841: Level: developer
2843: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2844: @*/
2845: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2846: {
2847: PetscSectionSym_Label *sl;
2848: const char *name;
2849: PetscInt i;
2851: PetscFunctionBegin;
2853: sl = (PetscSectionSym_Label *)sym->data;
2854: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2855: for (i = 0; i <= sl->numStrata; i++) {
2856: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2858: if (stratum == value) break;
2859: }
2860: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2861: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2862: if (size) {
2863: PetscAssertPointer(size, 3);
2864: *size = sl->sizes[i];
2865: }
2866: if (minOrient) {
2867: PetscAssertPointer(minOrient, 4);
2868: *minOrient = sl->minMaxOrients[i][0];
2869: }
2870: if (maxOrient) {
2871: PetscAssertPointer(maxOrient, 5);
2872: *maxOrient = sl->minMaxOrients[i][1];
2873: }
2874: if (perms) {
2875: PetscAssertPointer(perms, 6);
2876: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2877: }
2878: if (rots) {
2879: PetscAssertPointer(rots, 7);
2880: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2881: }
2882: PetscFunctionReturn(PETSC_SUCCESS);
2883: }
2885: /*@C
2886: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2888: Logically
2890: Input Parameters:
2891: + sym - the section symmetries
2892: . stratum - the stratum value in the label that we are assigning symmetries for
2893: . size - the number of dofs for points in the `stratum` of the label
2894: . minOrient - the smallest orientation for a point in this `stratum`
2895: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2896: . mode - how `sym` should copy the `perms` and `rots` arrays
2897: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2898: - 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
2900: Level: developer
2902: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2903: @*/
2904: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2905: {
2906: PetscSectionSym_Label *sl;
2907: const char *name;
2908: PetscInt i, j, k;
2910: PetscFunctionBegin;
2912: sl = (PetscSectionSym_Label *)sym->data;
2913: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2914: for (i = 0; i <= sl->numStrata; i++) {
2915: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2917: if (stratum == value) break;
2918: }
2919: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2920: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2921: sl->sizes[i] = size;
2922: sl->modes[i] = mode;
2923: sl->minMaxOrients[i][0] = minOrient;
2924: sl->minMaxOrients[i][1] = maxOrient;
2925: if (mode == PETSC_COPY_VALUES) {
2926: if (perms) {
2927: PetscInt **ownPerms;
2929: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2930: for (j = 0; j < maxOrient - minOrient; j++) {
2931: if (perms[j]) {
2932: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2933: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2934: }
2935: }
2936: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2937: }
2938: if (rots) {
2939: PetscScalar **ownRots;
2941: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2942: for (j = 0; j < maxOrient - minOrient; j++) {
2943: if (rots[j]) {
2944: PetscCall(PetscMalloc1(size, &ownRots[j]));
2945: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2946: }
2947: }
2948: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2949: }
2950: } else {
2951: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2952: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2953: }
2954: PetscFunctionReturn(PETSC_SUCCESS);
2955: }
2957: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2958: {
2959: PetscInt i, j, numStrata;
2960: PetscSectionSym_Label *sl;
2961: DMLabel label;
2963: PetscFunctionBegin;
2964: sl = (PetscSectionSym_Label *)sym->data;
2965: numStrata = sl->numStrata;
2966: label = sl->label;
2967: for (i = 0; i < numPoints; i++) {
2968: PetscInt point = points[2 * i];
2969: PetscInt ornt = points[2 * i + 1];
2971: for (j = 0; j < numStrata; j++) {
2972: if (label->validIS[j]) {
2973: PetscInt k;
2975: PetscCall(ISLocate(label->points[j], point, &k));
2976: if (k >= 0) break;
2977: } else {
2978: PetscBool has;
2980: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2981: if (has) break;
2982: }
2983: }
2984: 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],
2985: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2986: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2987: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2988: }
2989: PetscFunctionReturn(PETSC_SUCCESS);
2990: }
2992: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2993: {
2994: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2995: IS valIS;
2996: const PetscInt *values;
2997: PetscInt Nv, v;
2999: PetscFunctionBegin;
3000: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
3001: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
3002: PetscCall(ISGetIndices(valIS, &values));
3003: for (v = 0; v < Nv; ++v) {
3004: const PetscInt val = values[v];
3005: PetscInt size, minOrient, maxOrient;
3006: const PetscInt **perms;
3007: const PetscScalar **rots;
3009: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
3010: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
3011: }
3012: PetscCall(ISDestroy(&valIS));
3013: PetscFunctionReturn(PETSC_SUCCESS);
3014: }
3016: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3017: {
3018: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
3019: DMLabel dlabel;
3021: PetscFunctionBegin;
3022: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
3023: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
3024: PetscCall(DMLabelDestroy(&dlabel));
3025: PetscCall(PetscSectionSymCopy(sym, *dsym));
3026: PetscFunctionReturn(PETSC_SUCCESS);
3027: }
3029: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
3030: {
3031: PetscSectionSym_Label *sl;
3033: PetscFunctionBegin;
3034: PetscCall(PetscNew(&sl));
3035: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
3036: sym->ops->distribute = PetscSectionSymDistribute_Label;
3037: sym->ops->copy = PetscSectionSymCopy_Label;
3038: sym->ops->view = PetscSectionSymView_Label;
3039: sym->ops->destroy = PetscSectionSymDestroy_Label;
3040: sym->data = (void *)sl;
3041: PetscFunctionReturn(PETSC_SUCCESS);
3042: }
3044: /*@
3045: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
3047: Collective
3049: Input Parameters:
3050: + comm - the MPI communicator for the new symmetry
3051: - label - the label defining the strata
3053: Output Parameter:
3054: . sym - the section symmetries
3056: Level: developer
3058: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3059: @*/
3060: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
3061: {
3062: PetscFunctionBegin;
3063: PetscCall(DMInitializePackage());
3064: PetscCall(PetscSectionSymCreate(comm, sym));
3065: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
3066: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
3067: PetscFunctionReturn(PETSC_SUCCESS);
3068: }