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: /*@C
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;
34: PetscCall(DMInitializePackage());
36: PetscCall(PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView));
38: (*label)->numStrata = 0;
39: (*label)->defaultValue = -1;
40: (*label)->stratumValues = NULL;
41: (*label)->validIS = NULL;
42: (*label)->stratumSizes = NULL;
43: (*label)->points = NULL;
44: (*label)->ht = NULL;
45: (*label)->pStart = -1;
46: (*label)->pEnd = -1;
47: (*label)->bt = NULL;
48: PetscCall(PetscHMapICreate(&(*label)->hmap));
49: PetscCall(PetscObjectSetName((PetscObject)*label, name));
50: PetscCall(DMLabelSetType(*label, DMLABELCONCRETE));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@C
55: DMLabelSetUp - SetUp a `DMLabel` object
57: Collective
59: Input parameters:
60: . label - The `DMLabel`
62: Level: intermediate
64: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
65: @*/
66: PetscErrorCode DMLabelSetUp(DMLabel label)
67: {
68: PetscFunctionBegin;
70: PetscTryTypeMethod(label, setup);
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*
75: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
77: Not collective
79: Input parameter:
80: + label - The `DMLabel`
81: - v - The stratum value
83: Output parameter:
84: . label - The `DMLabel` with stratum in sorted list format
86: Level: developer
88: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
89: */
90: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
91: {
92: IS is;
93: PetscInt off = 0, *pointArray, p;
95: if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) return PETSC_SUCCESS;
96: PetscFunctionBegin;
97: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
98: PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
99: PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
100: PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
101: PetscCall(PetscHSetIClear(label->ht[v]));
102: PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
103: if (label->bt) {
104: for (p = 0; p < label->stratumSizes[v]; ++p) {
105: const PetscInt point = pointArray[p];
106: 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);
107: PetscCall(PetscBTSet(label->bt, point - label->pStart));
108: }
109: }
110: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
111: PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
112: PetscCall(PetscFree(pointArray));
113: } else {
114: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
115: }
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: if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) return PETSC_SUCCESS;
169: PetscFunctionBegin;
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;
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: PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
448: {
449: PetscBool iascii;
451: PetscFunctionBegin;
452: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
453: if (iascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@C
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: DMLabelReset - Destroys internal data structures in a `DMLabel`
484: Not Collective
486: Input Parameter:
487: . label - The `DMLabel`
489: Level: beginner
491: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`, `DMLabelCreate()`
492: @*/
493: PetscErrorCode DMLabelReset(DMLabel label)
494: {
495: PetscInt v;
497: PetscFunctionBegin;
499: for (v = 0; v < label->numStrata; ++v) {
500: if (label->ht[v]) PetscCall(PetscHSetIDestroy(&label->ht[v]));
501: PetscCall(ISDestroy(&label->points[v]));
502: }
503: label->numStrata = 0;
504: PetscCall(PetscFree(label->stratumValues));
505: PetscCall(PetscFree(label->stratumSizes));
506: PetscCall(PetscFree(label->ht));
507: PetscCall(PetscFree(label->points));
508: PetscCall(PetscFree(label->validIS));
509: label->stratumValues = NULL;
510: label->stratumSizes = NULL;
511: label->ht = NULL;
512: label->points = NULL;
513: label->validIS = NULL;
514: PetscCall(PetscHMapIReset(label->hmap));
515: label->pStart = -1;
516: label->pEnd = -1;
517: PetscCall(PetscBTDestroy(&label->bt));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@
522: DMLabelDestroy - Destroys a `DMLabel`
524: Collective
526: Input Parameter:
527: . label - The `DMLabel`
529: Level: beginner
531: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
532: @*/
533: PetscErrorCode DMLabelDestroy(DMLabel *label)
534: {
535: PetscFunctionBegin;
536: if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
538: if (--((PetscObject)(*label))->refct > 0) {
539: *label = NULL;
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
542: PetscCall(DMLabelReset(*label));
543: PetscCall(PetscHMapIDestroy(&(*label)->hmap));
544: PetscCall(PetscHeaderDestroy(label));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
549: {
550: PetscFunctionBegin;
551: for (PetscInt v = 0; v < label->numStrata; ++v) {
552: PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
553: PetscCall(PetscObjectReference((PetscObject)(label->points[v])));
554: (*labelnew)->points[v] = label->points[v];
555: }
556: PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
557: PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /*@
562: DMLabelDuplicate - Duplicates a `DMLabel`
564: Collective
566: Input Parameter:
567: . label - The `DMLabel`
569: Output Parameter:
570: . labelnew - new label
572: Level: intermediate
574: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
575: @*/
576: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
577: {
578: const char *name;
580: PetscFunctionBegin;
582: PetscCall(DMLabelMakeAllValid_Private(label));
583: PetscCall(PetscObjectGetName((PetscObject)label, &name));
584: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));
586: (*labelnew)->numStrata = label->numStrata;
587: (*labelnew)->defaultValue = label->defaultValue;
588: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
589: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
590: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
591: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
592: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
593: for (PetscInt v = 0; v < label->numStrata; ++v) {
594: (*labelnew)->stratumValues[v] = label->stratumValues[v];
595: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
596: (*labelnew)->validIS[v] = PETSC_TRUE;
597: }
598: (*labelnew)->pStart = -1;
599: (*labelnew)->pEnd = -1;
600: (*labelnew)->bt = NULL;
601: PetscUseTypeMethod(label, duplicate, labelnew);
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@C
606: DMLabelCompare - Compare two `DMLabel` objects
608: Collective; No Fortran Support
610: Input Parameters:
611: + comm - Comm over which to compare labels
612: . l0 - First `DMLabel`
613: - l1 - Second `DMLabel`
615: Output Parameters
616: + equal - (Optional) Flag whether the two labels are equal
617: - message - (Optional) Message describing the difference
619: Level: intermediate
621: Notes:
622: The output flag equal is the same on all processes.
623: If it is passed as `NULL` and difference is found, an error is thrown on all processes.
624: Make sure to pass `NULL` on all processes.
626: The output message is set independently on each rank.
627: It is set to `NULL` if no difference was found on the current rank. It must be freed by user.
628: If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner.
629: Make sure to pass `NULL` on all processes.
631: For the comparison, we ignore the order of stratum values, and strata with no points.
633: The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel.
635: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
636: @*/
637: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
638: {
639: const char *name0, *name1;
640: char msg[PETSC_MAX_PATH_LEN] = "";
641: PetscBool eq;
642: PetscMPIInt rank;
644: PetscFunctionBegin;
649: PetscCallMPI(MPI_Comm_rank(comm, &rank));
650: PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
651: PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
652: {
653: PetscInt v0, v1;
655: PetscCall(DMLabelGetDefaultValue(l0, &v0));
656: PetscCall(DMLabelGetDefaultValue(l1, &v1));
657: eq = (PetscBool)(v0 == v1);
658: 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));
659: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
660: if (!eq) goto finish;
661: }
662: {
663: IS is0, is1;
665: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
666: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
667: PetscCall(ISEqual(is0, is1, &eq));
668: PetscCall(ISDestroy(&is0));
669: PetscCall(ISDestroy(&is1));
670: if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
671: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
672: if (!eq) goto finish;
673: }
674: {
675: PetscInt i, nValues;
677: PetscCall(DMLabelGetNumValues(l0, &nValues));
678: for (i = 0; i < nValues; i++) {
679: const PetscInt v = l0->stratumValues[i];
680: PetscInt n;
681: IS is0, is1;
683: PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
684: if (!n) continue;
685: PetscCall(DMLabelGetStratumIS(l0, v, &is0));
686: PetscCall(DMLabelGetStratumIS(l1, v, &is1));
687: PetscCall(ISEqualUnsorted(is0, is1, &eq));
688: PetscCall(ISDestroy(&is0));
689: PetscCall(ISDestroy(&is1));
690: if (!eq) {
691: 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));
692: break;
693: }
694: }
695: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
696: }
697: finish:
698: /* If message output arg not set, print to stderr */
699: if (message) {
700: *message = NULL;
701: if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
702: } else {
703: if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
704: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
705: }
706: /* If same output arg not ser and labels are not equal, throw error */
707: if (equal) *equal = eq;
708: else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: /*@
713: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
715: Not Collective
717: Input Parameter:
718: . label - The `DMLabel`
720: Level: intermediate
722: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
723: @*/
724: PetscErrorCode DMLabelComputeIndex(DMLabel label)
725: {
726: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
728: PetscFunctionBegin;
730: PetscCall(DMLabelMakeAllValid_Private(label));
731: for (v = 0; v < label->numStrata; ++v) {
732: const PetscInt *points;
733: PetscInt i;
735: PetscCall(ISGetIndices(label->points[v], &points));
736: for (i = 0; i < label->stratumSizes[v]; ++i) {
737: const PetscInt point = points[i];
739: pStart = PetscMin(point, pStart);
740: pEnd = PetscMax(point + 1, pEnd);
741: }
742: PetscCall(ISRestoreIndices(label->points[v], &points));
743: }
744: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
745: label->pEnd = pEnd;
746: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: /*@
751: DMLabelCreateIndex - Create an index structure for membership determination
753: Not Collective
755: Input Parameters:
756: + label - The `DMLabel`
757: . pStart - The smallest point
758: - pEnd - The largest point + 1
760: Level: intermediate
762: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
763: @*/
764: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
765: {
766: PetscInt v;
768: PetscFunctionBegin;
770: PetscCall(DMLabelDestroyIndex(label));
771: PetscCall(DMLabelMakeAllValid_Private(label));
772: label->pStart = pStart;
773: label->pEnd = pEnd;
774: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
775: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
776: for (v = 0; v < label->numStrata; ++v) {
777: IS pointIS;
778: const PetscInt *points;
779: PetscInt i;
781: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
782: PetscCall(ISGetIndices(pointIS, &points));
783: for (i = 0; i < label->stratumSizes[v]; ++i) {
784: const PetscInt point = points[i];
786: 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);
787: PetscCall(PetscBTSet(label->bt, point - pStart));
788: }
789: PetscCall(ISRestoreIndices(label->points[v], &points));
790: PetscCall(ISDestroy(&pointIS));
791: }
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }
795: /*@
796: DMLabelDestroyIndex - Destroy the index structure
798: Not Collective
800: Input Parameter:
801: . label - the `DMLabel`
803: Level: intermediate
805: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
806: @*/
807: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
808: {
809: PetscFunctionBegin;
811: label->pStart = -1;
812: label->pEnd = -1;
813: PetscCall(PetscBTDestroy(&label->bt));
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /*@
818: DMLabelGetBounds - Return the smallest and largest point in the label
820: Not Collective
822: Input Parameter:
823: . label - the `DMLabel`
825: Output Parameters:
826: + pStart - The smallest point
827: - pEnd - The largest point + 1
829: Level: intermediate
831: Note:
832: This will compute an index for the label if one does not exist.
834: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
835: @*/
836: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
837: {
838: PetscFunctionBegin;
840: if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
841: if (pStart) {
843: *pStart = label->pStart;
844: }
845: if (pEnd) {
847: *pEnd = label->pEnd;
848: }
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: DMLabelHasValue - Determine whether a label assigns the value to any point
855: Not Collective
857: Input Parameters:
858: + label - the `DMLabel`
859: - value - the value
861: Output Parameter:
862: . contains - Flag indicating whether the label maps this value to any point
864: Level: developer
866: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
867: @*/
868: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
869: {
870: PetscInt v;
872: PetscFunctionBegin;
875: PetscCall(DMLabelLookupStratum(label, value, &v));
876: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
877: PetscFunctionReturn(PETSC_SUCCESS);
878: }
880: /*@
881: DMLabelHasPoint - Determine whether a label assigns a value to a point
883: Not Collective
885: Input Parameters:
886: + label - the `DMLabel`
887: - point - the point
889: Output Parameter:
890: . contains - Flag indicating whether the label maps this point to a value
892: Level: developer
894: Note:
895: The user must call `DMLabelCreateIndex()` before this function.
897: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
898: @*/
899: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
900: {
901: PetscFunctionBeginHot;
904: PetscCall(DMLabelMakeAllValid_Private(label));
905: if (PetscDefined(USE_DEBUG)) {
906: PetscCheck(label->bt, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
907: 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);
908: }
909: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@
914: DMLabelStratumHasPoint - Return true if the stratum contains a point
916: Not Collective
918: Input Parameters:
919: + label - the `DMLabel`
920: . value - the stratum value
921: - point - the point
923: Output Parameter:
924: . contains - true if the stratum contains the point
926: Level: intermediate
928: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
929: @*/
930: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
931: {
932: PetscFunctionBeginHot;
935: if (value == label->defaultValue) {
936: PetscInt pointVal;
938: PetscCall(DMLabelGetValue(label, point, &pointVal));
939: *contains = (PetscBool)(pointVal == value);
940: } else {
941: PetscInt v;
943: PetscCall(DMLabelLookupStratum(label, value, &v));
944: if (v >= 0) {
945: if (label->validIS[v] || label->readonly) {
946: IS is;
947: PetscInt i;
949: PetscUseTypeMethod(label, getstratumis, v, &is);
950: PetscCall(ISLocate(is, point, &i));
951: PetscCall(ISDestroy(&is));
952: *contains = (PetscBool)(i >= 0);
953: } else {
954: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
955: }
956: } else { // value is not present
957: *contains = PETSC_FALSE;
958: }
959: }
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: /*@
964: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
965: When a label is created, it is initialized to -1.
967: Not Collective
969: Input parameter:
970: . label - a `DMLabel` object
972: Output parameter:
973: . defaultValue - the default value
975: Level: beginner
977: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
978: @*/
979: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
980: {
981: PetscFunctionBegin;
983: *defaultValue = label->defaultValue;
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: DMLabelSetDefaultValue - Set 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`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1002: @*/
1003: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1004: {
1005: PetscFunctionBegin;
1007: label->defaultValue = defaultValue;
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@
1012: 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
1013: `DMLabelSetDefaultValue()`)
1015: Not Collective
1017: Input Parameters:
1018: + label - the `DMLabel`
1019: - point - the point
1021: Output Parameter:
1022: . value - The point value, or the default value (-1 by default)
1024: Level: intermediate
1026: Note:
1027: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1028: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1030: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1031: @*/
1032: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1033: {
1034: PetscInt v;
1036: PetscFunctionBeginHot;
1039: *value = label->defaultValue;
1040: for (v = 0; v < label->numStrata; ++v) {
1041: if (label->validIS[v] || label->readonly) {
1042: IS is;
1043: PetscInt i;
1045: PetscUseTypeMethod(label, getstratumis, v, &is);
1046: PetscCall(ISLocate(label->points[v], point, &i));
1047: PetscCall(ISDestroy(&is));
1048: if (i >= 0) {
1049: *value = label->stratumValues[v];
1050: break;
1051: }
1052: } else {
1053: PetscBool has;
1055: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1056: if (has) {
1057: *value = label->stratumValues[v];
1058: break;
1059: }
1060: }
1061: }
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: /*@
1066: 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
1067: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1069: Not Collective
1071: Input Parameters:
1072: + label - the `DMLabel`
1073: . point - the point
1074: - value - The point value
1076: Level: intermediate
1078: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1079: @*/
1080: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1081: {
1082: PetscInt v;
1084: PetscFunctionBegin;
1086: /* Find label value, add new entry if needed */
1087: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1088: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1089: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1090: /* Set key */
1091: PetscCall(DMLabelMakeInvalid_Private(label, v));
1092: PetscCall(PetscHSetIAdd(label->ht[v], point));
1093: PetscFunctionReturn(PETSC_SUCCESS);
1094: }
1096: /*@
1097: DMLabelClearValue - Clear the value a label assigns to a point
1099: Not Collective
1101: Input Parameters:
1102: + label - the `DMLabel`
1103: . point - the point
1104: - value - The point value
1106: Level: intermediate
1108: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1109: @*/
1110: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1111: {
1112: PetscInt v;
1114: PetscFunctionBegin;
1116: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1117: /* Find label value */
1118: PetscCall(DMLabelLookupStratum(label, value, &v));
1119: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1121: if (label->bt) {
1122: 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);
1123: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1124: }
1126: /* Delete key */
1127: PetscCall(DMLabelMakeInvalid_Private(label, v));
1128: PetscCall(PetscHSetIDel(label->ht[v], point));
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }
1132: /*@
1133: DMLabelInsertIS - Set all points in the `IS` to a value
1135: Not Collective
1137: Input Parameters:
1138: + label - the `DMLabel`
1139: . is - the point `IS`
1140: - value - The point value
1142: Level: intermediate
1144: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1145: @*/
1146: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1147: {
1148: PetscInt v, n, p;
1149: const PetscInt *points;
1151: PetscFunctionBegin;
1154: /* Find label value, add new entry if needed */
1155: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1156: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1157: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1158: /* Set keys */
1159: PetscCall(DMLabelMakeInvalid_Private(label, v));
1160: PetscCall(ISGetLocalSize(is, &n));
1161: PetscCall(ISGetIndices(is, &points));
1162: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1163: PetscCall(ISRestoreIndices(is, &points));
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: /*@
1168: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1170: Not Collective
1172: Input Parameter:
1173: . label - the `DMLabel`
1175: Output Parameter:
1176: . numValues - the number of values
1178: Level: intermediate
1180: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1181: @*/
1182: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1183: {
1184: PetscFunctionBegin;
1187: *numValues = label->numStrata;
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: /*@
1192: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1194: Not Collective
1196: Input Parameter:
1197: . label - the `DMLabel`
1199: Output Parameter:
1200: . is - the value `IS`
1202: Level: intermediate
1204: Notes:
1205: The output `IS` should be destroyed when no longer needed.
1206: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1207: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1209: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1210: @*/
1211: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1212: {
1213: PetscFunctionBegin;
1216: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: /*@
1221: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1223: Not Collective
1225: Input Parameter:
1226: . label - the `DMLabel`
1228: Output Parameter:
1229: . is - the value `IS`
1231: Level: intermediate
1233: Notes:
1234: The output `IS` should be destroyed when no longer needed.
1235: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1237: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1238: @*/
1239: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1240: {
1241: PetscInt i, j;
1242: PetscInt *valuesArr;
1244: PetscFunctionBegin;
1247: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1248: for (i = 0, j = 0; i < label->numStrata; i++) {
1249: PetscInt n;
1251: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1252: if (n) valuesArr[j++] = label->stratumValues[i];
1253: }
1254: if (j == label->numStrata) {
1255: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1256: } else {
1257: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1258: }
1259: PetscCall(PetscFree(valuesArr));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*@
1264: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1266: Not Collective
1268: Input Parameters:
1269: + label - the `DMLabel`
1270: - value - the value
1272: Output Parameter:
1273: . index - the index of value in the list of values
1275: Level: intermediate
1277: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1278: @*/
1279: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1280: {
1281: PetscInt v;
1283: PetscFunctionBegin;
1286: /* Do not assume they are sorted */
1287: for (v = 0; v < label->numStrata; ++v)
1288: if (label->stratumValues[v] == value) break;
1289: if (v >= label->numStrata) *index = -1;
1290: else *index = v;
1291: PetscFunctionReturn(PETSC_SUCCESS);
1292: }
1294: /*@
1295: DMLabelHasStratum - Determine whether points exist with the given value
1297: Not Collective
1299: Input Parameters:
1300: + label - the `DMLabel`
1301: - value - the stratum value
1303: Output Parameter:
1304: . exists - Flag saying whether points exist
1306: Level: intermediate
1308: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1309: @*/
1310: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1311: {
1312: PetscInt v;
1314: PetscFunctionBegin;
1317: PetscCall(DMLabelLookupStratum(label, value, &v));
1318: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1319: PetscFunctionReturn(PETSC_SUCCESS);
1320: }
1322: /*@
1323: DMLabelGetStratumSize - Get the size of a stratum
1325: Not Collective
1327: Input Parameters:
1328: + label - the `DMLabel`
1329: - value - the stratum value
1331: Output Parameter:
1332: . size - The number of points in the stratum
1334: Level: intermediate
1336: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1337: @*/
1338: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1339: {
1340: PetscInt v;
1342: PetscFunctionBegin;
1345: PetscCall(DMLabelLookupStratum(label, value, &v));
1346: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }
1350: /*@
1351: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1353: Not Collective
1355: Input Parameters:
1356: + label - the `DMLabel`
1357: - value - the stratum value
1359: Output Parameters:
1360: + start - the smallest point in the stratum
1361: - end - the largest point in the stratum
1363: Level: intermediate
1365: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1366: @*/
1367: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1368: {
1369: IS is;
1370: PetscInt v, min, max;
1372: PetscFunctionBegin;
1374: if (start) {
1376: *start = -1;
1377: }
1378: if (end) {
1380: *end = -1;
1381: }
1382: PetscCall(DMLabelLookupStratum(label, value, &v));
1383: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1384: PetscCall(DMLabelMakeValid_Private(label, v));
1385: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1386: PetscUseTypeMethod(label, getstratumis, v, &is);
1387: PetscCall(ISGetMinMax(is, &min, &max));
1388: PetscCall(ISDestroy(&is));
1389: if (start) *start = min;
1390: if (end) *end = max + 1;
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1395: {
1396: PetscFunctionBegin;
1397: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1398: *pointIS = label->points[v];
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: /*@
1403: DMLabelGetStratumIS - Get an `IS` with the stratum points
1405: Not Collective
1407: Input Parameters:
1408: + label - the `DMLabel`
1409: - value - the stratum value
1411: Output Parameter:
1412: . points - The stratum points
1414: Level: intermediate
1416: Notes:
1417: The output `IS` should be destroyed when no longer needed.
1418: Returns `NULL` if the stratum is empty.
1420: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1421: @*/
1422: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1423: {
1424: PetscInt v;
1426: PetscFunctionBegin;
1429: *points = NULL;
1430: PetscCall(DMLabelLookupStratum(label, value, &v));
1431: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1432: PetscCall(DMLabelMakeValid_Private(label, v));
1433: PetscUseTypeMethod(label, getstratumis, v, points);
1434: PetscFunctionReturn(PETSC_SUCCESS);
1435: }
1437: /*@
1438: DMLabelSetStratumIS - Set the stratum points using an `IS`
1440: Not Collective
1442: Input Parameters:
1443: + label - the `DMLabel`
1444: . value - the stratum value
1445: - points - The stratum points
1447: Level: intermediate
1449: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1450: @*/
1451: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1452: {
1453: PetscInt v;
1455: PetscFunctionBegin;
1458: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1459: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1460: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1461: PetscCall(DMLabelClearStratum(label, value));
1462: PetscCall(ISGetLocalSize(is, &(label->stratumSizes[v])));
1463: PetscCall(PetscObjectReference((PetscObject)is));
1464: PetscCall(ISDestroy(&(label->points[v])));
1465: label->points[v] = is;
1466: label->validIS[v] = PETSC_TRUE;
1467: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1468: if (label->bt) {
1469: const PetscInt *points;
1470: PetscInt p;
1472: PetscCall(ISGetIndices(is, &points));
1473: for (p = 0; p < label->stratumSizes[v]; ++p) {
1474: const PetscInt point = points[p];
1476: 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);
1477: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1478: }
1479: }
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: /*@
1484: DMLabelClearStratum - Remove a stratum
1486: Not Collective
1488: Input Parameters:
1489: + label - the `DMLabel`
1490: - value - the stratum value
1492: Level: intermediate
1494: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1495: @*/
1496: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1497: {
1498: PetscInt v;
1500: PetscFunctionBegin;
1502: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1503: PetscCall(DMLabelLookupStratum(label, value, &v));
1504: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1505: if (label->validIS[v]) {
1506: if (label->bt) {
1507: PetscInt i;
1508: const PetscInt *points;
1510: PetscCall(ISGetIndices(label->points[v], &points));
1511: for (i = 0; i < label->stratumSizes[v]; ++i) {
1512: const PetscInt point = points[i];
1514: 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);
1515: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1516: }
1517: PetscCall(ISRestoreIndices(label->points[v], &points));
1518: }
1519: label->stratumSizes[v] = 0;
1520: PetscCall(ISDestroy(&label->points[v]));
1521: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1522: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1523: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1524: } else {
1525: PetscCall(PetscHSetIClear(label->ht[v]));
1526: }
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: /*@
1531: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1533: Not Collective
1535: Input Parameters:
1536: + label - The `DMLabel`
1537: . value - The label value for all points
1538: . pStart - The first point
1539: - pEnd - A point beyond all marked points
1541: Level: intermediate
1543: Note:
1544: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1546: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1547: @*/
1548: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1549: {
1550: IS pIS;
1552: PetscFunctionBegin;
1553: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1554: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1555: PetscCall(ISDestroy(&pIS));
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: /*@
1560: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1562: Not Collective
1564: Input Parameters:
1565: + label - The `DMLabel`
1566: . value - The label value
1567: - p - A point with this value
1569: Output Parameter:
1570: . 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
1572: Level: intermediate
1574: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1575: @*/
1576: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1577: {
1578: IS pointIS;
1579: const PetscInt *indices;
1580: PetscInt v;
1582: PetscFunctionBegin;
1585: *index = -1;
1586: PetscCall(DMLabelLookupStratum(label, value, &v));
1587: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1588: PetscCall(DMLabelMakeValid_Private(label, v));
1589: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1590: PetscCall(ISGetIndices(pointIS, &indices));
1591: PetscCall(PetscFindInt(p, label->stratumSizes[v], indices, index));
1592: PetscCall(ISRestoreIndices(pointIS, &indices));
1593: PetscCall(ISDestroy(&pointIS));
1594: PetscFunctionReturn(PETSC_SUCCESS);
1595: }
1597: /*@
1598: DMLabelFilter - Remove all points outside of [`start`, `end`)
1600: Not Collective
1602: Input Parameters:
1603: + label - the `DMLabel`
1604: . start - the first point kept
1605: - end - one more than the last point kept
1607: Level: intermediate
1609: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1610: @*/
1611: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1612: {
1613: PetscInt v;
1615: PetscFunctionBegin;
1617: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1618: PetscCall(DMLabelDestroyIndex(label));
1619: PetscCall(DMLabelMakeAllValid_Private(label));
1620: for (v = 0; v < label->numStrata; ++v) {
1621: PetscCall(ISGeneralFilter(label->points[v], start, end));
1622: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1623: }
1624: PetscCall(DMLabelCreateIndex(label, start, end));
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: /*@
1629: DMLabelPermute - Create a new label with permuted points
1631: Not Collective
1633: Input Parameters:
1634: + label - the `DMLabel`
1635: - permutation - the point permutation
1637: Output Parameter:
1638: . labelnew - the new label containing the permuted points
1640: Level: intermediate
1642: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1643: @*/
1644: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1645: {
1646: const PetscInt *perm;
1647: PetscInt numValues, numPoints, v, q;
1649: PetscFunctionBegin;
1652: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1653: PetscCall(DMLabelMakeAllValid_Private(label));
1654: PetscCall(DMLabelDuplicate(label, labelNew));
1655: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1656: PetscCall(ISGetLocalSize(permutation, &numPoints));
1657: PetscCall(ISGetIndices(permutation, &perm));
1658: for (v = 0; v < numValues; ++v) {
1659: const PetscInt size = (*labelNew)->stratumSizes[v];
1660: const PetscInt *points;
1661: PetscInt *pointsNew;
1663: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1664: PetscCall(PetscCalloc1(size, &pointsNew));
1665: for (q = 0; q < size; ++q) {
1666: const PetscInt point = points[q];
1668: 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);
1669: pointsNew[q] = perm[point];
1670: }
1671: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1672: PetscCall(PetscSortInt(size, pointsNew));
1673: PetscCall(ISDestroy(&((*labelNew)->points[v])));
1674: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1675: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1676: PetscCall(PetscFree(pointsNew));
1677: } else {
1678: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1679: }
1680: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1681: }
1682: PetscCall(ISRestoreIndices(permutation, &perm));
1683: if (label->bt) {
1684: PetscCall(PetscBTDestroy(&label->bt));
1685: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1686: }
1687: PetscFunctionReturn(PETSC_SUCCESS);
1688: }
1690: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1691: {
1692: MPI_Comm comm;
1693: PetscInt s, l, nroots, nleaves, offset, size;
1694: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1695: PetscSection rootSection;
1696: PetscSF labelSF;
1698: PetscFunctionBegin;
1699: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1700: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1701: /* Build a section of stratum values per point, generate the according SF
1702: and distribute point-wise stratum values to leaves. */
1703: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1704: PetscCall(PetscSectionCreate(comm, &rootSection));
1705: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1706: if (label) {
1707: for (s = 0; s < label->numStrata; ++s) {
1708: const PetscInt *points;
1710: PetscCall(ISGetIndices(label->points[s], &points));
1711: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1712: PetscCall(ISRestoreIndices(label->points[s], &points));
1713: }
1714: }
1715: PetscCall(PetscSectionSetUp(rootSection));
1716: /* Create a point-wise array of stratum values */
1717: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1718: PetscCall(PetscMalloc1(size, &rootStrata));
1719: PetscCall(PetscCalloc1(nroots, &rootIdx));
1720: if (label) {
1721: for (s = 0; s < label->numStrata; ++s) {
1722: const PetscInt *points;
1724: PetscCall(ISGetIndices(label->points[s], &points));
1725: for (l = 0; l < label->stratumSizes[s]; l++) {
1726: const PetscInt p = points[l];
1727: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1728: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1729: }
1730: PetscCall(ISRestoreIndices(label->points[s], &points));
1731: }
1732: }
1733: /* Build SF that maps label points to remote processes */
1734: PetscCall(PetscSectionCreate(comm, leafSection));
1735: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1736: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1737: PetscCall(PetscFree(remoteOffsets));
1738: /* Send the strata for each point over the derived SF */
1739: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1740: PetscCall(PetscMalloc1(size, leafStrata));
1741: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1742: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1743: /* Clean up */
1744: PetscCall(PetscFree(rootStrata));
1745: PetscCall(PetscFree(rootIdx));
1746: PetscCall(PetscSectionDestroy(&rootSection));
1747: PetscCall(PetscSFDestroy(&labelSF));
1748: PetscFunctionReturn(PETSC_SUCCESS);
1749: }
1751: /*@
1752: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1754: Collective
1756: Input Parameters:
1757: + label - the `DMLabel`
1758: - sf - the map from old to new distribution
1760: Output Parameter:
1761: . labelnew - the new redistributed label
1763: Level: intermediate
1765: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1766: @*/
1767: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1768: {
1769: MPI_Comm comm;
1770: PetscSection leafSection;
1771: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1772: PetscInt *leafStrata, *strataIdx;
1773: PetscInt **points;
1774: const char *lname = NULL;
1775: char *name;
1776: PetscInt nameSize;
1777: PetscHSetI stratumHash;
1778: size_t len = 0;
1779: PetscMPIInt rank;
1781: PetscFunctionBegin;
1783: if (label) {
1785: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1786: PetscCall(DMLabelMakeAllValid_Private(label));
1787: }
1788: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1789: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1790: /* Bcast name */
1791: if (rank == 0) {
1792: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1793: PetscCall(PetscStrlen(lname, &len));
1794: }
1795: nameSize = len;
1796: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm));
1797: PetscCall(PetscMalloc1(nameSize + 1, &name));
1798: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1799: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1800: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1801: PetscCall(PetscFree(name));
1802: /* Bcast defaultValue */
1803: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1804: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
1805: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1806: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
1807: /* Determine received stratum values and initialise new label*/
1808: PetscCall(PetscHSetICreate(&stratumHash));
1809: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
1810: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
1811: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
1812: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
1813: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1814: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
1815: /* Turn leafStrata into indices rather than stratum values */
1816: offset = 0;
1817: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
1818: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
1819: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
1820: for (p = 0; p < size; ++p) {
1821: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1822: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1823: leafStrata[p] = s;
1824: break;
1825: }
1826: }
1827: }
1828: /* Rebuild the point strata on the receiver */
1829: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
1830: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1831: for (p = pStart; p < pEnd; p++) {
1832: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1833: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1834: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1835: }
1836: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
1837: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
1838: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
1839: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1840: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
1841: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s])));
1842: }
1843: /* Insert points into new strata */
1844: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
1845: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1846: for (p = pStart; p < pEnd; p++) {
1847: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1848: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1849: for (s = 0; s < dof; s++) {
1850: stratum = leafStrata[offset + s];
1851: points[stratum][strataIdx[stratum]++] = p;
1852: }
1853: }
1854: for (s = 0; s < (*labelNew)->numStrata; s++) {
1855: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &(points[s][0]), PETSC_OWN_POINTER, &((*labelNew)->points[s])));
1856: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
1857: }
1858: PetscCall(PetscFree(points));
1859: PetscCall(PetscHSetIDestroy(&stratumHash));
1860: PetscCall(PetscFree(leafStrata));
1861: PetscCall(PetscFree(strataIdx));
1862: PetscCall(PetscSectionDestroy(&leafSection));
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1866: /*@
1867: DMLabelGather - Gather all label values from leafs into roots
1869: Collective
1871: Input Parameters:
1872: + label - the `DMLabel`
1873: - sf - the `PetscSF` communication map
1875: Output Parameter:
1876: . labelNew - the new `DMLabel` with localised leaf values
1878: Level: developer
1880: Note:
1881: This is the inverse operation to `DMLabelDistribute()`.
1883: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
1884: @*/
1885: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1886: {
1887: MPI_Comm comm;
1888: PetscSection rootSection;
1889: PetscSF sfLabel;
1890: PetscSFNode *rootPoints, *leafPoints;
1891: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1892: const PetscInt *rootDegree, *ilocal;
1893: PetscInt *rootStrata;
1894: const char *lname;
1895: char *name;
1896: PetscInt nameSize;
1897: size_t len = 0;
1898: PetscMPIInt rank, size;
1900: PetscFunctionBegin;
1903: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1904: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1905: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1906: PetscCallMPI(MPI_Comm_size(comm, &size));
1907: /* Bcast name */
1908: if (rank == 0) {
1909: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1910: PetscCall(PetscStrlen(lname, &len));
1911: }
1912: nameSize = len;
1913: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm));
1914: PetscCall(PetscMalloc1(nameSize + 1, &name));
1915: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1916: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1917: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1918: PetscCall(PetscFree(name));
1919: /* Gather rank/index pairs of leaves into local roots to build
1920: an inverse, multi-rooted SF. Note that this ignores local leaf
1921: indexing due to the use of the multiSF in PetscSFGather. */
1922: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
1923: PetscCall(PetscMalloc1(nroots, &leafPoints));
1924: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1925: for (p = 0; p < nleaves; p++) {
1926: PetscInt ilp = ilocal ? ilocal[p] : p;
1928: leafPoints[ilp].index = ilp;
1929: leafPoints[ilp].rank = rank;
1930: }
1931: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
1932: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
1933: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1934: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
1935: PetscCall(PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints));
1936: PetscCall(PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints));
1937: PetscCall(PetscSFCreate(comm, &sfLabel));
1938: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
1939: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1940: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
1941: /* Rebuild the point strata on the receiver */
1942: for (p = 0, idx = 0; p < nroots; p++) {
1943: for (d = 0; d < rootDegree[p]; d++) {
1944: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
1945: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
1946: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
1947: }
1948: idx += rootDegree[p];
1949: }
1950: PetscCall(PetscFree(leafPoints));
1951: PetscCall(PetscFree(rootStrata));
1952: PetscCall(PetscSectionDestroy(&rootSection));
1953: PetscCall(PetscSFDestroy(&sfLabel));
1954: PetscFunctionReturn(PETSC_SUCCESS);
1955: }
1957: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
1958: {
1959: const PetscInt *degree;
1960: const PetscInt *points;
1961: PetscInt Nr, r, Nl, l, val, defVal;
1963: PetscFunctionBegin;
1964: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1965: /* Add in leaves */
1966: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
1967: for (l = 0; l < Nl; ++l) {
1968: PetscCall(DMLabelGetValue(label, points[l], &val));
1969: if (val != defVal) valArray[points[l]] = val;
1970: }
1971: /* Add in shared roots */
1972: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
1973: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
1974: for (r = 0; r < Nr; ++r) {
1975: if (degree[r]) {
1976: PetscCall(DMLabelGetValue(label, r, &val));
1977: if (val != defVal) valArray[r] = val;
1978: }
1979: }
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1984: {
1985: const PetscInt *degree;
1986: const PetscInt *points;
1987: PetscInt Nr, r, Nl, l, val, defVal;
1989: PetscFunctionBegin;
1990: PetscCall(DMLabelGetDefaultValue(label, &defVal));
1991: /* Read out leaves */
1992: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
1993: for (l = 0; l < Nl; ++l) {
1994: const PetscInt p = points[l];
1995: const PetscInt cval = valArray[p];
1997: if (cval != defVal) {
1998: PetscCall(DMLabelGetValue(label, p, &val));
1999: if (val == defVal) {
2000: PetscCall(DMLabelSetValue(label, p, cval));
2001: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2002: }
2003: }
2004: }
2005: /* Read out shared roots */
2006: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2007: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2008: for (r = 0; r < Nr; ++r) {
2009: if (degree[r]) {
2010: const PetscInt cval = valArray[r];
2012: if (cval != defVal) {
2013: PetscCall(DMLabelGetValue(label, r, &val));
2014: if (val == defVal) {
2015: PetscCall(DMLabelSetValue(label, r, cval));
2016: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2017: }
2018: }
2019: }
2020: }
2021: PetscFunctionReturn(PETSC_SUCCESS);
2022: }
2024: /*@
2025: DMLabelPropagateBegin - Setup a cycle of label propagation
2027: Collective
2029: Input Parameters:
2030: + label - The `DMLabel` to propagate across processes
2031: - sf - The `PetscSF` describing parallel layout of the label points
2033: Level: intermediate
2035: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2036: @*/
2037: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2038: {
2039: PetscInt Nr, r, defVal;
2040: PetscMPIInt size;
2042: PetscFunctionBegin;
2043: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2044: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2045: if (size > 1) {
2046: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2047: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2048: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2049: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2050: }
2051: PetscFunctionReturn(PETSC_SUCCESS);
2052: }
2054: /*@
2055: DMLabelPropagateEnd - Tear down a cycle of label propagation
2057: Collective
2059: Input Parameters:
2060: + label - The `DMLabel` to propagate across processes
2061: - sf - The `PetscSF` describing parallel layout of the label points
2063: Level: intermediate
2065: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2066: @*/
2067: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2068: {
2069: PetscFunctionBegin;
2070: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2071: PetscCall(PetscFree(label->propArray));
2072: label->propArray = NULL;
2073: PetscFunctionReturn(PETSC_SUCCESS);
2074: }
2076: /*@C
2077: DMLabelPropagatePush - Tear down a cycle of label propagation
2079: Collective
2081: Input Parameters:
2082: + label - The `DMLabel` to propagate across processes
2083: . sf - The `PetscSF` describing parallel layout of the label points
2084: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2085: - ctx - An optional user context for the callback, or `NULL`
2087: Calling sequence of `markPoint`:
2088: $ PetscErrorCode markPoint(DMLabel label, PetscInt p, PetscInt val, void *ctx);
2089: + label - The `DMLabel`
2090: . p - The point being marked
2091: . val - The label value for `p`
2092: - ctx - An optional user context
2094: Level: intermediate
2096: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2097: @*/
2098: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
2099: {
2100: PetscInt *valArray = label->propArray, Nr;
2101: PetscMPIInt size;
2103: PetscFunctionBegin;
2104: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2105: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2106: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2107: if (size > 1 && Nr >= 0) {
2108: /* Communicate marked edges
2109: The current implementation allocates an array the size of the number of root. We put the label values into the
2110: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2112: TODO: We could use in-place communication with a different SF
2113: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2114: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2116: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2117: 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
2118: edge to the queue.
2119: */
2120: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2121: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2122: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2123: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2124: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2125: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2126: }
2127: PetscFunctionReturn(PETSC_SUCCESS);
2128: }
2130: /*@
2131: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2133: Not Collective
2135: Input Parameter:
2136: . label - the `DMLabel`
2138: Output Parameters:
2139: + section - the section giving offsets for each stratum
2140: - is - An `IS` containing all the label points
2142: Level: developer
2144: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2145: @*/
2146: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2147: {
2148: IS vIS;
2149: const PetscInt *values;
2150: PetscInt *points;
2151: PetscInt nV, vS = 0, vE = 0, v, N;
2153: PetscFunctionBegin;
2155: PetscCall(DMLabelGetNumValues(label, &nV));
2156: PetscCall(DMLabelGetValueIS(label, &vIS));
2157: PetscCall(ISGetIndices(vIS, &values));
2158: if (nV) {
2159: vS = values[0];
2160: vE = values[0] + 1;
2161: }
2162: for (v = 1; v < nV; ++v) {
2163: vS = PetscMin(vS, values[v]);
2164: vE = PetscMax(vE, values[v] + 1);
2165: }
2166: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2167: PetscCall(PetscSectionSetChart(*section, vS, vE));
2168: for (v = 0; v < nV; ++v) {
2169: PetscInt n;
2171: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2172: PetscCall(PetscSectionSetDof(*section, values[v], n));
2173: }
2174: PetscCall(PetscSectionSetUp(*section));
2175: PetscCall(PetscSectionGetStorageSize(*section, &N));
2176: PetscCall(PetscMalloc1(N, &points));
2177: for (v = 0; v < nV; ++v) {
2178: IS is;
2179: const PetscInt *spoints;
2180: PetscInt dof, off, p;
2182: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2183: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2184: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2185: PetscCall(ISGetIndices(is, &spoints));
2186: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2187: PetscCall(ISRestoreIndices(is, &spoints));
2188: PetscCall(ISDestroy(&is));
2189: }
2190: PetscCall(ISRestoreIndices(vIS, &values));
2191: PetscCall(ISDestroy(&vIS));
2192: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2193: PetscFunctionReturn(PETSC_SUCCESS);
2194: }
2196: /*@C
2197: DMLabelRegister - Adds a new label component implementation
2199: Not Collective
2201: Input Parameters:
2202: + name - The name of a new user-defined creation routine
2203: - create_func - The creation routine itself
2205: Notes:
2206: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2208: Sample usage:
2209: .vb
2210: DMLabelRegister("my_label", MyLabelCreate);
2211: .ve
2213: Then, your label type can be chosen with the procedural interface via
2214: .vb
2215: DMLabelCreate(MPI_Comm, DMLabel *);
2216: DMLabelSetType(DMLabel, "my_label");
2217: .ve
2218: or at runtime via the option
2219: .vb
2220: -dm_label_type my_label
2221: .ve
2223: Level: advanced
2225: .seealso: `DMLabel`, `DM`, `DMLabel`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2226: @*/
2227: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2228: {
2229: PetscFunctionBegin;
2230: PetscCall(DMInitializePackage());
2231: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2232: PetscFunctionReturn(PETSC_SUCCESS);
2233: }
2235: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2236: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2238: /*@C
2239: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2241: Not Collective
2243: Level: advanced
2245: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2246: @*/
2247: PetscErrorCode DMLabelRegisterAll(void)
2248: {
2249: PetscFunctionBegin;
2250: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2251: DMLabelRegisterAllCalled = PETSC_TRUE;
2253: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2254: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2255: PetscFunctionReturn(PETSC_SUCCESS);
2256: }
2258: /*@C
2259: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2261: Level: developer
2263: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2264: @*/
2265: PetscErrorCode DMLabelRegisterDestroy(void)
2266: {
2267: PetscFunctionBegin;
2268: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2269: DMLabelRegisterAllCalled = PETSC_FALSE;
2270: PetscFunctionReturn(PETSC_SUCCESS);
2271: }
2273: /*@C
2274: DMLabelSetType - Sets the particular implementation for a label.
2276: Collective
2278: Input Parameters:
2279: + label - The label
2280: - method - The name of the label type
2282: Options Database Key:
2283: . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType`
2285: Level: intermediate
2287: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2288: @*/
2289: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2290: {
2291: PetscErrorCode (*r)(DMLabel);
2292: PetscBool match;
2294: PetscFunctionBegin;
2296: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2297: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2299: PetscCall(DMLabelRegisterAll());
2300: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2301: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2303: PetscTryTypeMethod(label, destroy);
2304: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2305: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2306: PetscCall((*r)(label));
2307: PetscFunctionReturn(PETSC_SUCCESS);
2308: }
2310: /*@C
2311: DMLabelGetType - Gets the type name (as a string) from the label.
2313: Not Collective
2315: Input Parameter:
2316: . label - The `DMLabel`
2318: Output Parameter:
2319: . type - The `DMLabel` type name
2321: Level: intermediate
2323: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2324: @*/
2325: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2326: {
2327: PetscFunctionBegin;
2330: PetscCall(DMLabelRegisterAll());
2331: *type = ((PetscObject)label)->type_name;
2332: PetscFunctionReturn(PETSC_SUCCESS);
2333: }
2335: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2336: {
2337: PetscFunctionBegin;
2338: label->ops->view = DMLabelView_Concrete;
2339: label->ops->setup = NULL;
2340: label->ops->duplicate = DMLabelDuplicate_Concrete;
2341: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2342: PetscFunctionReturn(PETSC_SUCCESS);
2343: }
2345: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2346: {
2347: PetscFunctionBegin;
2349: PetscCall(DMLabelInitialize_Concrete(label));
2350: PetscFunctionReturn(PETSC_SUCCESS);
2351: }
2353: /*@
2354: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2355: the local section and an `PetscSF` describing the section point overlap.
2357: Collective
2359: Input Parameters:
2360: + s - The `PetscSection` for the local field layout
2361: . sf - The `PetscSF` describing parallel layout of the section points
2362: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2363: . label - The label specifying the points
2364: - labelValue - The label stratum specifying the points
2366: Output Parameter:
2367: . gsection - The `PetscSection` for the global field layout
2369: Level: developer
2371: Note:
2372: This gives negative sizes and offsets to points not owned by this process
2374: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2375: @*/
2376: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2377: {
2378: PetscInt *neg = NULL, *tmpOff = NULL;
2379: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2381: PetscFunctionBegin;
2385: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2386: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2387: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2388: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2389: if (nroots >= 0) {
2390: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2391: PetscCall(PetscCalloc1(nroots, &neg));
2392: if (nroots > pEnd - pStart) {
2393: PetscCall(PetscCalloc1(nroots, &tmpOff));
2394: } else {
2395: tmpOff = &(*gsection)->atlasDof[-pStart];
2396: }
2397: }
2398: /* Mark ghost points with negative dof */
2399: for (p = pStart; p < pEnd; ++p) {
2400: PetscInt value;
2402: PetscCall(DMLabelGetValue(label, p, &value));
2403: if (value != labelValue) continue;
2404: PetscCall(PetscSectionGetDof(s, p, &dof));
2405: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2406: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2407: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2408: if (neg) neg[p] = -(dof + 1);
2409: }
2410: PetscCall(PetscSectionSetUpBC(*gsection));
2411: if (nroots >= 0) {
2412: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2413: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2414: if (nroots > pEnd - pStart) {
2415: for (p = pStart; p < pEnd; ++p) {
2416: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2417: }
2418: }
2419: }
2420: /* Calculate new sizes, get process offset, and calculate point offsets */
2421: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2422: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2423: (*gsection)->atlasOff[p] = off;
2424: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2425: }
2426: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2427: globalOff -= off;
2428: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2429: (*gsection)->atlasOff[p] += globalOff;
2430: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2431: }
2432: /* Put in negative offsets for ghost points */
2433: if (nroots >= 0) {
2434: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2435: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2436: if (nroots > pEnd - pStart) {
2437: for (p = pStart; p < pEnd; ++p) {
2438: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2439: }
2440: }
2441: }
2442: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2443: PetscCall(PetscFree(neg));
2444: PetscFunctionReturn(PETSC_SUCCESS);
2445: }
2447: typedef struct _n_PetscSectionSym_Label {
2448: DMLabel label;
2449: PetscCopyMode *modes;
2450: PetscInt *sizes;
2451: const PetscInt ***perms;
2452: const PetscScalar ***rots;
2453: PetscInt (*minMaxOrients)[2];
2454: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2455: } PetscSectionSym_Label;
2457: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2458: {
2459: PetscInt i, j;
2460: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2462: PetscFunctionBegin;
2463: for (i = 0; i <= sl->numStrata; i++) {
2464: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2465: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2466: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2467: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2468: }
2469: if (sl->perms[i]) {
2470: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2472: PetscCall(PetscFree(perms));
2473: }
2474: if (sl->rots[i]) {
2475: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2477: PetscCall(PetscFree(rots));
2478: }
2479: }
2480: }
2481: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2482: PetscCall(DMLabelDestroy(&sl->label));
2483: sl->numStrata = 0;
2484: PetscFunctionReturn(PETSC_SUCCESS);
2485: }
2487: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2488: {
2489: PetscFunctionBegin;
2490: PetscCall(PetscSectionSymLabelReset(sym));
2491: PetscCall(PetscFree(sym->data));
2492: PetscFunctionReturn(PETSC_SUCCESS);
2493: }
2495: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2496: {
2497: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2498: PetscBool isAscii;
2499: DMLabel label = sl->label;
2500: const char *name;
2502: PetscFunctionBegin;
2503: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2504: if (isAscii) {
2505: PetscInt i, j, k;
2506: PetscViewerFormat format;
2508: PetscCall(PetscViewerGetFormat(viewer, &format));
2509: if (label) {
2510: PetscCall(PetscViewerGetFormat(viewer, &format));
2511: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2512: PetscCall(PetscViewerASCIIPushTab(viewer));
2513: PetscCall(DMLabelView(label, viewer));
2514: PetscCall(PetscViewerASCIIPopTab(viewer));
2515: } else {
2516: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2517: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2518: }
2519: } else {
2520: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2521: }
2522: PetscCall(PetscViewerASCIIPushTab(viewer));
2523: for (i = 0; i <= sl->numStrata; i++) {
2524: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2526: if (!(sl->perms[i] || sl->rots[i])) {
2527: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2528: } else {
2529: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2530: PetscCall(PetscViewerASCIIPushTab(viewer));
2531: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2532: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2533: PetscCall(PetscViewerASCIIPushTab(viewer));
2534: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2535: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2536: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2537: } else {
2538: PetscInt tab;
2540: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2541: PetscCall(PetscViewerASCIIPushTab(viewer));
2542: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2543: if (sl->perms[i] && sl->perms[i][j]) {
2544: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2545: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2546: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2547: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2548: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2549: }
2550: if (sl->rots[i] && sl->rots[i][j]) {
2551: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2552: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2553: #if defined(PETSC_USE_COMPLEX)
2554: 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])));
2555: #else
2556: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2557: #endif
2558: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2559: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2560: }
2561: PetscCall(PetscViewerASCIIPopTab(viewer));
2562: }
2563: }
2564: PetscCall(PetscViewerASCIIPopTab(viewer));
2565: }
2566: PetscCall(PetscViewerASCIIPopTab(viewer));
2567: }
2568: }
2569: PetscCall(PetscViewerASCIIPopTab(viewer));
2570: }
2571: PetscFunctionReturn(PETSC_SUCCESS);
2572: }
2574: /*@
2575: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2577: Logically
2579: Input parameters:
2580: + sym - the section symmetries
2581: - label - the `DMLabel` describing the types of points
2583: Level: developer:
2585: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2586: @*/
2587: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2588: {
2589: PetscSectionSym_Label *sl;
2591: PetscFunctionBegin;
2593: sl = (PetscSectionSym_Label *)sym->data;
2594: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2595: if (label) {
2596: sl->label = label;
2597: PetscCall(PetscObjectReference((PetscObject)label));
2598: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2599: 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));
2600: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2601: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2602: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2603: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2604: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2605: }
2606: PetscFunctionReturn(PETSC_SUCCESS);
2607: }
2609: /*@C
2610: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2612: Logically Collective
2614: Input Parameters:
2615: + sym - the section symmetries
2616: - stratum - the stratum value in the label that we are assigning symmetries for
2618: Output Parameters:
2619: + size - the number of dofs for points in the `stratum` of the label
2620: . minOrient - the smallest orientation for a point in this `stratum`
2621: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2622: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2623: - 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
2625: Level: developer
2627: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2628: @*/
2629: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2630: {
2631: PetscSectionSym_Label *sl;
2632: const char *name;
2633: PetscInt i;
2635: PetscFunctionBegin;
2637: sl = (PetscSectionSym_Label *)sym->data;
2638: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2639: for (i = 0; i <= sl->numStrata; i++) {
2640: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2642: if (stratum == value) break;
2643: }
2644: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2645: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2646: if (size) {
2648: *size = sl->sizes[i];
2649: }
2650: if (minOrient) {
2652: *minOrient = sl->minMaxOrients[i][0];
2653: }
2654: if (maxOrient) {
2656: *maxOrient = sl->minMaxOrients[i][1];
2657: }
2658: if (perms) {
2660: *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;
2661: }
2662: if (rots) {
2664: *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;
2665: }
2666: PetscFunctionReturn(PETSC_SUCCESS);
2667: }
2669: /*@C
2670: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2672: Logically
2674: InputParameters:
2675: + sym - the section symmetries
2676: . stratum - the stratum value in the label that we are assigning symmetries for
2677: . size - the number of dofs for points in the `stratum` of the label
2678: . minOrient - the smallest orientation for a point in this `stratum`
2679: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2680: . mode - how `sym` should copy the `perms` and `rots` arrays
2681: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2682: - 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
2684: Level: developer
2686: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2687: @*/
2688: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2689: {
2690: PetscSectionSym_Label *sl;
2691: const char *name;
2692: PetscInt i, j, k;
2694: PetscFunctionBegin;
2696: sl = (PetscSectionSym_Label *)sym->data;
2697: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2698: for (i = 0; i <= sl->numStrata; i++) {
2699: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2701: if (stratum == value) break;
2702: }
2703: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2704: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2705: sl->sizes[i] = size;
2706: sl->modes[i] = mode;
2707: sl->minMaxOrients[i][0] = minOrient;
2708: sl->minMaxOrients[i][1] = maxOrient;
2709: if (mode == PETSC_COPY_VALUES) {
2710: if (perms) {
2711: PetscInt **ownPerms;
2713: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2714: for (j = 0; j < maxOrient - minOrient; j++) {
2715: if (perms[j]) {
2716: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2717: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2718: }
2719: }
2720: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2721: }
2722: if (rots) {
2723: PetscScalar **ownRots;
2725: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2726: for (j = 0; j < maxOrient - minOrient; j++) {
2727: if (rots[j]) {
2728: PetscCall(PetscMalloc1(size, &ownRots[j]));
2729: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2730: }
2731: }
2732: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2733: }
2734: } else {
2735: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2736: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2737: }
2738: PetscFunctionReturn(PETSC_SUCCESS);
2739: }
2741: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2742: {
2743: PetscInt i, j, numStrata;
2744: PetscSectionSym_Label *sl;
2745: DMLabel label;
2747: PetscFunctionBegin;
2748: sl = (PetscSectionSym_Label *)sym->data;
2749: numStrata = sl->numStrata;
2750: label = sl->label;
2751: for (i = 0; i < numPoints; i++) {
2752: PetscInt point = points[2 * i];
2753: PetscInt ornt = points[2 * i + 1];
2755: for (j = 0; j < numStrata; j++) {
2756: if (label->validIS[j]) {
2757: PetscInt k;
2759: PetscCall(ISLocate(label->points[j], point, &k));
2760: if (k >= 0) break;
2761: } else {
2762: PetscBool has;
2764: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2765: if (has) break;
2766: }
2767: }
2768: 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],
2769: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2770: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2771: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2772: }
2773: PetscFunctionReturn(PETSC_SUCCESS);
2774: }
2776: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2777: {
2778: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2779: IS valIS;
2780: const PetscInt *values;
2781: PetscInt Nv, v;
2783: PetscFunctionBegin;
2784: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2785: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2786: PetscCall(ISGetIndices(valIS, &values));
2787: for (v = 0; v < Nv; ++v) {
2788: const PetscInt val = values[v];
2789: PetscInt size, minOrient, maxOrient;
2790: const PetscInt **perms;
2791: const PetscScalar **rots;
2793: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2794: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
2795: }
2796: PetscCall(ISDestroy(&valIS));
2797: PetscFunctionReturn(PETSC_SUCCESS);
2798: }
2800: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2801: {
2802: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2803: DMLabel dlabel;
2805: PetscFunctionBegin;
2806: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
2807: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
2808: PetscCall(DMLabelDestroy(&dlabel));
2809: PetscCall(PetscSectionSymCopy(sym, *dsym));
2810: PetscFunctionReturn(PETSC_SUCCESS);
2811: }
2813: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2814: {
2815: PetscSectionSym_Label *sl;
2817: PetscFunctionBegin;
2818: PetscCall(PetscNew(&sl));
2819: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2820: sym->ops->distribute = PetscSectionSymDistribute_Label;
2821: sym->ops->copy = PetscSectionSymCopy_Label;
2822: sym->ops->view = PetscSectionSymView_Label;
2823: sym->ops->destroy = PetscSectionSymDestroy_Label;
2824: sym->data = (void *)sl;
2825: PetscFunctionReturn(PETSC_SUCCESS);
2826: }
2828: /*@
2829: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2831: Collective
2833: Input Parameters:
2834: + comm - the MPI communicator for the new symmetry
2835: - label - the label defining the strata
2837: Output Parameter:
2838: . sym - the section symmetries
2840: Level: developer
2842: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2843: @*/
2844: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2845: {
2846: PetscFunctionBegin;
2847: PetscCall(DMInitializePackage());
2848: PetscCall(PetscSectionSymCreate(comm, sym));
2849: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
2850: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
2851: PetscFunctionReturn(PETSC_SUCCESS);
2852: }